对应原版教程第11章: http:///books/release/OSCA/overview.html
从基因表达水平来看,主要是哪些基因造成了每个cluster得以区分的结果是这一步的主要目的。实现这一步后,对于cluster的生物学意义解读将有很大帮助。 笔记要点
1、cluster marker gene的意义- 物以类聚,人以群分。细胞因gene表达相似而归为一类;因不同的表达特征而形成多个cluster。因此,这一步的目的就是找出造成每个cluster与其它clusters分离的marker gene---每个cluster的特征基因集(高表达/低表达)。可以为下一步的cluster细胞类型注释提供主要参考信息。
- cluster marker gene的鉴定思路类似于 Bulk RNA-seq的差异分析。不同之处在于往往会有十多个cluster,因此定义cluster的marker gene有蛮多细节值得考虑。
- 并且由于scRNA-seq与Bulk RNA-seq的不同,一般不推荐使用诸如DESeq2、edgeR等差异分析思路(原因在后面会提到)。本篇笔记主要学习了使用
scran 包的findMarkers() 函数提供的基于不同method的鉴定思路。
#示例数据 sce.pbmc #获取方式见原教程 #class: SingleCellExperiment #dim: 33694 3985
table(sce.pbmc$label) # 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 #205 508 541 56 374 125 46 432 302 867 47 155 166 61 84 16
2、基于t检验- Welch t-test是基于两组的均值是否相同的假设检验思路;为
findMarkers() 函数的默认方法 - 比较思路是首先是选择其中的一个cluster与其余所有的cluster分别两两比较。就比较结果而言,定义cluster的marker gene有如下三个角度可供选择(对应的零假设也稍有不同)
2.1 any differences(generous)- 零假设:对于基因A,cluster X与其余全部cluster的表达均值均相同
- 换句话说:只要cluster X的基因A的表达水平与其余任意一个cluster的t检验结果显著,就认为是该cluster的marker gene。
library(scran) markers.pbmc <- findMarkers(sce.pbmc) markers.pbmc #分别储存每个cluster的分析结果 #List of length 16 #names(16): 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
chosen <- "7" interesting <- markers.pbmc[[chosen]] colnames(interesting) # [1] "Top" "p.value" "FDR" "summary.logFC" "logFC.1" # [6] "logFC.2" "logFC.3" "logFC.4" "logFC.5" "logFC.6" # [11] "logFC.8" "logFC.9" "logFC.10" "logFC.11" "logFC.12" # [16] "logFC.13" "logFC.14" "logFC.15" "logFC.16"
head(interesting[,1:5],15) # DataFrame with 15 rows and 5 columns # Top p.value FDR summary.logFC logFC.1 # <integer> <numeric> <numeric> <numeric> <numeric> # S100A4 1 2.59737e-38 1.27018e-36 -4.27560 -2.1597143 # TAGLN2 1 8.65033e-28 2.44722e-26 5.07327 4.7714408 # FCGR3A 1 8.84356e-63 1.15048e-60 -3.07121 -1.2852520 # GZMA 1 1.15392e-120 7.20000e-118 -1.92877 -2.5023377 # HLA-DQA1 1 3.43640e-83 8.90663e-81 -3.54890 -0.0220385 # ... ... ... ... ... ... # CTSS 2 1.58027e-33 6.15557e-32 -3.51820 -0.242616 # HOPX 2 4.60496e-78 1.05551e-75 -1.99060 -1.990602 # PF4 2 9.57857e-35 3.97464e-33 6.71811 6.862880 # PRELID1 2 1.10561e-107 5.47832e-105 -1.61240 -0.761206 # AC090498.1 2 1.19635e-156 1.55038e-153 -1.93799 -0.872609
a7=as.data.frame(interesting) dim(a7) #[1] 33694 19 View(a7)
如下图,不同列的释义 Top :cluster 7 分别与其它所有cluster的t-test结果里,P值最显著的基因。按理,应当有15个Top1 gene。但如图可以看到只有10个Top1,我认为是由于存在重复的gene结果导致。并且在前面出现过的基因,再后面的TopX就不再出现了。
p.value :偏离零假设的程度,具体计算是结合特定cluster与其余所有cluster的两两t检验p值的combined p value(Simon method 多重检验)
summary.logFC :即最显著的那次cluster两两比较的logFoldChange
logFC.X :该基因分别与其它所有cluster的logFC结果
这种思路定义的cluster marker gene标准是比较宽松的。因为只有一个基因在任意一次cluster间比较p值显著,就会认为是marker gene
2.2 all difference(stringent)- 零假设:对于基因A,cluster X与其余任一cluster的表达均值均相同
- 换句话说:cluster X的基因A的表达水平与其余所有cluster的t检验结果均显著,才能认为是该cluster的marker gene。可以理解为cluster-specific gene。
markers.pbmc2 <- findMarkers(sce.pbmc, pval.type="all") interesting2 <- markers.pbmc2[[chosen]] colnames(interesting2) # [1] "p.value" "FDR" "summary.logFC" "logFC.1" "logFC.2" # [6] "logFC.3" "logFC.4" "logFC.5" "logFC.6" "logFC.8" # [11] "logFC.9" "logFC.10" "logFC.11" "logFC.12" "logFC.13" # [16] "logFC.14" "logFC.15" "logFC.16
a7_2=as.data.frame(interesting2) View(a7_2)
p.value 为该基因的15次t检验的p值结果中的最大值;summary.logFC 同样与之对应;其余列含义可参考2.1
这种思路定义的cluster marker gene标准是比较严苛的,即当某个cluster的一个基因表达水平与其它所有cluster都具有显著差异时,才会被认为是marker gene。
2.3 middle difference(middle)- 零假设:对于基因A,cluster X与最多其余一半的cluster表达均值均相同
- 换句话说:cluster X的基因A的表达水平与至少其余一般的cluster的t检验结果都显著,才能认为是该cluster的marker gene。
markers.pbmc3 <- findMarkers(sce.pbmc, pval.type="some") interesting3 <- markers.pbmc3[[chosen]] colnames(interesting3) a7_3=as.data.frame(interesting3) View(a7_3)
p.value 为该基因的15次t检验的p值结果排名中间的结果;summary.logFC 同样与之对应;其余列含义可参考2.1
相较于2.1过于宽松和2.2过于严苛的方法,这种思路则采取了折中的标准进行筛选。
2.4 其它参数direction 上调?下调?对于marker gene的生物意义来说,表达上调的marker gene更具有可解释性(细胞类型注释)以及可实验验证性。可通过设置direction 参数,只筛选出上调的marker gene
markers.pbmc.up <- findMarkers(sce.pbmc, direction="up")
但这种筛选方法对于那些由相对下调的特征基因集特征细胞组成的cluster可能就没有预期的结果。 lfc 差异倍数 一般t-test的零假设为均值相同,也可以设置参数lfc 设定差异倍数为零假设阈值,并配合direction 参数可筛选出差异倍数大于指定阈值的高表达基因。
markers.pbmc.up2 <- findMarkers(sce.pbmc, direction="up", lfc=1)
- 批次效应
block /design 对于存在batch的测序情况,可以通过block /design 指定batch,进行分析(每个batch单单独差异分析,然后计算合并 combined pvalue)
以sce.416b 为例
#方式一:block参数 m.out <- findMarkers(sce.416b, block=sce.416b$block, direction="up")
#方式二:design参数 design <- model.matrix(~sce.416b$block) design <- design[,-1,drop=FALSE]
m.alt <- findMarkers(sce.416b, design=design, direction="up")
3、其它检验方法- 之间介绍的
pval.type 、direction 、lfc 等参数也都适用于下面的方法。
3.1 Wilcoxon rank sum test
- 又称为 Wilcoxon-Mann-Whitney test、WMW test,用于评价同一指标在两组中的分离度(separation)。
- 计算思路简单来说就是:对于基因A在cluster1与cluster2的表达情况,任意抽取cluster1的一个细胞的基因A表达量是否高于cluster2的任意一个细胞的基因A表达量。反映在结果中,用AUC(area-under-the-curve)指标反应。该值越接近1,则表示cluster1的基因A表达量显著高于cluster2;越接近0,则表示cluster1的基因A表达量显著低于cluster2。
markers.pbmc.wmw <- findMarkers(sce.pbmc, test="wilcox", direction="up") names(markers.pbmc.wmw) # [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15" "16"
interesting.wmw <- markers.pbmc.wmw[[chosen]] interesting.wmw[1:10,1:4]
WMW方法区别于t-test的优势在于不用考虑两个cluster的细胞数量差异,但同时缺点是计算量比较大。
3.2 binomial test- 这种分析方法重新设置了表达矩阵,凡是表达量大于0的都视为1。即把count表达矩阵转化为0/1矩阵(0表示不表达、1代表表达)。
- 此时的零假设就是基因A在两个cluster的active(有表达)的程度相同。(至于表达量的高低就不管了)
- foldchange的含义:对于基因A,cluster1的
1 的比例与cluster2的1 的比例的foldchange
markers.pbmc.binom <- findMarkers(sce.pbmc, test="binom", direction="up") names(markers.pbmc.binom)
interesting.binom <- markers.pbmc.binom[[chosen]] colnames(interesting.binom) interesting.binom[1:10,1:4]
这种定义marker gene的标准是十分严苛的。筛选出来的marker基因直接反映了cluster对于该基因的表达有无情况。
3.3 关于传统Bulk RNA-seq方法- 对于Bulk RNA-seq的差异分析方法(DESeq2、voom-limma pipeline),还是不太适合适用于scRNA-seq的clust间的两两比较。因为一般将cluster的一个细胞视为1 replication的基础上,每个cluster往往会有数百次重复,不符合DESeq2、voom-limma pipeline等适合有限重复的情况;另一方面DESeq2、voom-limma pipeline等假设每次重复的基因表达情况相似,而scRNA-seq测序过程中的每个细胞的表达水平分布可能受各种因素不相一致,即使在同一个cluster里。
- 即使不做推荐,教程里还是演示了采用voom-limma pipeline循环比较所有cluster间两两差异分析的代码,这里就不展示了。
- 换一个角度思考,
findMarkers() 函数提供的分析方法可以适用于涉及多组的Bulk RNA-seq差异分析方法。
关于replication 这个问题,在上述所有方法中是把每个cluster的cell视为一次重复,但这还是一个值得细思的问题。即使这些细胞来自同一个取样组织,但本质上还是不符合生物学重复的概念。换句话说,一次测序的结果应当视为一次重复,对于Bulk RNA-seq比较容易理解,对于scRNA-seq往往会忽略这个因素。
|