前言
我们在前两期的推文：muscat：专注于多样本多分组的单细胞差异分析；muscat代码实操（一）：模拟复杂的 scRNA-seq 数据以及评估不同的差异表达分析方法的性能中，分别介绍了muscat的功能框架以及如何使用muscat模拟复杂的scRNA-seq数据和评估不同的差异表达分析方法的性能。在本期的推文中，我们将深入探讨muscat包的核心功能——差异状态分析（DS分析）。DS分析的目标是揭示细胞在不同条件或时间点下的功能状态及其行为的变化。
为了让大家更系统地理解，我们将按照以下五个步骤进行介绍：
- 单细胞数据的质控
- 将数据调整为muscat包所需的格式
- DS分析
- 结果的筛选和过滤
- 结果可视化
希望通过这样的步骤，大家能够更加清晰地掌握DS分析的全过程，并能够在实际研究中进行有效应用。接下来，让我们来一起学习一下这项分析吧。
代码流程
0. 加载R包和数据
数据来源：GSE96583；包含8名狼疮患者在IFN-β治疗前后获得的scRNA-seq PBCM数据，数据集包含大约35000个基因和29000个细胞。
library(dplyr)
library(ggplot2)
library(limma)
library(muscat)
library(purrr)
library(ExperimentHub)
eh <- ExperimentHub()
query(eh, "Kang")
(sce <- eh[["EH2259"]])
1. 数据质控和格式准备
数据质控
# remove undetected genes
sce <- sce[rowSums(counts(sce) > 0) > 0, ]
dim(sce)
## [1] 18890 29065
# calculate per-cell quality control (QC) metrics
library(scater)
qc <- perCellQCMetrics(sce)
# remove cells with few or many detected genes
ol <- isOutlier(metric = qc$detected, nmads = 2, log = TRUE)
sce <- sce[, !ol]
dim(sce)
## [1] 18890 26820
# remove lowly expressed genes
sce <- sce[rowSums(counts(sce) > 1) >= 10, ]
dim(sce)
## [1] 7118 26820
# compute sum-factors & normalize
sce <- computeLibraryFactors(sce)
sce <- logNormCounts(sce)
library(sctransform)
assays(sce)$vstresiduals <- vst(counts(sce), verbosity = FALSE)$y
接下来，我们来看一下muscat要求输入的数据格式：
- sample_id: 唯一的样本标识符，例如PeterPan_ref1、Nautilus_trt3等。
- cluster_id: 子群体（或称为簇）的分配，例如T cells、monocytes等。
- group_id: 实验组/条件，例如control/treatment、healthy/diseased等。
sce$id <- paste0(sce$stim, sce$ind)
(sce <- prepSCE(sce,
kid = "cell", # subpopulation assignments
gid = "stim", # group IDs (ctrl/stim)
sid = "id", # sample IDs (ctrl/stim.1234)
drop = TRUE)) # drop all other colData columns
存储簇和样本ID以及它们的数量
nk <- length(kids <- levels(sce$cluster_id))
ns <- length(sids <- levels(sce$sample_id))
names(kids) <- kids
names(sids) <- sids
2. 数据概览
2.1 簇-样本大小
计算每个簇-样本组合中的细胞数量，并将结果转置。
# nb. of cells per cluster-sample
t(table(sce$cluster_id, sce$sample_id))
2.2 降维
# compute UMAP using 1st 20 PCs
sce <- runUMAP(sce, pca = 20)
接下来，定义了一个简单的包装函数.plot_dr()来美化降维图的输出。
.plot_dr <- function(sce, dr, col)
plotReducedDim(sce, dimred = dr, colour_by = col) +
guides(fill = guide_legend(override.aes = list(alpha = 1, size = 3))) +
theme_minimal() + theme(aspect.ratio = 1)
# downsample to max. 100 cells per cluster
cs_by_k <- split(colnames(sce), sce$cluster_id)
cs100 <- unlist(sapply(cs_by_k, function(u)
sample(u, min(length(u), 100))))
# plot t-SNE & UMAP colored by cluster & group ID
for (dr in c("TSNE", "UMAP"))
for (col in c("cluster_id", "group_id"))
.plot_dr(sce[, cs100], dr, col)
对每个簇进行下采样，最多选择100个细胞。然后，它使用前面定义的.plot_dr()函数绘制t-SNE和UMAP图，这些图由簇ID和组ID分别着色。
image.png
3. 差异状态分析
作者主要考虑两种方法进行差异状态（DS）分析：
- 直接作用于细胞水平的混合模型进行DS分析；
- 基于聚合的方法对伪批量数据进行DS分析。
3.1 聚合单细胞数据为伪批量数据
为了利用现有的稳健的批量RNA-seq DE框架，如edgeR、DESeq2和limma，作者首先聚合每个样本（在每个簇中）的测量值以获得伪批量数据。
pb <- aggregateData(sce,
assay = "counts", fun = "sum",
by = c("cluster_id", "sample_id"))
#使用aggregateData函数将数据聚合为伪批量数据。
# one sheet per subpopulation
assayNames(pb)
# pseudobulks for 1st subpopulation
t(head(assay(pb)))
3.2 伪批量级别的MDS图
在进行任何正式的测试之前，我们可以计算聚合信号的多维缩放（MDS）图来探索样本之间的整体相似性。
(pb_mds <- pbMDS(pb))
#使用pbMDS函数计算伪批量数据的MDS图。
image.png
对上图进行美化
# use very distinctive shaping of groups & change cluster colors
pb_mds <- pb_mds +
scale_shape_manual(values = c(17, 4)) +
scale_color_manual(values = RColorBrewer::brewer.pal(8, "Set2"))
# change point size & alpha
pb_mds$layers[[1]]$aes_params$size <- 5
pb_mds$layers[[1]]$aes_params$alpha <- 0.6
pb_mds
image.png
3.3 样本级分析：伪批量方法
# run DS analysis
res <- pbDS(pb, verbose = FALSE)
# access results table for 1st comparison
tbl <- res$table[[1]]
# one data.frame per cluster
names(tbl)
# view results for 1st cluster
k1 <- tbl[[1]]
head(format(k1[, -ncol(k1)], digits = 2))
考虑到实验设计的复杂性（例如，当存在两个以上的组时），可能需要明确指定感兴趣的比较。在这里作者提供了model.matrix函数为pbDS提供一个捕获实验设计的设计矩阵和limma包中的makeContrasts函数指定我们感兴趣的比较。
# construct design & contrast matrix
ei <- metadata(sce)$experiment_info
mm <- model.matrix(~ 0 + ei$group_id)
dimnames(mm) <- list(ei$sample_id, levels(ei$group_id))
library(limma)
contrast <- makeContrasts("stim-ctrl", levels = mm)
# run DS analysis
pbDS(pb, design = mm, contrast = contrast)
3.4 细胞级分析：混合模型
muscat提供了使用3种MM方法的实现。
# 1st approach
mm <- mmDS(sce, method = "dream",
n_cells = 10, n_samples = 2,
min_counts = 1, min_cells = 20)
# 2nd & 3rd approach
mm <- mmDS(sce, method = "vst", vst = "sctransform")
mm <- mmDS(sce, method = "nbinom")
4. 结果
接下来，我们介绍一下如何查看差异测试结果，以及计算和使用基因的表达频率进行进一步的结果过滤。
4.1 结果过滤
tbl_fil <- lapply(tbl, function(u) {
u <- dplyr::filter(u, p_adj.loc < 0.05, abs(logFC) > 1)
dplyr::arrange(u, p_adj.loc)
})
首先过滤出FDR小于5%和logFC大于1的结果，并按调整后的p值排序。
4.2 计算每个簇的DS基因数量和百分比
n_de <- vapply(tbl_fil, nrow, numeric(1))
p_de <- format(n_de / nrow(sce) * 100, digits = 3)
data.frame("#DS" = n_de, "%DS" = p_de, check.names = FALSE)
4.3 查看每个簇的前两个结果
top2 <- bind_rows(lapply(tbl_fil, top_n, 2, p_adj.loc))
format(top2[, -ncol(top2)], digits = 2)
4.4 计算表达频率
每个基因的表达频率，即每个样品和/或组中表达给定基因的细胞比例。frq <- calcExprFreqs(sce, assay = "counts", th = 0) 使用calcExprFreqs函数计算每个基因在每个样本和/或组中的表达频率。
# one sheet per cluster
assayNames(frq)
## [1] "B cells" "CD14+ Monocytes" "CD4 T cells"
## [4] "CD8 T cells" "Dendritic cells" "FCGR3A+ Monocytes"
## [7] "Megakaryocytes" "NK cells"
# expression frequencies in each
# sample & group; 1st cluster
t(head(assay(frq), 5))
## NOC2L HES4 ISG15 TNFRSF18 TNFRSF4
## ctrl101 0.09734513 0.000000000 0.08849558 0.061946903 0.008849558
## ctrl1015 0.07773109 0.008403361 0.18067227 0.090336134 0.021008403
## ctrl1016 0.08333333 0.013888889 0.20833333 0.048611111 0.000000000
## ctrl1039 0.06666667 0.033333333 0.33333333 0.100000000 0.033333333
## ctrl107 0.09803922 0.000000000 0.07843137 0.058823529 0.019607843
## ctrl1244 0.12686567 0.037313433 0.10447761 0.126865672 0.029850746
## ctrl1256 0.08750000 0.004166667 0.18750000 0.091666667 0.025000000
## ctrl1488 0.06837607 0.004273504 0.11965812 0.076923077 0.012820513
## stim101 0.08333333 0.055555556 0.98611111 0.020833333 0.006944444
## stim1015 0.08963585 0.072829132 0.99719888 0.044817927 0.005602241
## stim1016 0.03875969 0.054263566 1.00000000 0.007751938 0.023255814
## stim1039 0.05128205 0.051282051 1.00000000 0.025641026 0.025641026
## stim107 0.05357143 0.035714286 1.00000000 0.017857143 0.000000000
## stim1244 0.11702128 0.095744681 0.98936170 0.053191489 0.021276596
## stim1256 0.06635071 0.037914692 0.99526066 0.037914692 0.004739336
## stim1488 0.05653710 0.014134276 0.99646643 0.031802120 0.003533569
## ctrl 0.08509142 0.009845288 0.15963432 0.084388186 0.018284107
## stim 0.07235339 0.050266565 0.99543031 0.033511043 0.008377761
4.5 基于表达频率过滤结果
gids <- levels(sce$group_id)
frq10 <- vapply(as.list(assays(frq)),
function(u) apply(u[, gids] > 0.1, 1, any),
logical(nrow(sce)))
t(head(frq10))
## NOC2L HES4 ISG15 TNFRSF18 TNFRSF4 SDF4
## B cells FALSE FALSE TRUE FALSE FALSE FALSE
## CD14+ Monocytes FALSE TRUE TRUE FALSE FALSE TRUE
## CD4 T cells FALSE FALSE TRUE FALSE FALSE FALSE
## CD8 T cells FALSE FALSE TRUE FALSE FALSE FALSE
## Dendritic cells FALSE TRUE TRUE FALSE TRUE TRUE
## FCGR3A+ Monocytes FALSE TRUE TRUE FALSE FALSE FALSE
## Megakaryocytes FALSE FALSE TRUE FALSE FALSE FALSE
## NK cells FALSE FALSE TRUE TRUE FALSE TRUE
tbl_fil2 <- lapply(kids, function(k)
dplyr::filter(tbl_fil[[k]],
gene %in% names(which(frq10[, k]))))
#首先确定哪些基因在至少一个组中的平均表达频率超过10%。随后过滤出这些基因的DS结果。
4.6 再次计算每个簇的DS基因数量和百分比
n_de <- vapply(tbl_fil2, nrow, numeric(1))
p_de <- format(n_de / nrow(sce) * 100, digits = 3)
data.frame("#DS" = n_de, "%DS" = p_de, check.names = FALSE)
## #DS %DS
## B cells 226 3.175
## CD14+ Monocytes 622 8.738
## CD4 T cells 154 2.164
## CD8 T cells 94 1.321
## Dendritic cells 117 1.644
## FCGR3A+ Monocytes 382 5.367
## Megakaryocytes 28 0.393
## NK cells 160 2.248
再次计算每个簇中DS基因的数量和百分比，但这次是基于表达频率过滤后的结果。
4.7 结果格式化
当进行多个对比或系数测试时，runDS返回的结果可能会变得非常复杂，不便于探索或导出。resDS函数可以用来格式化这些结果，并提供了两种不同的格式化模式：bind = "row" 或 bind = "col"。
resDS(sce, res, bind = "row", frq = frq)
#所有比较的结果都会垂直合并（类似于do.call("rbind", ...)）成一个整洁格式的表格，其中contrast/coef列指定了比较。这里还附加了预先计算的表达频率。
resDS(sce, res, bind = "col", cpm = TRUE)
#结果会水平合并成一个单一的宽表，其中给定基因和簇的所有结果都保留在一行中。相应的对比或系数的标识符会追加到列名中。这种格式在想要查看特定基因在多种处理中的行为时很有用，但在包含许多比较时会变得混乱。
4.8 直接计算表达频率
如果表达频率尚未使用calcExprFreqs预先计算，可以通过指定frq = TRUE直接将其添加到结果表中：
resDS(sce, res, frq = TRUE)
5. 可视化结果
差异性分析的目的是在不同条件下识别特定细胞群体的状态（或表达）的变化。在这种情境下，关键的问题是：哪些基因只在一个（或很少）细胞群体中显示差异表达？有多少差异表达的基因在细胞群体之间是共享的？也就是不同细胞群体之间的差异性发现的一致性是怎样的？
5.1 UpSet图
为了了解不同细胞群体之间关于差异表达基因的一致性，作者使用UpSet图（UpSet图是一种可视化集合交集的方法）可视化了在某些细胞群体中共享或独特的差异表达基因的数量。
library(UpSetR)
de_gs_by_k <- map(tbl_fil, "gene")
upset(fromList(de_gs_by_k))
image.png
5.2 DR colored by expression
一组由基因表达着色的t-SNE图，这些基因是所有细胞群体中的前8个DS基因。
top8 <- bind_rows(tbl_fil) %>%
top_n(8, dplyr::desc(p_adj.loc)) %>%
pull("gene")
#这里使用dplyr包的函数来提取所有细胞群体中的前8个DS基因。
ps <- lapply(top8, function(g)
.plot_dr(sce[, cs100], "TSNE", g) +
ggtitle(g) + theme(legend.position = "none"))
#对于top8中的每个基因，都会生成一个由该基因的表达着色的t-SNE图。
plot_grid(plotlist = ps, ncol = 4, align = "vh")#排列并显示图形
image.png
5.3 细胞水平的小提琴图
plotExpression(sce[, sce$cluster_id == "B cells"],
features = tbl_fil$`B cells`$gene[seq_len(6)],
x = "sample_id", colour_by = "group_id", ncol = 3) +
guides(fill = guide_legend(override.aes = list(size = 5, alpha = 1))) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
image.png
5.4 Pseudobulk 热图
pbHeatmap(sce, res, top_n = 5)#展示每个细胞群体的前5个DS基因。
image.png
为单个细胞群体生成热图:
pbHeatmap(sce, res, k = "B cells")#这里展示了B细胞群体的前20个DS基因。
image.png
为所有细胞群体生成单个基因的热图:
pbHeatmap(sce, res, g = "ISG20")
image.png
小结
截止到目前，有关muscat的使用方法就都介绍完了。muscat在样本级别上对scRNA-seq进行差异分析，在一定程度上克服的scRNA-seq数据高噪声和稀疏性导致的差异基因表达分析的假阳性，有利于我们在研究中得到更为可靠的结论。在实际使用muscat的过程中，我们不难发现，该R包的设计非常全面和细致。它涵盖了从数据的质量控制到结果的可视化的整个流程，为科研工作者提供了极大的便利。感兴趣的小伙伴赶紧用起来吧。
好啦，本期的分享到这里就结束了，我们下期再会~~
