Analysis methods

To build our longitudinally sampled PBMC scRNA-seq dataset from 95 healthy subjects, we utilized data and analysis environments within the Human Immune System Explorer (HISE) system to trace data processing, analytical code, and analysis environments from raw FASTQ data to a labeled, high-quality final dataset. A graph representation of the analysis trace for this project is available in the Certificate of Reproducibility.

This analysis consisted of the following major steps:

• Consistent pipeline processing of hashed 10x scRNA-seq data
• Sample selection from HISE data storage
• Labeling of cells using our CellTypist models and doublet scoring
• Separation of cell classes (AIFI_L2 labels) and subclustering to identify cross-class doublets and low-quality clusters
• Clustering and label refinement for high-resolution population subset labels (AIFI_L3 labels)
• Assembly of final datasets
• Differential expression tests

These steps and some additional analyses are available at: 
aifimmunology/sound-life-scrna-analysis

Pipeline processing

After sequencing, Gene Expression and Hash Tag Oligonucleotide libraries from pooled samples in our pipeline batches were demultiplexed and assembled as individual files for each biological sample. See our Multiplexed scRNA-seq pipeline documentation for additional details.

Data selection from HISE storage

View Notebook: 00-R_Sample_Selection.ipynb

To assemble the input data for our dataset, we selected all available samples from all subjects in our SoundLife Young Adult (age 25-35) and SoundLife Older Adult (age 55-65) cohorts. We then filtered to remove a batch of samples with previously identified quality problems (B004), as well as subjects with non-healthy or abnormal disease states recorded at the time of any of their sample collection visits. In cases where aliquots from the same blood draw were processed through our pipeline multiple times, we selected the most recently generated data (highest batch number).

In total, 868 samples were selected for use in our dataset (Young Adult, n = 49 subjects; n =  418 samples; Older Adult, n = 47 subjects, n = 450 samples), consisting of 15,781,886 cells before additional QC filtering.

View Notebook: 01-Python_retrieve_cmv_bmi.ipynb

In addition to processed data and sample metadata identification, we retrieved clinical lab results and clinical subject data to assemble CMV infection status for all subjects, as well as BMI for all adult subjects (Sound Life Cohorts).

Labeling and doublet detection

To add cell subset labels, we utilized CellTypist (v1.6.1, Dominguez Conde, et al.) and the CellTypist models we generated using our Immune Health PBMC Atlas dataset at 3 levels of subset label resolution (AIFI_L1, 9 broad cell classes; AIFI_L2, 29 intermediate cell classes, and AIFI_L3, 71 high-resolution cell classes) by following the approach described in the CellTypist reference documentation at https://celltypist.org.

View Notebook: 03-Python_doublet_detection.ipynb

While our cell hashing pipeline identifies doublets that occur between multiple samples, we are blind to doublets that occur with the same cell droplet. To assist in identifying these remaining doublets, we performed doublet detection using the scrublet package (v0.2.3, Wolock, et al.) as implemented in the scanpy.external module of scanpy (v1.9.6, Alexander Wolf, et al.).

View Notebook: 04-Python_assembly_of_subsets.ipynb

After subset labeling and doublet detection were performed on each sample, we assembled data, labels, scrublet calls, sample metadata, CMV status, and BMI data across 8 sets of samples defined by the cohort, biological sex, and CMV status of subjects to facilitate downstream analyses.

QC filtering

View Notebook: 05-Python_subset_qc_filtering.ipynb

After assembly of all cells across our samples, we filtered our data against a set of QC criteria designed to remove possible doublets (based on scrublet and high gene detection, > 5,000 genes), low-quality cells (based on low gene detection, < 200 genes), and dying or dead cells (based on mitochondrial UMI content, > 10%). Mitochondrial content was assessed using scanpy by identifying mitochondrial genes (prefixed with "MT-"), and using the scanpy function calculate_qc_metrics, as demonstrated in the scanpy tutorials (https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html). QC criteria were applied to the dataset in series, as shown in the table below. In total, 951,864 cells were removed (6.03%), and 14,830,022 cells passed QC filtering (93.97%).

Clustering, class subsetting, and marker-based doublet removal

After QC filtering, all remaining cells were clustered within AIFI_L2 labels using a scanpy workflow that consisted of normalization using scanpy.pp.normalize_total with the parameter target_sum = 1e4; log transformation using scanpy.pp.log1p; highly variable gene detection using scanpy.pp.highly_variable_genes; scaling of highly variable gene data using sc.pp.scale; Dimensionality reduction using scanpy.tl.pca with the parameter svd_solver = 'arpack'; Cohort integration using the Harmony method (harmonypy v0.0.9; Korsunsky, et al.) implemented in scanpy.external.harmony_integrate with the parameter key = 'cohort.cohortGuid'; Nearest neighbor calculation with scanpy.pp.neighbors and parameters n_neighbors = 50, use_rep = 'X_pca_harmony', and n_pcs = 30; UMAP embedding with scanpy.tl.umap; and Leiden clustering (leidenalg v0.10.1; Traag, et al.) using scanpy.tl.leiden. Unless specified, default parameters provided by scanpy v1.9.6 were used. 

For most AIFI_L2 cell classes, the processing and clustering described above was performed using cells from all samples. However, 5 cell classes consisting of > 1 million cells were processed separately as 8 sets of cells based on cohort, sex, and CMV status to make file handling and clustering processes more manageable: CD4 naive T cells, CD4 memory T cells, CD8 memory T cells, CD56dim NK cells, and CD14 Monocytes.

View Notebook: 07-Python_filter_cell_classes.ipynb
View Notebook: 07a-R_filter_review_plots.ipynb

After clustering, clusters were filtered to remove doublets based on marker gene expression. For example, CD4 naive T cell clusters with a high fraction of MS4A1 expression were removed as B cell doublets. For each cluster, the fraction of cells with detected gene expression for each marker was computed, and clusters meeting a set of cell class-specific criteria were removed. These thresholds are shown in the table below.

In addition to these filters based on fraction of gene detection, we included an additional filter for specific subsets based on mean gene expression of HBB where this better distinguished doublet populations that included Erythrocyte contamination:

Finally, we also removed clusters with low gene detection based on a single threshold: If > 30% of cells in a cluster had < 750 genes detected, the cluster was removed as low quality. Due to their generally lower gene detection, this threshold was not applied to Erythrocyte or Platetelet populations.

The table below presents the total number of cells removed for each reason listed above. In total, 800,059 cells (5.39%) were removed by this filtering process, and 14,029,963 cells (94.61%) were retained for final clean up based on high-resolution AIFI_L3 labels.

Clustering and refinement of high-resolution labels

View Notebook: 09-Python_partition_L3_data.ipynb

After removal of doublets at the AIFI_L2 label resolution, we assembled cells from all samples for each subset in the high-resolution AIFI_L3 labeling results. Separately for each subset, we performed the scanpy data processing procedure described above to enable a final review of high resolution labels. Each subset was examined to identify remaining doublet clusters or clusters of cells that appeared to be mislabeled based on marker gene expression in collaboration with immunological domain experts.

During this process, the following additional cells were removed from downstream analysis:

In total, 234,117 cells were removed at this stage of label refinement, and 13,795,846 cells were retained for downstream analyses. In addition to cells that were removed from the dataset, we reassigned the following cells based on gene expression to generate the final label set:

In total, 182,331 cells were reassigned in this final pass of cell type refinement.

Population subset frequency analyses

View Notebook: 15-R_assemble_type_frequencies.ipynb

To analyze changes in subset abundance, we tabulated subset frequency for each sample in our dataset for high resolution (L3) label assignments. We also compute proportions of each subset per sample, which were used for Centered Log Ratio (CLR) transformations. To compare abundance between sets of subjects, we used rank-sum Wilcoxon tests on CLR-transformed values. We compared abundance differences between Age Groups (Young Adult vs Older Adult), CMV status (Negative vs Positive), and Biological Sex (Female vs Male) using the first Flu Vaccine time point (Flu Year 1 Day 0) for each subject. For comparisons between longitudinal time points in the flu vaccine series (Day 7 vs Day 0), we used paired signed-rank Wilcoxon tests.

Pseudobulk data assembly

We generated Pseudobulk data by aggregating results for cells with the same high resolution (L3) label within each sample. This aggregation was performed using multiple methods for downstream analysis: summing UMI counts per gene (used for DESeq2 analysis, below), taking the mean of normalized and log transformed values (for visualization).

Differential gene expression tests

We performed differential gene expression tests using the DESeq2 v1.42.0 (Love, Huber, and Anders, 2014) package in R v4.3.2. Prior to differential tests, genes were filtered based on a minimum of 10% gene detection across all single cells used for the comparison. Summed UMI counts from single cells were used for the comparisons. DESeq2 analyses used default parameters. For comparisons of groups of subjects at baseline, we used pseudobulk data from Flu Year 1 Day 0 samples. The DESeq2 design formula for these comparisons was `Aggregated Counts ~ Age + CMV + Sex`, and then contrasts were made for each factor to identify differential results. Genes were considered differentially expressed if their adjusted p-value was below 0.05 and their absolute fold change was greater than 0.1.
 

References

Domínguez Conde C, Xu C, Jarvis LB, Rainbow DB, Wells SB, Gomes T, et al. Cross-tissue immune cell analysis reveals tissue-specific features in humans. Science. 2022;376: eabl5197.
doi:10.1126/science.abl5197

Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15: 550.
doi:10.1186/s13059-014-0550-8

Traag V, Waltman L, van Eck NJ. From Louvain to Leiden: guaranteeing well-connected communities. arXiv [cs.SI]. 2018. Available: http://arxiv.org/abs/1810.08473
doi:10.48550/arXiv.1810.08473

Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data analysis. Genome Biol. 2018;19: 15.
doi:10.1186/s13059-017-1382-0

Wolock SL, Lopez R, Klein AM. Scrublet: Computational Identification of Cell Doublets in Single-Cell Transcriptomic Data. Cell Syst. 2019;8: 281–291.e9.
doi:10.1016/j.cels.2018.11.005