# 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 ``` 这个参数对上面的三种统计方法都适用 --- # Agent Instructions: Querying This Documentation If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question. Perform an HTTP GET request on the current page URL with the `ask` query parameter: ``` GET https://jieandze1314.osca.top/03/03-6.md?ask= ``` The question should be specific, self-contained, and written in natural language. The response will contain a direct answer to the question and relevant excerpts and sources from the documentation. Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.