3.6 Marker/标记基因检测
刘小泽写于2020.7.3+7.5
1 前言
在上一章聚类分群的结尾，为了解释分群的结果，指定了几个基因进行区分，这几个基因就属于marker基因或者叫标志基因，它们是经过反复验证得到的。也就是说，一般看到相关的marker基因，就可以把某个cluster与某种细胞类型对应起来；另外这个思路还可以探索亚群之间发生的微小差异（例如通路激活、分化状态）与基因表达的联系
上面👆说的，主要还是验证，就是看看某个cluster到底是不是这个细胞类型，做这个的前提是你已经了解了你的细胞是什么类型，只是不确定。
对于常规流程来讲，更多的是探索，就是我们得到了分群的结果，然后再怎么分析？怎么和marker基因联系起来？数万个基因，哪些才是这个cluster的marker基因？
我们认为，导致cluster出现差异的那部分基因中，尤其是那些差异最显著的基因中，最可能包含marker基因。因此一个首要任务就是去进行cluster之间的差异分析，最后把每个cluster中最显著的前多少基因拿出来，放在一起（有没有想起来Seurat的那个热图？黑黄相间的那个）
下面还是使用10X PBMC数据
这次就不从头开始处理了，如果想知道从数据读取到聚类之间的步骤，可以去看之前写的，里面有详细处理代码。这次直接加载聚类后的数据即可：clustered.sce.pbmc.RData 链接：https://share.weiyun.com/4fOP3IDw 密码：xswtct
load('clustered.sce.pbmc.RData')
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):2 检测方法
2.1 常规方法-t检验
findMarkers()提供了三种检验方法，分别是：t、wilcox、binom。默认使用t检验
针对大量细胞的检验，Welch t-test计算速度很快，并且有不错的统计学意义。使用findMarkers() 可以对每个基因在clusters之间进行两两比较，返回的列表中包括了每个cluster中排过序的候选marker基因
findMarkers() 默认根据sce.pbmc中的colLabels()信息提取各个cluster，使用的就是groups参数，因此也可以写成：
有了这个参数，就可以根据其他方法的分群结果再去探索不一样的marker基因，而不用修改原始的sce.pbmc结果，只要groups参数满足与sce.pbmc列数一致就可以
怎么解释这个结果呢？--以cluster9为例
可以看到其中有很多logFC值，表示log2(cluster9/other_cluster)
findMarkers() 默认会根据第一列Top进行排序，而其中第一列”Top“指的是：cluster9中差异表达的前多少位基因。但Top1不代表只有1个基因，而是可能有很多并列第一，但为了区分，还是用pvalue给各个并列第一又排了顺序
知道了这个概念，提取cluster9中的Top6的所有基因也不是难事
看到其中platelet factor 4 (PF4) 与 pro-platelet basic protein (PPBP)基因表达量高，推测cluster9含有血小板或者它的前体
除了Top这一列，还有一列叫”summary.logFC“，它可以帮我们快速判断基因在cluster9中是上调还是下调，例如看PF4在这里的logFC值就很高
不同的比较方法
注意到，这里使用的差异比较方法是”两两比较“，即将一个cluster和其余各个cluster进行比较。
当然还有其他方法是将一个cluster与其他剩余cluster的均值进行比较。如果是与均值比较，那么就会细胞群体组成比较敏感，如果其中出现一个”支配欲超强的“cluster，那么它会把整个平均的表达水平带偏，结果看到的差异分析也并不准确。
另外，使用两两比较的方法还会提供有关marker基因更多的信息，比如能看到哪些clusters是由某个marker基因区分的
2.2 如果只关注上调基因
之前findMarkers目的是选取上调、下调基因作为marker基因的候选，但其实下调基因很难吸引我们的注意力，我们会首先关注图上红色的，也就是上调的基因们。另一方面，相比于上调，下调基因也很难通过实验去验证。因此，如果只关心上调的基因，可以使用单边t检验，将某个cluster和其他clusters进行比较，设置参数direction="up"
另外还可以根据logFC阈值过滤，其实这些都可以自己手动过滤，多个参数只是更方便一点。缺点就是需要记住这个函数有这个参数
根据两重过滤的结果（direction + logFC），再看cluster9的marker基因热图
注意
我们这里只是探索了上调基因。至于有一些亚群中可能部分基因下调才导致这个亚群与众不同，这样的亚群这里检测不到
这里的阈值可能会漏掉一些改变幅度不大，但依然重要的基因
2.3 寻找cluster特异的marker基因
上面提到，
findMarkers()会对两两比较结果做一个排名，然后选择p值比较显著的一些作为Top基因，返回的结果包括了所有cluster的logFC情况，比如这里有18个cluster，那么返回的结果也包含18个cluster的logFC。以上调差异表达为例：某个基因在cluster9与cluster1相比下上调，但这个基因在cluster9和cluster2中不上调，依然会把这个基因列出来
有一种更严格的过滤机制，就是只选在某个cluster与其他各个clusters相比都差异表达的基因，结果返回17个cluster（不包括自己）。以上调差异表达为例：结果得到的每个基因，不管是cluster9与cluster1比较，还是与cluster2比较，都是上调的；并且，这个基因仅仅在cluster9中是上调的
这种做法过于严格，以至于如果一个基因除了在cluster9中上调，还在cluster4中上调，这个基因就不会写进结果。【它的想法很简单，就是单纯针对cluster9，只找在它里面上调的】
另外，如果细胞分群效果不好，这样的寻找方法会过滤掉太多的潜在的感兴趣基因
举个例子：如果一群细胞混杂了单纯CD4+、单纯CD8+、二者都有、二者都无这四种情况。如果设置
pval.type="all"，那么Cd4或Cd8基因都不会列入marker基因结果，因为它们都会在两个亚群有差异表达情况
还有另一种方法：pval.type="any" ，就是只要在一个cluster中出现差异表达，就写入结果。但这个设置有有点过于宽松，因此一个折中的方法：pval.type="some" ，各个方法得到的summary.logFC这一列都是不同的，基因排名前后也略有变化。
我相信大家可能并关心每个方法是怎么去比较，怎么去做统计的，更关心的是，我用了这个方法，得到的结果可用性有多高？这个问题可以思考一下，如果真的有一种普适性的方法，开发者又怎么会去设置这么多参数呢？直接一个参数到底，不是更方便？
数据分析就是这样，存在太多的不确定性。但只有我们使用一种方法感觉走不通了，才会意识到多个参数多条路的滋味。
2.4 除了t检验，还有Wilcoxon、binomial检验
首先来了解一下各种统计检验函数
参考：https://shixiangwang.github.io/home/cn/post/2019-12-25-r-stats-funs-summary/
对于连续型数据
基于正态分布的检验：
均值检验：
t.test()两总体方差检验：
var.test()多个组间均值的比较（ANOVA）：
aov()多组样本的配对 t 检验：
pairwise.t.test()正态性检验：
shapiro.test()分布的对称性检验：
ks.test()检验两个向量是否服从同一分布：
ks.test()相关性检验：
cor.test()
不依赖分布的检验:
均值检验：Wilcoxon ，也就是说，Wilcoxon 检验是 t 检验的非参数版本。默认是秩和检验
wilcox.test多均值比较：
kruskal.test()Kruskal-Wallis 秩和检验方差检验：
fligner.test()Fligner-Killeen（中位数）检验完成不同组别的方差比较尺度参数差异：
ansari.test()针对两样本尺度参数差异的 Ansari-Bradley 检验；mood.test()使用 Mood 两样本检验
对于离散型数据
比例检验：
prop.test()比较两组观测值发生的概率是否有差异二项式检验：
binom.test()列联表检验：
fisher.test()用来确定两个分类变量是否相关； 针对小的列联表，Fisher 精确检验效果不错；大列联表可以用卡方检验代替；检验三个变量的混合影响，可以用Cochran-Mantel-Haenszel 检验mantelhaen.test()；McNemar 卡方可以检验二维列联表的对称性mcnemar.test()列联表非参数检验：Friedman 秩和检验（非参数的双边 ANOVA 检验）
friedman.test()
找marker基因的Wilcoxon方法
找marker基因的binomial方法
2.5 整合多个统计分析的结果
这样可以检验哪些基因是强有力的marker基因（能撑得住三大检验方法的考验）
当然也可以重点关注其中的某一个指标，比如wilcox的结果AUC越大，就表示基因表达量在各个cluster之间的分布越分散；t-test的logFC结果越大，表示更容易解释一个基因在两个cluster之间的变化幅度
3 另外，可以封锁一些不重要因素
大型数据集一般会涉及许多变化因素（比如批次、性别差异、年龄差异等等），它们会对数据波动产生影响，但又不是感兴趣的生物因素。如果在检测marker基因的时候带着它们，可能会结果产生干扰。
因此有一个参数block可以帮我们”锁住“它们，之后检测的时候就不会把这些因素纳入考量
这个参数对上面的三种统计方法都适用
最后更新于
这有帮助吗？