查看原文
其他

单细胞交响乐8-marker基因检测

豆豆花花 单细胞天地 2022-06-07
以下文章源自于生信星球 ,作者:刘小泽

原书名为:Orchestrating Single-Cell Analysis

我给它取名叫“单细胞交响乐”
因为单细胞分析就像一个大乐团,需要各个流程的协同配合

1 前言

在上一章聚类分群的结尾,为了解释分群的结果,指定了几个基因进行区分,这几个基因就属于marker基因或者叫标志基因,它们是经过反复验证得到的。也就是说,一般看到相关的marker基因,就可以把某个cluster与某种细胞类型对应起来;另外这个思路还可以探索亚群之间发生的微小差异(例如通路激活、分化状态)与基因表达的联系

上面👆说的,主要还是验证,就是看看某个cluster到底是不是这个细胞类型,做这个的前提是你已经了解了你的细胞是什么类型,只是不确定。

对于常规流程来讲,更多的是探索,就是我们得到了分群的结果,然后再怎么分析?怎么和marker基因联系起来?数万个基因,哪些才是这个cluster的marker基因?

我们认为,导致cluster出现差异的那部分基因中,尤其是那些差异最显著的基因中,最可能包含marker基因。因此一个首要任务就是去进行cluster之间的差异分析,最后把每个cluster中最显著的前多少基因拿出来,放在一起

下面还是使用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基因

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参数,因此也可以写成:

same.markers <- findMarkers(sce.pbmc, groups=colLabels(sce.pbmc))

有了这个参数,就可以根据其他方法的分群结果再去探索不一样的marker基因,而不用修改原始的sce.pbmc结果,只要groups参数满足与sce.pbmc列数一致就可以

怎么解释这个结果呢?--以cluster9为例
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给各个并列第一又排了顺序

# 比如这里看到的,并列第一就不止10个
interesting[1:10,1:4
## DataFrame with 10 rows and 4 columns
##                Top      p.value          FDR summary.logFC
##          <integer>    <numeric>    <numeric>     <numeric>
## 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的所有基因也不是难事

best.set <- interesting[interesting$Top <= 6,]
# 看到Top6其实总共有53个基因
> dim(best.set)
[153 21
# 提取这些基因在与各个cluster比较时的logFC
logFCs <- getMarkerEffects(best.set)
# 总共18个cluster,除去自己还剩17个
> dim(logFCs)
[153 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(-55, length.out=101))

看到其中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"

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
##           <integer>   <numeric>   <numeric>     <numeric>
## 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阈值过滤,其实这些都可以自己手动过滤,多个参数只是更方便一点。缺点就是需要记住这个函数有这个参数

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基因热图

best.set <- interesting.up2[interesting.up2$Top <= 5,]
logFCs <- getMarkerEffects(best.set)

library(pheatmap)
pheatmap(logFCs, breaks=seq(-55, length.out=101))
注意
  • 我们这里只是探索了上调基因。至于有一些亚群中可能部分基因下调才导致这个亚群与众不同,这样的亚群这里检测不到

  • 这里的阈值可能会漏掉一些改变幅度不大,但依然重要的基因

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中是上调的

# 设置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检验

首先来了解一下各种统计检验函数

参考: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方法
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
##           <integer>    <numeric>    <numeric>   <numeric>
## 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方法
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)
image-20200705135623191

2.5 整合多个统计分析的结果

这样可以检验哪些基因是强有力的marker基因(能撑得住三大检验方法的考验)

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])
image-20200705141801833

当然也可以重点关注其中的某一个指标,比如wilcox的结果AUC越大,就表示基因表达量在各个cluster之间的分布越分散;t-test的logFC结果越大,表示更容易解释一个基因在两个cluster之间的变化幅度

3 封锁一些不重要因素

大型数据集一般会涉及许多变化因素(比如批次、性别差异、年龄差异等等),它们会对数据波动产生影响,但又不是感兴趣的生物因素。如果在检测marker基因的时候带着它们,可能会结果产生干扰。

因此有一个参数block可以帮我们”锁住“它们,之后检测的时候就不会把这些因素纳入考量

# 举个例子,在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
##                    <integer>   <numeric>   <numeric>     <numeric>
## 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

这个参数对上面的三种统计方法都适用


点击底部的“阅读原文”,获得更好的阅读体验哦😻


往期回顾

Seurat教程 || 分析Cell Hashing数据

线上讲座 | 单细胞空间转录组专题——实验技术和前沿应用

scRNA-seq Clustering(二)

参考基因组用错了单细胞转录组流程照样可以走通

CellChat作者开讲啦!

Seurat小提琴图为什么有的只有点儿?

单细胞交响乐6-降维

单细胞数据中到底应该如何处理线粒体基因

以复现图表的方式来学习一篇文章

这可能是我见过的最草率的单细胞分群了






如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程



看完记得顺手点个“在看”哦!


生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存