1 前言
scRNA数据分析具有挑战性的工作一般就是对结果的解读。获得不同亚群的细胞很简单，但如何给每个亚群赋予生物学意义就比较困难，既需要对数据掌握，有需要有生物背景，而生物背景知识更新并不是那么快。我们现在掌握的，可能和真实水平还有很大偏差。
为了减弱自身生物背景的不足，可以用一些先验知识帮助判断。最常见的就是通过去了解参与特定生物学过程的经过验证的基因，利用GO、KEGG数据库；另外还可以将我们的表达谱与之前发布的参考数据集做对比，因为参考数据集中样本或细胞的生物学状态都经过了专家的证实。这两种形式（参考数据集、基因集）接下来都会探讨。
准备数据
下面2.2部分会用到
# 准备数据
library(scRNAseq)
sce.seger <- SegerstolpePancreasData()
# 基因注释
library(AnnotationHub)
edb <- AnnotationHub()[["AH73881"]]
symbols <- rowData(sce.seger)$symbol
ens.id <- mapIds(edb, keys=symbols, keytype="SYMBOL", column="GENEID")
ens.id <- ifelse(is.na(ens.id), symbols, ens.id)
# 去除重复行
keep <- !duplicated(ens.id)
sce.seger <- sce.seger[keep,]
rownames(sce.seger) <- ens.id[keep]
# 样本注释
emtab.meta <- colData(sce.seger)[,c("cell type",
"individual", "single cell well quality")]
colnames(emtab.meta) <- c("CellType", "Donor", "Quality")
colData(sce.seger) <- emtab.meta
sce.seger$CellType <- gsub(" cell", "", sce.seger$CellType)
sce.seger$CellType <- paste0(
toupper(substr(sce.seger$CellType, 1, 1)),
substring(sce.seger$CellType, 2))
# 质控
low.qual <- sce.seger$Quality == "low quality cell"
library(scater)
stats <- perCellQCMetrics(sce.seger)
qc <- quickPerCellQC(stats, percent_subsets="altexps_ERCC_percent",
batch=sce.seger$Donor,
subset=!sce.seger$Donor %in% c("HP1504901", "HP1509101"))
sce.seger <- sce.seger[,!(qc$discard | low.qual)]
# 归一化
# 看到quickCluster和computeSumFactors就知道使用的是去卷积化方法
library(scran)
clusters <- quickCluster(sce.seger)
sce.seger <- computeSumFactors(sce.seger, clusters=clusters)
sce.seger <- logNormCounts(sce.seger)
2 使用参考数据集
非常直接的方法，将我们的表达量与参考数据集的表达量去比对，看看和参考的样本有多少类似
这个牵扯到一个分类的问题，可以辅助以标准的机器学习（比如随机森林、支持向量机），任何之前发表过的RNA-Seq数据（bulk或single cell）中对细胞的划分都可以提供帮助。
这里主要是用了2019年发表的细胞注释R包SingleR (Aran et al. 2019) 。它会计算我们数据中的细胞与参考数据集中的细胞之间Spearman相关性，找到最相关的属性赋予我们的细胞。为了减少数据的噪音，这个R包先鉴定marker基因，然后根据marker基因的表达量去计算相关性。
2.1 使用内置的参考注释数据
SingleR中就包含了一些内置数据集，大部分是bulk RNA-Seq或芯片数据中经过筛选的细胞类型。
下面使用的参考数据集就来自Blueprint和ENCODE计划：(Martens and Stunnenberg 2013; The ENCODE Project Consortium 2012).
# 我们的数据
sce.pbmc
## class: SingleCellExperiment
## dim: 33694 3922
## metadata(1): Samples
## assays(2): counts logcounts
## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
## rowData names(2): ID Symbol
## colnames(3922): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ...
## TTTGTCACAGGTCCAC-1 TTTGTCATCCCAAGAT-1
## colData names(4): Sample Barcode sizeFactor label
## reducedDimNames(3): PCA TSNE UMAP
## altExpNames(0):
# 参考数据集
library(SingleR)
ref <- BlueprintEncodeData()
# 然后把我们的细胞在参考数据集中找对应的细胞类型
# 返回的pred结果是一个数据框，每行是我们自己数据的一个细胞
pred <- SingleR(test=sce.pbmc, ref=ref, labels=ref$label.main)
table(pred$labels)
##
## B-cells CD4+ T-cells CD8+ T-cells DC Eosinophils Erythrocytes
## 525 755 1254 1 1 5
## HSC Monocytes NK cells
## 14 1116 251
# 可视化
plotScoreHeatmap(pred)
横坐标是各种细胞类型，纵坐标是我们数据中的细胞；Score也是经过归一化后的，介于0-1
理论上，我们数据中的每个细胞都会有自己对应的细胞类型，并且在所属的细胞类型中，比其他细胞相关性更高（也就是得分更高、颜色差异更明显）。看到只有B-cells、monocytes表现突出，与其他细胞类型区分明显。而像CD4+、CD8+细胞，颜色相近，区分不明显
注意
当然，会存在一些细胞归属不是很清楚的，好像一个细胞对两种细胞类型都有点沾边，但又不好判断。函数会通过一个阈值将这些划分不清楚的细胞定为”NA”
sum(is.na(pred$pruned.labels))
## [1] 76
plotScoreDistribution(pred)
这个函数会把三种情况一起画出来，并且顺序依次是：匹配上的（黄色 assigned）、根本就不匹配的（灰色 other）、不清不楚的（橙色 pruned）
可以看到大部分都是灰色，没有匹配上的
另外，还可以看cluster对label的对应关系。理想情况是，一个cluster只对应一个label，但现实是：可能好几个clusters属于一个label；还有两个label的clusters很相似（例如CD4+、CD8+）
tab <- table(Assigned=pred$pruned.labels, Cluster=colLabels(sce.pbmc))
library(pheatmap)
pheatmap(log2(tab+10), color=colorRampPalette(c("white", "blue"))(101))
可以看到，使用已知的细胞注释可以简化细胞的生物学解释过程，但同时也会将细胞们限定在已知的框架中。如果发现我们聚类分群的结果与注释后的结果相差很大，不用着急去掉某种结果，留着下游再继续探讨差异来源
2.2 使用自定义的注释数据
比如手上哟一个之前注释过的数据，现在想用其中的注释信息对另一个数据进行注释
例如，之前有一个注释好的人类胰腺数据集 sce.muraro Muraro et al. (2016)
## class: SingleCellExperiment
## dim: 16940 2299
## metadata(0):
## assays(2): counts logcounts
## rownames(16940): ENSG00000268895 ENSG00000121410 ... ENSG00000159840
## ENSG00000074755
## rowData names(2): symbol chr
## colnames(2299): D28-1_1 D28-1_2 ... D30-8_93 D30-8_94
## colData names(4): label donor plate sizeFactor
## reducedDimNames(0):
## altExpNames(1): ERCC
去掉那些不明确的细胞类型
sce.muraro <- sce.muraro[,!is.na(sce.muraro$label) &
sce.muraro$label!="unclear"]
table(sce.muraro$label)
##
## acinar alpha beta delta duct endothelial
## 217 795 442 189 239 18
## epsilon mesenchymal pp
## 3 80 96
现在有一个新的数据集sce.seger Segerstolpe et al. (2016)，它本身也包括细胞注释，在数据准备代码中有提及：在sce.seger$CellType，这个注释是领域专家设定的
# 以这个新数据集作为测试数据，sce.muraro作为参考数据，看看注释结果如何
pred.seger <- SingleR(test=sce.seger, ref=sce.muraro,
labels=sce.muraro$label, de.method="wilcox")
table(pred.seger$labels)
##
## acinar alpha beta delta duct endothelial
## 188 889 279 105 385 17
## epsilon mesenchymal pp
## 5 53 169
用这个结果与测试数据本身带的细胞注释结果进行比较：
tab <- table(pred.seger$pruned.labels, sce.seger$CellType)
library(pheatmap)
pheatmap(log2(tab+10), color=colorRampPalette(c("white", "blue"))(101))
可以看到这两种注释结果比较相近
2.3 使用marker基因集
这部分可以作为2.2的延伸 它基于的想法是：高表达的marker基因与每个细胞的生物学功能相关，因此可以用marker基因作为桥梁，连接测试数据集与参考数据集
基于marker基因的方法会更快捷，同时也不受限于参考数据集是否有专家定义好的细胞类型
参考数据集 sce.zeisel Zeisel et al. (2015)
目的是将数据集的marker基因于对应的细胞类型联系起来
需要注意的是：下面的pairwiseWilcox()代码需要基于R 4.0，需要加载归一化后的的数据，前面的数据处理在之前章节都有提及 数据在：https://share.weiyun.com/fDePnIGC
load('normalized.sce.zeisel.RData')
> table(sce.zeisel$level1class)
astrocytes_ependymal endothelial-mural
179 160
interneurons microglia
290 78
oligodendrocytes pyramidal CA1
774 938
pyramidal SS
397
library(scran)
wilcox.z <- pairwiseWilcox(sce.zeisel, sce.zeisel$level1class,
lfc=1, direction="up")
markers.z <- getTopMarkers(wilcox.z$statistics, wilcox.z$pairs,
pairwise=FALSE, n=50)
> lengths(markers.z)
astrocytes_ependymal endothelial-mural
79 83
interneurons microglia
118 69
oligodendrocytes pyramidal CA1
81 125
pyramidal SS
149
测试数据集 sce.tasic Tasic et al. (2016)
library(scRNAseq)
sce.tasic <- TasicBrainData()
sce.tasic
## class: SingleCellExperiment
## dim: 24058 1809
## metadata(0):
## assays(1): counts
## rownames(24058): 0610005C13Rik 0610007C21Rik ... mt_X57780 tdTomato
## rowData names(0):
## colnames(1809): Calb2_tdTpositive_cell_1 Calb2_tdTpositive_cell_2 ...
## Rbp4_CTX_250ng_2 Trib2_CTX_250ng_1
## colData names(13): sample_title mouse_line ... secondary_type
## aibs_vignette_id
## reducedDimNames(0):
## altExpNames(1): ERCC
接下来使用AUCell 这个包
首先拿到参考数据集中的marker基因集markers.z，并构建一个GeneSetCollection对象
然后对测试数据集中基因表达量进行排名，最后根据排名和之前的基因集进行统计，得到每个marker基因集在每个细胞中的area under the curve (AUC)指标，就能得出哪一个marker基因集在测试数据集的细胞中占比最高
> names(markers.z)
[1] "astrocytes_ependymal" "endothelial-mural"
[3] "interneurons" "microglia"
[5] "oligodendrocytes" "pyramidal CA1"
[7] "pyramidal SS"
# 对marker基因集进行操作
library(GSEABase)
all.sets <- lapply(names(markers.z), function(x) {
GeneSet(markers.z[[x]], setName=x)
})
all.sets <- GeneSetCollection(all.sets)
> all.sets
GeneSetCollection
names: astrocytes_ependymal, endothelial-mural, ..., pyramidal SS (7 total)
unique identifiers: Apoe, Clu, ..., Gabrb2 (555 total)
types in collection:
geneIdType: NullIdentifier (1 total)
collectionType: NullCollection (1 total)
# 对测试数据集进行操作
library(AUCell)
# 基因表达量排名
rankings <- AUCell_buildRankings(counts(sce.tasic),
plotStats=FALSE, verbose=FALSE)
# 计算AUC(结果见下图)
cell.aucs <- AUCell_calcAUC(all.sets, rankings)
# 转置
results <- t(assay(cell.aucs))
这样就能清楚看到，每个细胞对应最高的AUC值属于哪个marker基因集，因此可以作为这个细胞的label
结果的检查 ——方法一：与参考数据中的cluster类型作比较
new.labels <- colnames(results)[max.col(results)]
table(new.labels, sce.tasic$broad_type)
可以看到一个cluster中的大部分细胞，都被赋予了一个相同的label
结果的检查 ——方法二：计算AUC最大值减去中位数
先说判断方法：如果注释的很明确，这个差值会较大；较小的差值说明结果不理想；当然我们可以自己定义一个阈值，过滤掉这些无法注释的细胞
下面的例子中，就是假定得到的大多都是注释较好的细胞，然后得到这个阈值
library(scater)
library(DelayedArray)
deltas <- rowMaxs(results) - rowMedians(results)
discard <- isOutlier(deltas, type="lower", batch=new.labels)
table(new.labels[discard])
##
## astrocytes_ependymal endothelial-mural interneurons
## 24 1 7
## oligodendrocytes pyramidal SS
## 10 16
看看这些被过滤掉的，在整体分布中处于什么位置
par(mar=c(10,4,1,1))
boxplot(split(deltas, new.labels), las=2)
points(attr(discard, "thresholds")[1,], col="red", pch=4, cex=2)
AUCell 包的好处在于不依赖表达量
只需要获得基因列表（可以基于文献或其他背景知识获得），例如从MsigDb上下载
library(BiocFileCache)
bfc <- BiocFileCache(ask=FALSE)
scsig.path <- bfcrpath(bfc, file.path("http://software.broadinstitute.org",
"gsea/msigdb/supplemental/scsig.all.v1.0.symbols.gmt"))
scsigs <- getGmt(scsig.path)
> scsigs
GeneSetCollection
names: Zheng_Cord_Blood_C1_Putative_Megakaryocyte_Progenitor, Zheng_Cord_Blood_C2_Putative_Basophil_Eosinophil_Mast_cell_Progenitor, ..., Muraro_Pancreas_Endothelial_Cell (256 total)
unique identifiers: ABCC3, ABCC4, ..., SEC31A (18656 total)
types in collection:
geneIdType: NullIdentifier (1 total)
collectionType: NullCollection (1 total)
不基于表达量有好有坏，好处就是避免了表达量出现的噪音干扰，处理也更快；坏处是可能会丢掉发现更精细的亚群可能性。
再看一眼另一个数据，不需要运行重点看看差异即可
# 先对测试数据集基因排名
muraro.mat <- counts(sce.muraro)
rownames(muraro.mat) <- rowData(sce.muraro)$symbol
muraro.rankings <- AUCell_buildRankings(muraro.mat,
plotStats=FALSE, verbose=FALSE)
# 然后把MsigDb的scsigs应用在测试数据上
scsig.aucs <- AUCell_calcAUC(scsigs, muraro.rankings)
scsig.results <- t(assay(scsig.aucs))
# 准备作图数据1
full.labels <- colnames(scsig.results)[max.col(scsig.results)]
tab <- table(full.labels, sce.muraro$label)
fullheat <- pheatmap(log10(tab+10), color=viridis::viridis(100), silent=TRUE)
# 然后把MsigDb的一个子集（Pancreas相关）的scsigs应用在测试数据上
scsigs.sub <- scsigs[grep("Pancreas", names(scsigs))]
sub.aucs <- AUCell_calcAUC(scsigs.sub, muraro.rankings)
sub.results <- t(assay(sub.aucs))
# 准备作图数据2
sub.labels <- colnames(sub.results)[max.col(sub.results)]
tab <- table(sub.labels, sce.muraro$label)
subheat <- pheatmap(log10(tab+10), color=viridis::viridis(100), silent=TRUE)
# 组合
gridExtra::grid.arrange(fullheat[[4]], subheat[[4]])
看到一个明显的变化就是：取了MsigDb关于Pancreas的子集后，效果变好了；而使用全部的基因集会导致混淆。因此使用哪些基因还是要基于我们对自己测试数据集的认识
3 基于marker基因的富集分析
思路：还是先得到每个cluster的标志性marker基因，然后对这些代表该cluster的基因们进行富集分析，看看它们集中在哪些通路
示例数据
再以一个新的小鼠乳腺数据 Bach et al. (2017) 为例，还是经历了：基因ID注释=》质控=》归一化=》找HVGs=》降维=》聚类=》找marker基因
数据可以直接下载
#--- loading ---#
library(scRNAseq)
sce.mam <- BachMammaryData(samples="G_1")
#--- gene-annotation ---#
library(scater)
rownames(sce.mam) <- uniquifyFeatureNames(
rowData(sce.mam)$Ensembl, rowData(sce.mam)$Symbol)
library(AnnotationHub)
ens.mm.v97 <- AnnotationHub()[["AH73905"]]
rowData(sce.mam)$SEQNAME <- mapIds(ens.mm.v97, keys=rowData(sce.mam)$Ensembl,
keytype="GENEID", column="SEQNAME")
#--- quality-control ---#
is.mito <- rowData(sce.mam)$SEQNAME == "MT"
stats <- perCellQCMetrics(sce.mam, subsets=list(Mito=which(is.mito)))
qc <- quickPerCellQC(stats, percent_subsets="subsets_Mito_percent")
sce.mam <- sce.mam[,!qc$discard]
#--- normalization ---#
library(scran)
set.seed(101000110)
clusters <- quickCluster(sce.mam)
sce.mam <- computeSumFactors(sce.mam, clusters=clusters)
sce.mam <- logNormCounts(sce.mam)
#--- variance-modelling ---#
set.seed(00010101)
dec.mam <- modelGeneVarByPoisson(sce.mam)
top.mam <- getTopHVGs(dec.mam, prop=0.1)
#--- dimensionality-reduction ---#
library(BiocSingular)
set.seed(101010011)
sce.mam <- denoisePCA(sce.mam, technical=dec.mam, subset.row=top.mam)
sce.mam <- runTSNE(sce.mam, dimred="PCA")
#--- clustering ---#
snn.gr <- buildSNNGraph(sce.mam, use.dimred="PCA", k=25)
colLabels(sce.mam) <- factor(igraph::cluster_walktrap(snn.gr)$membership)
#--- find marker ---#
markers.mam <- findMarkers(sce.mam, direction="up", lfc=1)
然后我们想看cluster2是什么细胞类型
首先得到cluster2的marker基因，并且还挑差异显著的基因
chosen <- "2"
cur.markers <- markers.mam[[chosen]]
is.de <- cur.markers$FDR <= 0.05
summary(is.de)
## Mode FALSE TRUE
## logical 27819 179
进行GO富集分析
library(org.Mm.eg.db)
library(clusterProfiler)
gn=unique(rownames(cur.markers)[is.de])
trans=bitr(gn,fromType = "SYMBOL",
toType = 'ENTREZID',
OrgDb = org.Mm.eg.db)
go.out=enrichGO(gene = trans$ENTREZID,
keyType = 'ENTREZID',
OrgDb = org.Mm.eg.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.05,
readable = TRUE)
View(head(summary(go.out),20))
接下来就是考验自己对数据的认知了，假设我们这里对这个数据很熟悉。那么看到这里基因主要参与了lipid synthesis、 cell adhesion 、 tube formation。再结合乳腺的研究背景，推测cluster2可能包含与乳汁产生和分泌相关的管腔上皮细胞，进一步检查一下相关的基因Csn2和Csn3，也是上调的
plotExpression(sce.mam, features=c("Csn2", "Csn3"),
x="label", colour_by="label")
如果想看cluster2中某个感兴趣通路中的基因表现情况
# 比如想看与cell adhesion相关的通路
library(stringr)
adhesion <- go.out["GO:0022408",]
adhesion_gn=str_split(adhesion$geneID,"/",simplify = T)
head(cur.markers[rownames(cur.markers) %in% adhesion_gn,1:4], 10)
## DataFrame with 10 rows and 4 columns
## Top p.value FDR summary.logFC
## <integer> <numeric> <numeric> <numeric>
## Spint2 11 3.28234e-34 1.37163e-31 2.39280
## Epcam 17 8.86978e-94 7.09531e-91 2.32968
## Cebpb 21 6.76957e-16 2.03800e-13 1.80192
## Cd24a 21 3.24195e-33 1.29669e-30 1.72318
## Btn1a1 24 2.16574e-13 6.12488e-11 1.26343
## Cd9 51 1.41373e-11 3.56592e-09 2.73785
## Ceacam1 52 1.66948e-38 7.79034e-36 1.56912
## Sdc4 59 9.15001e-07 1.75467e-04 1.84014
## Anxa1 68 2.58840e-06 4.76777e-04 1.29724
## Cdh1 69 1.73658e-07 3.54897e-05 1.31265
TODO：补充12.5 Computing gene set activities
