# 3.6 Marker/标记基因检测 ## 1 前言 在上一章聚类分群的结尾，为了解释分群的结果，指定了几个基因进行区分，这几个基因就属于**marker基因或者叫标志基因**，它们是经过反复验证得到的。也就是说，一般看到相关的marker基因，就可以把某个cluster与某种细胞类型对应起来；另外这个思路还可以探索亚群之间发生的微小差异（例如通路激活、分化状态）与基因表达的联系 上面👆说的，主要还是验证，就是看看某个cluster到底是不是这个细胞类型，做这个的前提是你已经了解了你的细胞是什么类型，只是不确定。 对于常规流程来讲，更多的是探索，就是我们**得到了分群的结果，然后再怎么分析？怎么和marker基因联系起来？数万个基因，哪些才是这个cluster的marker基因？** 我们认为，导致cluster出现差异的那部分基因中，尤其是那些差异最显著的基因中，最可能包含marker基因。因此一个首要任务就是去进行cluster之间的差异分析，最后把每个cluster中最显著的前多少基因拿出来，放在一起（有没有想起来Seurat的那个热图？黑黄相间的那个） ### **下面还是使用10X PBMC数据** 这次就不从头开始处理了，如果想知道从数据读取到聚类之间的步骤，可以去看之前写的，里面有详细处理代码。这次直接加载聚类后的数据即可：`clustered.sce.pbmc.RData` 链接： 密码：xswtct ```r 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基因 ```r library(scran) markers.pbmc <- findMarkers(sce.pbmc) markers.pbmc ## List of length 18 ## names(18): 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ``` `findMarkers()` 默认根据`sce.pbmc`中的`colLabels()`信息提取各个cluster，使用的就是`groups`参数，因此也可以写成： ```r same.markers <- findMarkers(sce.pbmc, groups=colLabels(sce.pbmc)) ``` 有了这个参数，就可以根据其他方法的分群结果再去探索不一样的marker基因，而不用修改原始的`sce.pbmc`结果，只要`groups`参数满足与`sce.pbmc`列数一致就可以 **怎么解释这个结果呢？--以cluster9为例** ```r chosen <- "9" interesting <- markers.pbmc[[chosen]] colnames(interesting) ## [1] "Top" "p.value" "FDR" "summary.logFC" ## [5] "logFC.1" "logFC.2" "logFC.3" "logFC.4" ## [9] "logFC.5" "logFC.6" "logFC.7" "logFC.8" ## [13] "logFC.10" "logFC.11" "logFC.12" "logFC.13" ## [17] "logFC.14" "logFC.15" "logFC.16" "logFC.17" ## [21] "logFC.18" ``` 可以看到其中有很多logFC值，表示`log2(cluster9/other_cluster)` `findMarkers()` 默认会根据第一列Top进行排序，而其中第一列”Top“指的是：cluster9中差异表达的前多少位基因。**但Top1不代表只有1个基因，而是可能有很多并列第一**，但为了区分，还是用pvalue给各个并列第一又排了顺序 ```r # 比如这里看到的，并列第一就不止10个 interesting[1:10,1:4] ## DataFrame with 10 rows and 4 columns ## Top p.value FDR summary.logFC ## ## S100A4 1 3.29706e-57 3.05195e-55 -4.52198 ## TAGLN2 1 1.65522e-24 3.58425e-23 4.83531 ## PF4 1 2.54870e-35 9.99719e-34 5.91366 ## GZMA 1 1.41952e-120 7.71441e-118 -1.95444 ## HLA-DQA1 1 1.79189e-88 4.75402e-86 -3.64622 ## FCN1 1 1.13468e-246 4.77901e-243 -2.81179 ## SERPINA1 1 1.12795e-68 1.72751e-66 -2.43278 ## RPL23A 1 2.42151e-37 1.04737e-35 -4.07367 ## RPL17 1 0.00000e+00 0.00000e+00 -2.82856 ## RPS21 1 1.08454e-56 9.90316e-55 -3.99499 ``` 知道了这个概念，提取cluster9中的Top6的所有基因也不是难事 ```r best.set <- interesting[interesting$Top <= 6,] # 看到Top6其实总共有53个基因 > dim(best.set) [1] 53 21 # 提取这些基因在与各个cluster比较时的logFC logFCs <- getMarkerEffects(best.set) # 总共18个cluster，除去自己还剩17个 > dim(logFCs) [1] 53 17 > logFCs[1:4,1:4] 1 2 3 4 S100A4 -2.7287830 -0.89602358 -2.887000 -4.33065808 TAGLN2 3.8690448 3.72075587 3.955002 4.14646142 PF4 6.0988107 6.10290759 6.103301 6.07864054 GZMA -0.4340303 -0.09184649 -1.954439 -0.07511017 library(pheatmap) pheatmap(logFCs, breaks=seq(-5, 5, length.out=101)) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-03-091314.png) 看到其中platelet factor 4 (PF4) 与 pro-platelet basic protein (PPBP)基因表达量高，推测cluster9含有血小板或者它的前体 **除了Top这一列，还有一列叫”summary.logFC“**，它可以帮我们快速判断基因在cluster9中是上调还是下调，例如看PF4在这里的logFC值就很高 ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-03-091940.png) **不同的比较方法** * 注意到，**这里使用的差异比较方法是”两两比较**“，即将一个cluster和其余各个cluster进行比较。 * 当然还有**其他方法是将一个cluster与其他剩余cluster的均值进行比较**。如果是与均值比较，那么就会细胞群体组成比较敏感，如果其中出现一个”支配欲超强的“cluster，那么它会把整个平均的表达水平带偏，结果看到的差异分析也并不准确。 * 另外，使用两两比较的方法还会提供有关marker基因更多的信息，比如能看到哪些clusters是由某个marker基因区分的 ### **2.2 如果只关注上调基因** 之前`findMarkers`目的是选取上调、下调基因作为marker基因的候选，但其实下调基因很难吸引我们的注意力，**我们会首先关注图上红色的，也就是上调的基因们**。另一方面，相比于上调，下调基因也很难通过实验去验证。因此，如果只关心上调的基因，可以**使用单边t检验，**将某个cluster和其他clusters进行比较，设置参数`direction="up"` ```r markers.pbmc.up <- findMarkers(sce.pbmc, direction="up") interesting.up <- markers.pbmc.up[[chosen]] interesting.up[1:10,1:4] ## DataFrame with 10 rows and 4 columns ## Top p.value FDR summary.logFC ## ## TAGLN2 1 8.27609e-25 9.29516e-21 4.83531 ## PF4 1 1.27435e-35 4.29379e-31 5.91366 ## SDPR 2 2.26416e-21 1.90722e-17 4.72820 ## GPX1 2 1.79269e-20 1.00671e-16 4.83143 ## TMSB4X 2 1.61389e-31 2.71891e-27 3.71343 ## PPBP 3 2.67043e-20 1.28539e-16 5.54885 ## NRGN 3 1.41986e-20 9.56813e-17 4.18416 ## CCL5 5 2.55331e-18 9.55903e-15 4.62327 ## GNG11 6 2.06623e-18 8.70243e-15 4.73606 ## HIST1H2AC 7 1.05437e-17 3.55260e-14 4.76160 ``` 另外还可以根据logFC阈值过滤，其实这些都可以自己手动过滤，多个参数只是更方便一点。缺点就是需要记住这个函数有这个参数 ```r markers.pbmc.up2 <- findMarkers(sce.pbmc, direction="up", lfc=1) interesting.up2 <- markers.pbmc.up2[[chosen]] interesting.up2[1:10,1:4] ``` 根据两重过滤的结果（direction + logFC），再看cluster9的marker基因热图 ```r best.set <- interesting.up2[interesting.up2$Top <= 5,] logFCs <- getMarkerEffects(best.set) library(pheatmap) pheatmap(logFCs, breaks=seq(-5, 5, length.out=101)) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-03-111225.png) **注意** * 我们这里只是探索了上调基因。至于有一些亚群中可能部分基因下调才导致这个亚群与众不同，这样的亚群这里检测不到 * 这里的阈值可能会漏掉一些改变幅度不大，但依然重要的基因 ### **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中是上调的 ```r # 设置pval.type="all"就是做了这件事，并且还是选上调的基因 markers.pbmc.up3 <- findMarkers(sce.pbmc, pval.type="all", direction="up") # 然后就想看看感兴趣的cluster9与其他各个clusters比较后差异基因 interesting.up3 <- markers.pbmc.up3[[chosen]] interesting.up3[1:10,1:3] ``` * 这种做法过于严格，以至于如果一个基因除了在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检验** **首先来了解一下各种统计检验函数** > 参考： **对于连续型数据** 基于正态分布的检验： * 均值检验：`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方法** ```r 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] "16" "17" "18" interesting.wmw <- markers.pbmc.wmw[[chosen]] interesting.wmw[1:10,1:4] ## DataFrame with 10 rows and 4 columns ## Top p.value FDR summary.AUC ## ## PF4 1 3.13749e-164 1.05715e-159 0.988833 ## TMSB4X 1 5.07215e-27 2.05905e-24 0.992149 ## SDPR 2 2.12114e-145 3.57349e-141 0.955218 ## NRGN 2 1.18240e-131 7.96793e-128 0.966119 ## TAGLN2 3 1.55560e-28 6.98860e-26 0.967186 ## PPBP 3 3.57148e-134 4.01125e-130 0.932743 ## GNG11 3 2.46077e-126 1.38189e-122 0.932491 ## TUBB1 3 7.55573e-133 6.36457e-129 0.921632 ## HIST1H2AC 4 4.69094e-94 1.43688e-90 0.930973 ## ACTB 5 1.53723e-23 5.28523e-21 0.949431 ``` **找marker基因的binomial方法** ```r markers.pbmc.binom <- findMarkers(sce.pbmc, test="binom", direction="up") # 选出cluster9的marker基因 interesting.binom <- markers.pbmc.binom[[chosen]] colnames(interesting.binom) ## [1] "Top" "p.value" "FDR" "summary.logFC" ## [5] "logFC.1" "logFC.2" "logFC.3" "logFC.4" ## [9] "logFC.5" "logFC.6" "logFC.7" "logFC.8" ## [13] "logFC.10" "logFC.11" "logFC.12" "logFC.13" ## [17] "logFC.14" "logFC.15" "logFC.16" "logFC.17" ## [21] "logFC.18" library(scater) top.genes <- head(rownames(interesting.binom)) plotExpression(sce.pbmc, x="label", features=top.genes) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-05-055703.png) ### **2.5 整合多个统计分析的结果** 这样可以检验哪些基因是强有力的marker基因（能撑得住三大检验方法的考验） ```r combined <- multiMarkerStats( t=findMarkers(sce.pbmc, direction="up"), wilcox=findMarkers(sce.pbmc, test="wilcox", direction="up"), binom=findMarkers(sce.pbmc, test="binom", direction="up") ) > combined List of length 18 names(18): 1 2 3 4 5 6 7 8 ... 12 13 14 15 16 17 18 head(combined[["9"]][,1:15]) ``` ![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-05-061802.png) 当然也可以重点关注其中的某一个指标，比如wilcox的结果AUC越大，就表示基因表达量在各个cluster之间的分布越分散；t-test的logFC结果越大，表示更容易解释一个基因在两个cluster之间的变化幅度 ## 3 另外，可以封锁一些不重要因素 大型数据集一般会涉及许多变化因素（比如批次、性别差异、年龄差异等等），它们会对数据波动产生影响，但又不是感兴趣的生物因素。如果在检测marker基因的时候带着它们，可能会结果产生干扰。 因此**有一个参数`block`可以帮我们”锁住“它们**，之后检测的时候就不会把这些因素纳入考量 ```r # 举个例子，在416b数据集中有批次因素，可以锁定 m.out <- findMarkers(sce.416b, block=sce.416b$block, direction="up") demo <- m.out[["1"]] # 前Top5的marker基因 demo[demo$Top <= 5,1:4] ## DataFrame with 13 rows and 4 columns ## Top p.value FDR summary.logFC ## ## Foxs1 1 1.37387e-12 4.35563e-10 3.07058 ## Pirb 1 2.08277e-33 1.21332e-29 5.87820 ## Myh11 1 6.44327e-47 3.00282e-42 4.38182 ## Tmsb4x 2 3.22944e-44 7.52525e-40 1.47689 ## Ctsd 2 6.78109e-38 7.90065e-34 2.89152 ## ... ... ... ... ... ## Tob1 4 6.63870e-09 1.18088e-06 2.74161 ## Pi16 4 1.69247e-32 7.88758e-29 5.76914 ## Cd53 5 1.08574e-27 2.97646e-24 5.75200 ## Alox5ap 5 1.33791e-28 4.15679e-25 1.36676 ## CBFB-MYH11-mcherry 5 3.75556e-35 3.50049e-31 3.01677 ``` 这个参数对上面的三种统计方法都适用