To analyze our TEA-seq dataset, 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.
All code used for the analysis and figures in our manuscript, along with additional exploratory analyses, are available on Github:
aifimmunology/repro-vrd-tea-seq
TEA-seq data processing
After sequencing, scRNA-seq and scATAC-seq libraries were aligned and transcripts were counted using CellRanger Multi to generate transcript count matrices and unique fragment alignments for each cell barcode. ADT and HTO count matrices were generated using BarCounter v1.0 ( Swanson, E., et al., 2022 ). RNA, HTO, and ADT count matrices were combined and demultiplexed for each treatment condition using HTO count distributions with BarMixer. Additional scripts utilizing BedTools ( Quinlan, A., et al., 2014 ), awk, and ArchR ( Granja, J.M., et al., 2021 ) were used to demultiplex scATAC-seq fragments for each sample and perform initial quality control analysis.
View notebook: 03-R_qc_filtering.ipynb
Cell barcodes were filtered to select cells with the following characteristics: > 500 scRNA-seq genes per cell and > 1,000 unique scATAC-seq fragments per cell. Initial analysis of the scRNA-seq data was performed using Seurat to cluster cells and identify non-T cell contamination and doublets. In total 181,518 cells were considered, 37,690 cells (20.8%) were excluded based on QC criteria and clustering, and 143,828 cells (79.2%) were retained for downstream analysis.
Cell-type labeling
View notebook: 04-R_t_cell_class_filtering.ipynb
Normalized ADT count distributions for CD4 and CD8a antibodies were inspected from each sample to identify cutoffs for CD4-positive and CD8-positive populations. Cells were then separated into CD4 and CD8 T cell classes for downstream labeling steps.
View notebook: 01-R_cd4_type_labeling.ipynb
View notebook: 02-R_cd8_type_labeling.ipynb
For each class, scRNA-seq data was analyzed using Seurat ( Hao, Y., et al., 2022 ): counts were transformed using SCTransform, and label transfer was performed using the PBMC cell type reference (from Azimuth; Hao, Y., et al., 2020 ). To constrain label transfer, reference cells were filtered for those with a top-level cell type label of either CD4 T or CD8 T. Labels were transferred using Seurat's FindTransferAnchors and MapQuery functions for CD4 and CD8 cells and respective references. After label transfer, all CD4 or CD8 cells were integrated across samples using Seurat functions SelectIntegrationFeatures, PrepSCTIntegration, FindIntegrationAnchors, and IntegrateData. Dimensionality reduction was then performed using RunPCA and RunUMAP, and Leiden clustering was performed using FindNeighbors and FindClusters. Marker gene expression and transfer-based labels (above) were used to assign cell types to each of the CD4 and CD8 Leiden clusters.
MAST DEG analysis
View notebook: 01-R_treat-type-downsample_MAST.ipynb
For each cell type, we downsampled cells across all treatment conditions based on the lowest number of cells in any condition. In order to be considered for differential expression testing, a gene must have been detected in at least 10% of cells in either the drug or DMSO control to be considered for analysis. Then, using the downsampled cells and genes that passed detection cutoffs, we tested for differential expression using MAST ( Finak, G., et al., 2015 ) using the following formula:
Exp ~ CDR + Covariates + Treatment
which included inclusion of Cellular Detection Rate (CDR) to be able to jointly estimate nuisance factors and treatment effects.
DEG comparisons by hypergeometric tests
View notebook: 02-R_deg_result_overlaps.ipynb
To account for differences in the number of cells in different conditions, we compared the top 500 up- and down-regulated DEGs. The number of overlapping genes between each set of DEGs were counted, and we identified the number of common genes that were tested in both conditions. We then used the R function phyper to calculate the probability of overlap between gene sets using the parameters q = N overlapping genes - 1, m = 500 (number of genes in Set A), n = number of common genes not in Set A, k = 500 (number of genes in Set B), and lower.tail = FALSE. After computing all tests, P-values were adjusted for multiple hypothesis testing using the p.adjust function in R with the parameter method = "BH" for the Benjamini and Hochberg method .
Pathway enrichment analysis
View notebook: 01-R_mast-deg_hallmark_gsea.ipynb
View notebook: 02-R_mast-deg_reactome_gsea.ipynb
View notebook: 03-R_bortezomib_mast-deg_hallmark_no-PSM.ipynb
We used the fgsea package ( Korotkevich, G., et al., 2021 ) to perform Gene Set Enrichment Analysis (GSEA) for pathway sets obtained from MSigDB (Hallmark Pathways; Liberzon, A., et al., 2015 ) and Reactome (Reactome Pathways; Gillespie, M., et al., 2021 ). Prior to use of fgsea, DEG results from MAST were ranked based on nominal P-values and direction of effect, such that the most significantly up-regulated genes were first and most significantly down-regulated genes were last in order. Enrichment tests for Hallmark Pathways and Reactome Pathways were performed separately.
Differential epitope detection analysis
View notebook: 01-R_DEP_lm_by_treatment-type.ipynb
Differential epitope analysis was performed for each cell type, drug, and time point relative to the DMSO control using normalized ADT count data (Seurat NormalizeData, normalization.method = "CLR", margin = 2). As with DEG analysis, cells were downsampled evenly across timepoints to the minimum cell number by drug and cell type. Epitope expression was modeled linearly as exp ~ treatment with DMSO as the reference treatment level. Treatment coefficient P-values were adjusted using the Benjamini & Hochberg correction.
DAP analysis
View notebook: 02-R_peaks_and_annotation.ipynb
Cells were distributed into separate ArchR projects according to cell type. Peaks were identified using using ArchR's addGroupCoverages and addReproduciblePeakSet function with parameter groupBy = "Sample". Insertion site counts were tabulated using ArchR's addPeakMatrix function. Treatment-specific ArchR projects were generated using only cells from each treatment condition at each time point and the relevant timepoint-matched DMSO control cells, then DAPs were identified using ArchR's getMarkerFeatures function.
Motif enrichment analysis
View notebook: 02-R_peaks_and_annotation.ipynb
For each cell type, peaks were annotated with cisBP motif locations using ArchR's addMotifAnnotations with parameters motifSet = "cisbp" and the default parameter of version = 2 to use cisBP Version 2 motifs ( Weirauch, M., et al., 2014 ). To identify enriched motifs, DAPs were divided into peaks with increased accessibility (false discovery rate [FDR] < 0.05 and Log2FC > 0) and decreased accessibility (FDR < 0.05 and Log2FC < 0). For DAPs in each direction, hypergeometric tests for enrichment were performed using ArchR's peakAnnoEnrichment function.