Pseudobulk differential expression analysis

To create “pseudobulk” RNA libraries for differential expression analysis, all unsorted cells in a given sample were computationally “pooled” according to their automated cluster assignment by summing all counts for a given gene in a given cell-type in a given sample. In addition, this procedure was applied to specifically gated populations of interest (circulating T follicular helper (cTfh), Plasmablast, CD71 high B cells) from the sorted cells as these were not present at high enough quantities in the unsorted cells. Pseudobulk libraries made up by few cells and therefore likely not modeled properly by bulk differential expression methods were removed from analysis for each cell-type separately to remove samples that contained fewer than 8 cells and less than 40000 unique molecular identifier counts detected after pooling. Lowly expressed genes were removed for each cell type individually using the filterByExpr function from edgeR (^{McCarthy et al., 2012}). Differentially expressed genes were identified using the limma voom (^{Law et al., 2014}) workflow which models the log of the counts per million (cpm) of each gene. Scaling factors for library size normalization were calculated with the calcNormFactors function with method = “RLE.” In models utilizing repeated samples from the same subject, subject random effects were incorporated using the duplicateCorrelation function in limma. Models controlled for age, batch, and days since symptom onset, with the exception of models comparing to HC, where days since symptom onset was not controlled for as HC do not have a time since onset. Age was modeled as a continuous variable representing the age in years of the subject. The batch was modeled as factor variable representing the experimental batch (3 in total), and days since symptom onset was modeled as a continuous variable representing the number of days since self-reported symptom onset.