对应原版教程第11章http://bioconductor.org/books/release/OSCA/overview.html
从基因表达水平来看，主要是哪些基因造成了每个cluster得以区分的结果是这一步的主要目的。实现这一步后，对于cluster的生物学意义解读将有很大帮助。
源网，侵删~
笔记要点
1、cluster marker gene的意义
2、基于t检验
3、其它检验方法
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
- 以cluster 7的分析结果为例
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 多重检验） -
FDR：p.value的FDR值 -
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)
- 如下图结果，少了
top列，原因很容易理解。 -
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往往会忽略这个因素。
