README.md
In jasonratcliff/aggregateBioVar: Differential Gene Expression Analysis for Multi-subject scRNA-seq
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, aggregateBioVar
provides tools to summarize summarize single cell gene expression
profiles at the level of samples (i.e., subjects) or populations.
Given an input
SingleCellExperiment
object (Amezquita et al. 2020) with
pre-defined cell states, aggregateBioVar() stratifies data as a list
of
SummarizedExperiment
objects (Morgan et al. 2020). For each
cell type, gene counts are aggregated by subject into a
gene-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 from
GitHub 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 SingleCellExperiment
object 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 and
celltype 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 SummarizedExperiment
object 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 by
subjectVar.
assay(aggregate_counts$`Immune cell`, "counts")
#> DataFrame with 1339 rows and 7 columns
#> SWT1 SWT2 SWT3 SCF1 SCF2 SWT4 SCF3
#> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
#> 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
#> <character> <character> <character> <factor>
#> 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 using
DESeq2
(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 = "log<sub>2</sub> (fold change)",
y = "-log<sub>10</sub> (p<sub>adj</sub>)"
) +
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.
.
Bache, Stefan Milton, and Hadley Wickham. 2014. *magrittr: A
Forward-Pipe Operator for R*.
.
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.
.
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*.
.
Wilke, Claus O. 2019. *cowplot: Streamlined Plot Theme and Plot
Annotations for ’ggplot2’*.
.
———. 2020. *ggtext: Improved Text Rendering Support for ’ggplot2’*.
.
jasonratcliff/aggregateBioVar documentation built on Sept. 29, 2020, 11:53 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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, aggregateBioVar
provides tools to summarize summarize single cell gene expression
profiles at the level of samples (i.e., subjects) or populations.
Given an input
SingleCellExperiment
object (Amezquita et al. 2020) with
pre-defined cell states, aggregateBioVar() stratifies data as a list
of
SummarizedExperiment
objects (Morgan et al. 2020). For each
cell type, gene counts are aggregated by subject into a
gene-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 from
GitHub 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 SingleCellExperiment
object 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 and
celltype 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 SummarizedExperiment
object 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 by
subjectVar.
assay(aggregate_counts$`Immune cell`, "counts")
#> DataFrame with 1339 rows and 7 columns
#> SWT1 SWT2 SWT3 SCF1 SCF2 SWT4 SCF3
#> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
#> 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
#> <character> <character> <character> <factor>
#> 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 using
DESeq2
(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 = "log<sub>2</sub> (fold change)",
y = "-log<sub>10</sub> (p<sub>adj</sub>)"
) +
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.
.
Bache, Stefan Milton, and Hadley Wickham. 2014. *magrittr: A
Forward-Pipe Operator for R*.
.
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.
.
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*.
.
Wilke, Claus O. 2019. *cowplot: Streamlined Plot Theme and Plot
Annotations for ’ggplot2’*.
.
———. 2020. *ggtext: Improved Text Rendering Support for ’ggplot2’*.
.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.