Integration with dreamlet / SingleCellExperiment

Developed by Gabriel Hoffman

Run on 2024-06-17 13:36:07.828705

Load and process single cell data

Here we perform analysis of PBMCs from 8 individuals stimulated with interferon-β Kang, et al, 2018, Nature Biotech. We perform standard processing with dreamlet to compute pseudobulk before applying crumblr.

Here, single cell RNA-seq data is downloaded from ExperimentHub.

library(dreamlet)
library(muscat)
library(ExperimentHub)
library(zenith)
library(scater)

# Download data, specifying EH2259 for the Kang, et al study
eh <- ExperimentHub()
sce <- eh[["EH2259"]]

sce$ind = as.character(sce$ind)

# only keep singlet cells with sufficient reads
sce <- sce[rowSums(counts(sce) > 0) > 0, ]
sce <- sce[, colData(sce)$multiplets == "singlet"]

# compute QC metrics
qc <- perCellQCMetrics(sce)

# remove cells with few or many detected genes
ol <- isOutlier(metric = qc$detected, nmads = 2, log = TRUE)
sce <- sce[, !ol]

# set variable indicating stimulated (stim) or control (ctrl)
sce$StimStatus <- sce$stim

Aggregate to pseudobulk

Dreamlet creates the pseudobulk dataset:

# Since 'ind' is the individual and 'StimStatus' is the stimulus status,
# create unique identifier for each sample
sce$id <- paste0(sce$StimStatus, sce$ind)

# Create pseudobulk data by specifying cluster_id and sample_id for aggregating cells
pb <- aggregateToPseudoBulk(sce,
  assay = "counts",
  cluster_id = "cell",
  sample_id = "id",
  verbose = FALSE)

Process data

Here we evaluate whether the observed cell proportions change in response to interferon-β.

library(crumblr)

# use dreamlet::cellCounts() to extract data
cellCounts(pb)[1:3,1:3]

##          B cells CD14+ Monocytes CD4 T cells
## ctrl101      101             136         288
## ctrl1015     424             644         819
## ctrl1016     119             315         413

# Apply crumblr transformation
# cobj is an EList object compatable with limma workflow
# cobj$E stores transformed values
# cobj$weights stores precision weights
cobj <- crumblr(cellCounts(pb))

Analysis

Now continue on with the downstream analysis

library(variancePartition)

fit <- dream(cobj, ~ StimStatus + ind, colData(pb))
fit <- eBayes(fit)

topTable(fit, coef = "StimStatusstim", number = Inf)

##                         logFC    AveExpr          t     P.Value  adj.P.Val         B
## CD8 T cells       -0.25085170  0.0857175 -4.0787416 0.002436375 0.01949100 -1.279815
## Dendritic cells    0.37386979 -2.1849234  3.1619195 0.010692544 0.02738587 -2.638507
## CD14+ Monocytes   -0.10525402  1.2698117 -3.1226341 0.011413912 0.02738587 -2.709377
## B cells           -0.10478652  0.5516882 -3.0134349 0.013692935 0.02738587 -2.940542
## CD4 T cells       -0.07840101  2.0201947 -2.2318104 0.050869691 0.08139151 -4.128069
## FCGR3A+ Monocytes  0.07425165 -0.2567492  1.6647681 0.128337022 0.17111603 -4.935304
## NK cells           0.10270672  0.3797777  1.5181860 0.161321761 0.18436773 -5.247806
## Megakaryocytes     0.01377768 -1.8655172  0.1555131 0.879651456 0.87965146 -6.198336

Given the results here, we see that CD8 T cells at others change relative abundance following treatment with iterferon-β.