GitHub - jasonratcliff/aggregateBioVar: Aggregation of Multi-subject scRNA-seq Data (original) (raw)
aggregateBioVar
Single cell RNA sequencing (scRNA-seq) studies allow gene expression quantification at the level of individual cells. These studies introduce multiple layers of biological complexity, including variations in gene expression between cell states within a sample (e.g., T cells versus macrophages), between samples within a population (e.g., biological or technical replicates), and between populations (e.g., healthy versus diseased individuals). Many early scRNA-seq studies involved analysis of gene expression within cells from a single sample. For single cell RNA-seq data collected from more than one subject, aggregateBioVarprovides tools to summarize summarize single cell gene expression profiles at the level of samples (i.e., subjects) or populations. Given an inputSingleCellExperimentobject (Amezquita et al. 2020) with pre-defined cell states, aggregateBioVar() stratifies data as a list ofSummarizedExperimentobjects (Morgan et al. 2020). For each cell type, gene counts are aggregated by subject into agene-by-subject count matrix, and column metadata are summarized to retain inter-subject variation for downstream analysis with bulk RNA-seq tools.
Installation
Install the development version of aggregateBioVar fromGitHub with:
install.packages("devtools")
devtools::install_github("jasonratcliff/aggregateBioVar", build_vignettes=TRUE)
Multi-subject scRNA-seq
library(aggregateBioVar)
Bioconductor Packages
library(SummarizedExperiment, quietly = TRUE) library(SingleCellExperiment, quietly = TRUE) library(DESeq2, quietly = TRUE)
Data analysis and visualization
library(dplyr, quietly = TRUE) library(magrittr, quietly = TRUE) library(ggplot2, quietly = TRUE) library(ggtext, quietly = TRUE)
To illustrate the utility of biological replication for scRNA-seq sequencing experiments, consider a SingleCellExperiment object with scRNA-seq data from 7 subjects (SCF1, SCF2, SCF3, SWT1, SWT2, SWT3, SWT4) in the context of a cystic fibrosis phenotype. Samples were collected from small airway epithelium of newborn Sus scrofa with genotypes from wild type (CFTR+/+, n=4) and _CFTR_-knockout (CFTR-/-, n=3) individuals. Note the dimensions of this object, with 1339 genes from 2722 cells:
small_airway #> class: SingleCellExperiment #> dim: 1339 2722 #> metadata(0): #> assays(1): counts #> rownames(1339): MPC1 PRKN ... OTOP1 UNC80 #> rowData names(0): #> colnames(2722): SWT1_AAAGAACAGACATAAC SWT1_AAAGGATTCTCCGAAA ... #> SCF3_TTTGGAGGTAAGGCTG SCF3_TTTGGTTCAATAGTAG #> colData names(6): orig.ident nCount_RNA ... Region celltype #> reducedDimNames(0): #> altExpNames(0):
The primary function aggregateBioVar() takes a SingleCellExperimentobject with column metadata variables indicating subject identity (e.g., biological sample; subjectVar) and assigned cell type (cellVar). The column metadata of a SingleCellExperiment object can be obtained by SummarizedExperiment::colData(). Here, the metadata variable orig.ident indicates the biological sample identifier andcelltype the inferred cell type.
Perform aggregation of counts and metadata by subject and cell type.
aggregate_counts <- aggregateBioVar( scExp = small_airway, subjectVar = "orig.ident", cellVar = "celltype" ) #> Coercing metadata variable to character: celltype
Each element of the returned list contains a SummarizedExperimentobject with aggregated counts from cells in the assigned cell type (indicated by cellVar).
aggregate_counts
#> $AllCells
#> class: SummarizedExperiment
#> dim: 1339 7
#> metadata(0):
#> assays(1): counts
#> rownames(1339): MPC1 PRKN ... OTOP1 UNC80
#> rowData names(0):
#> colnames(7): SWT1 SWT2 ... SWT4 SCF3
#> colData names(3): orig.ident Genotype Region
#>
#> $Immune cell
#> class: SummarizedExperiment
#> dim: 1339 7
#> metadata(0):
#> assays(1): counts
#> rownames(1339): MPC1 PRKN ... OTOP1 UNC80
#> rowData names(0):
#> colnames(7): SWT1 SWT2 ... SWT4 SCF3
#> colData names(4): orig.ident Genotype Region celltype
#>
#> $Secretory cell
#> class: SummarizedExperiment
#> dim: 1339 7
#> metadata(0):
#> assays(1): counts
#> rownames(1339): MPC1 PRKN ... OTOP1 UNC80
#> rowData names(0):
#> colnames(7): SWT1 SWT2 ... SWT4 SCF3
#> colData names(4): orig.ident Genotype Region celltype
#>
#> $Endothelial cell
#> class: SummarizedExperiment
#> dim: 1339 7
#> metadata(0):
#> assays(1): counts
#> rownames(1339): MPC1 PRKN ... OTOP1 UNC80
#> rowData names(0):
#> colnames(7): SWT1 SWT2 ... SWT4 SCF3
#> colData names(4): orig.ident Genotype Region celltype
For each cell type subset, within-subject gene counts are aggregated and column metadata are summarized to exclude variables with _intercellular_variation. This effectively retains subject metadata and can be used for downstream analysis with bulk RNA-seq tools. After aggregation, the number of columns in the SingleCellExperiment object matches the number of unique values in the subject metadata variable indicated bysubjectVar.
assay(aggregate_counts$Immune cell, "counts")
#> DataFrame with 1339 rows and 7 columns
#> SWT1 SWT2 SWT3 SCF1 SCF2 SWT4 SCF3
#>
#> MPC1 257 209 56 315 33 93 42
#> PRKN 1 0 0 0 0 2 0
#> SLC22A3 0 0 0 1 0 0 0
#> CNKSR3 3 4 5 8 1 1 6
#> UTRN 166 155 32 233 27 116 20
#> ... ... ... ... ... ... ... ...
#> CNNM1 0 0 0 0 0 0 0
#> GABRQ 0 0 0 0 0 0 0
#> GIF 0 0 0 0 0 0 0
#> OTOP1 0 0 0 0 0 0 0
#> UNC80 0 0 0 0 0 0 0
colData(aggregate_counts$Immune cell)
#> DataFrame with 7 rows and 4 columns
#> orig.ident Genotype Region celltype
#>
#> SWT1 SWT1 WT Small Immune cell
#> SWT2 SWT2 WT Small Immune cell
#> SWT3 SWT3 WT Small Immune cell
#> SCF1 SCF1 CFTRKO Small Immune cell
#> SCF2 SCF2 CFTRKO Small Immune cell
#> SWT4 SWT4 WT Small Immune cell
#> SCF3 SCF3 CFTRKO Small Immune cell
Differential Gene Expression
The aggregate gene-by-subject matrix and subject metadata can be used as inputs for bulk RNA-seq tools to investigate gene expression. Here, an example is provided usingDESeq2(Love, Huber, and Anders 2014). A DESeqDataSet can be constructed from the aggregate gene-by-subject count matrix and summarized column metadata.
subj_dds_dataset <-
DESeqDataSetFromMatrix(
countData = assay(aggregate_counts$Secretory cell, "counts"),
colData = colData(aggregate_counts$Secretory cell),
design = ~ Genotype
)
#> converting counts to integer mode
#> Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
#> design formula are characters, converting to factors
subj_dds <- DESeq(subj_dds_dataset) #> estimating size factors #> estimating dispersions #> gene-wise dispersion estimates #> mean-dispersion relationship #> final dispersion estimates #> fitting model and testing
subj_dds_results <- results(subj_dds, contrast = c("Genotype", "WT", "CFTRKO"))
Add negative log10 adjusted P-values, then plot against log2 fold change. Genes with adjusted P-values < 0.05 and fold-change absolute values > 1.0 are highlighted in red and labeled by feature.
subj_dds_transf <- as.data.frame(subj_dds_results) %>% bind_cols(feature = rownames(subj_dds_results)) %>% mutate(log_padj = - log(.data$padj, base = 10))
ggplot(data = subj_dds_transf) + geom_point(aes(x = log2FoldChange, y = log_padj), na.rm = TRUE) + geom_point( data = filter( .data = subj_dds_transf, abs(.data$log2FoldChange) > 1, .data$padj < 0.05 ), aes(x = log2FoldChange, y = log_padj), color = "red" ) + geom_label( data = filter( .data = subj_dds_transf, abs(.data$log2FoldChange) > 1, .data$padj < 0.05 ), aes(x = log2FoldChange, y = log_padj + 0.4, label = feature) ) + theme_classic() + labs( x = "log2 (fold change)", y = "-log10 (padj)" ) + theme( axis.title.x = element_markdown(), axis.title.y = element_markdown())
Vignettes
For a detailed workflow and description of package components, see the package vignette:
vignette("multi-subject-scRNA-seq", package = "aggregateBioVar")
References
Amezquita, Robert, Aaron Lun, Etienne Becht, Vince Carey, Lindsay Carpp, Ludwig Geistlinger, Federico Marini, et al. 2020. “Orchestrating Single-Cell Analysis with Bioconductor.” Nature Methods 17: 137–45.https://www.nature.com/articles/s41592-019-0654-x.
Bache, Stefan Milton, and Hadley Wickham. 2014. magrittr: A Forward-Pipe Operator for R.https://CRAN.R-project.org/package=magrittr.
Love, Michael I., Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-seq Data with DESeq2.”Genome Biology 15 (12): 550.https://doi.org/10.1186/s13059-014-0550-8.
Morgan, Martin, Valerie Obenchain, Jim Hester, and Hervé Pagès. 2020.SummarizedExperiment: SummarizedExperiment Container.
Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2020.dplyr: A Grammar of Data Manipulation.https://CRAN.R-project.org/package=dplyr.
Wilke, Claus O. 2019. cowplot: Streamlined Plot Theme and Plot Annotations for ’ggplot2’.https://CRAN.R-project.org/package=cowplot.
———. 2020. ggtext: Improved Text Rendering Support for ’ggplot2’.https://CRAN.R-project.org/package=ggtext.