OSCA单细胞数据分析笔记10—Marker gene detection

OSCA单细胞数据分析笔记10—Marker gene detection by 单细胞天地

分享是一种态度

对应原版教程第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往往会忽略这个因素。