scRNA-seq生信数学技术学习

seurat-FindAllMarkers()源码解析

2021-12-13  本文已影响0人  whitebird

现在回想单细胞项目的一些常规分析，clusters/groups间的差异表达分析（differential expression analysis）几乎是分析必备项，FindAllMarkers()和FindMarkers()的使用频率非常高。做差异表达分析的工具很多，我们如何根据自己的数据和实验处理方案来选择合适的分析工具，还是有些不清晰。

这次想通过梳理FindAllMarkers源码，回答以下几个问题：
1.scRNA-seq数据的哪些方面使差异基因分析复杂化
2.FindAllMarkers算法原理
3.FindAllMarkers使用的注意事项

一、cluster差异表达分析

对cluster进行差异表达分析，使我们能够找到特定于每个cluster的标记基因，也有助于cluster的细胞类型鉴定。
根据cluster的定义，它受cluster本身定义方式的影响，因此在找到标记之前确定正确的cluster分辨率很重要。 如果某些cluster群缺少显著性标记，可以尝试调整聚类数目。

1.1 差异表达分析方法
1.2. seurat执行差异表达分析

根据单细胞数据预处理的方式不同（lognormalize和sctransform），执行差异表达分析代码有所不同。
经sctransform预处理的单细胞数据，进行差异表达分析：

# lognormalize预处理
DefaultAssay(pbmc) <- "RNA" 
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

# sctransform预处理
# method1
DefaultAssay(pbmc) <- "RNA"  
pbmc <- NormalizeData(pbmc)  
all.genes <- rownames(pbmc)  
pbmc <- ScaleData(pbmc, features = all.genes)  
all.markers.RNA <- FindAllMarkers(pbmc, only.pos = T, min.pct = 0.5, logfc.threshold = 0.5) 

# 使用SCT assay进行差异表达分析
# method2
DefaultAssay(srat) <- "SCT"  
all.markers.SCT <- FindAllMarkers(pbmc, only.pos = T, min.pct = 0.5, logfc.threshold = 0.5)
1.3 cluster差异基因列表

Differentially expressed gene (DEG)：差异表达基因
生成的差异基因列表列名说明：

名称 解释
gene 基因名称
cluster 该基因对应的cluster
pct.1 在当前cluster细胞中检测到该基因表达的细胞比例
pct.2 在其它cluster细胞中检测到该基因表达的细胞比例
avg_logFC 两组间平均logFC，Seurat v4默认log2。正值表示特征在第一组中表达得更高
p_val 未调整P-value，数值越小越显著
p_val_adj 基于使用数据集中所有特征的bonferroni校正，校正后的p值
image.png
1.4 为什么scRNA-seq的差异表达分析与bulk RNA-seq不同

由于单细胞测序技术的局限性，单细胞测序数据通常具有高噪音，有较高的dropout问题，即很多低表达或中度表达的基因无法有效检测到，造成基因表达矩阵大量的零值。所以，以前针对传统bulk RNA-seq测序数据开发的差异表达检测方法或软件不一定完全适用于单细胞测序数据。

1.5 常用的统计检验
  1. 非参数检验，如Wilcoxon rank sum test (Mann-Whitney U test)
    例如scRNA-seq 中存在dropouts(零值)的情况，可能会失败
  2. 特定于scRNA-seq的方法，如MAST
    充分利用每组的大量样本（细胞）
    MAST考虑了随机dropouts 丢失和双细胞的情况
  3. 特定于bulk RNA-seq，如DESeq2
    基于负二项分布，适用于UMI数据
    注意：不应该过滤基因，因为 DESeq2 通过借用其他具有相似表达水平的基因的信息来模拟离散度
    运行非常慢！ 仅用于比较2个cluster
算法 工具
Hurdle model combining tests for difference in mean and detection rate MAST Appropriate for read counts or molecule/UMI counts
Non-parametric, based on ranking expression Wilcoxon Robust and very easy to use
Negative binomial model GLM, parametric DESeq2/edgeR
1.6 统计学中的Wilcoxon rank-sum test

统计学中的假设检验：用来判断样本与样本，样本与总体的差异是由抽样误差引起还是本质差别造成的统计推断方法。基本是思路类似数学中的“反证法”，是先对总体做出某种假设，然后通过抽样研究的统计推理，对此假设应该是拒绝或者接受。
在假设检验中，我们设定原假设H0和备选假设H1：

在假设检验中，先设定原假设(H0)，再设定与其相反的备择假设(H1)。接下来随机抽取样本，若在原假设成立的情况下，样本发生的概率(P)非常小，说明原假设不成立，备择假设成立，则拒绝原假设。否则，接受原假设。

判断依据：依据小概率事件的原理，统计学上，把小概率事件在一次实验中看成是实际不可能发生的事件，一般认为等于或小于0.05或0.01的概率为小概率。当P值小于或等于显著性水平时，则拒绝原假设。

Wilcoxon检验（也称为 Mann-Withney-Wilcoxon 检验）：

Wilcoxon秩和检验：推断连续型变量资料或有序变量资料的两个独立样本代表的两个总体分布是否有差别。

  1. 从每个总体中选择随机样本，设 n1 和 n2 分别是较小样本和较大样本中的观察数；

  2. 编秩：将两组数据n1+n2个观测值混合，由小到大统一编秩；不同组遇到相同数据取平均秩次；

  3. 计算各组秩和；

  4. 建立双尾检验假设，确定检验水准

    • H0: 两组样本的总体中位数等于0
    • H1: 两组样本的样本的中位数不等于0
    • α =0.05

  5. 根据分组样本数、各组秩和查询临界值表格，得出统计学结论；

    • 利用统计软件或查Wilcoxon符号秩检验的分布表，

    • 如果P值较小（比如小于或等于给定的显著性水平，如α=0.05），则可以拒绝原假设；如果P值较大，则没有充分的证据来拒绝原假设。


      image.png
image.png
1.7 统计学中的Bonferroni校验

问题：为什么要进行P值校正？
当同一研究问题下进行多次假设检验时，不再符合小概率原理所说的“一次试验”，需要进行多重检验校正。如果在该研究问题下只要有检验是阳性的，就对该问题下阳性结论的话，对该问题的检验的犯一类错误的概率就会增大。如果同一问题下进行n次检验，每次的检验水准为α（每次假阳性概率为α），则n次检验至少出现一次假阳性的概率会比α大。换句话说，在同时对多组数据进行处理和比较的时候，很可能其中部分数据因为随机效应而超过阈值，造成假阳性结果。而检验的次数越多，出现假阳性的概率就越大。

Bonferroni校正：最严格的多重检验校正方法，该方法是由Carlo Emilio Bonferroni发展的，故称Bonferroni校正。
其校正原理为：在同一数据集上同时检验n个相互独立的假设，那么用于每一假设的统计显著水平，应为仅检验一个假设时的显著水平的1/n。举个例子：如要在同一数据集上检验两个独立的假设，显著水平设为常见的0.05。此时用于检验该两个假设应使用更严格的0.025。即0.05* (1/2)。假如有10个检验，经过Bonferroni 校正以后，新的p值可以这样计算1-(1-p)^(1/10)=p'，p'就是 新的显著性水平。

我们要对数千个基因进行统计测验，所以我们可能只是偶然获得p-values很小的基因（假阳性基因），但是实际，在样本间毫无差异。每个基因的p值只是单一检验(单基因)的结果。检验的基因越多，假阳性率就越高。所以，我们需要对基因的p值进行多重假设检验校正，从而减低假阳性结果在我们的检验中出现的次数。
FindAllMarkers函数默认使用Bonferroni校验。校验的次数取决于测试（基因）的数目，如果我们测试较少的基因，校正就没有那么苛刻。
Bonferroni correction: adjusted p-value = raw p-value * number of genes tested

二、FindMarkers函数解析

FindMarkers函数就是一个S3泛型函数；包含Seurat，DimReduc，Assay等数据类型。

# NAMESPACE
S3method(FindMarkers,Assay) 
S3method(FindMarkers,DimReduc) 
S3method(FindMarkers,Seurat) 
S3method(FindMarkers,default)

FindMarkers <- function(object, ...) {  
  UseMethod(generic = 'FindMarkers', object = object)  
}

# differential_expression.R
FindMarkers.default <- function(){...}
FindMarkers.Assay <- function(){...}
FindMarkers.DimReduc <- function(){...}
FindMarkers.Seurat  <- function(){...}

案例1：找出seurat对象seurat_obj中cluster0和cluster1的差异基因。

DefaultAssay(seurat_obj) <- "RNA" 
DEG_wilcox <- FindMarkers(seurat_obj, ident.1 = 0, ident.2= 1,  test.use = "wilcox", min.pct = 0.5, logfc.threshold = 0.5, only.pos = T) 
head(DEG_wilcox) 
#           p_val avg_log2FC pct.1 pct.2 p_val_adj 
# RBP7         0   1.698551 0.726 0.004         0 
# PGD          0   1.670035 0.823 0.046         0 
# AGTRAP       0   1.488062 0.833 0.136         0 
# TNFRSF1B     0   1.696910 0.858 0.119         0 
# EFHD2        0   2.097743 0.917 0.079         0 
# CDA          0   1.577383 0.760 0.001         0 
dim(DEG_wilcox) 
# [1] 789   5

2.1 在RNA assay上做差异表达分析

LogNormalize预处理后进行差异基因分析，FindMarkers函数的代码很多，我把主要步骤记录如下：
step1：将以上FindMarkers命令中的参数传至FindMarkers.Seurat()函数；

object <- seurat_obj 
ident.1 = 0 
ident.2= 1 
test.use = "wilcox" 
min.pct = 0.5 
logfc.threshold = 0.5 
only.pos = T

step2：获取归一化的基因表达矩阵seurat_obj[['RNA']]@data；在FindMarkers.Assay()函数会调用FoldChange()函数。

data.use <-  GetAssayData(object = object, slot = "data") 
fc.results <- FoldChange( 
  object = object, 
  slot = data.slot, 
  cells.1 = cells.1, 
  cells.2 = cells.2, 
  features = features, 
  pseudocount.use = pseudocount.use, 
  mean.fxn = mean.fxn, 
  fc.name = fc.name, 
  base = base 
)

step3：FoldChange.default()的计算非常关键。fold change是组间基因表达量的差异倍数，avg_log2FC是基因在组间的平均表达量取log2，Seurat v4默认log2。正值表示特征在第一组中表达得更高。
FoldChange计算方式如下：

# 基因表达均值计算方式
mean.fxn <- mean.fxn %||% switch( 
  EXPR = slot, 
  'data' = function(x) { 
    return(log(x = rowMeans(x = expm1(x = x)) + pseudocount.use, base = base)) 
  }, 
  'scale.data' = rowMeans, 
  function(x) { 
    return(log(x = rowMeans(x = x) + pseudocount.use, base = base)) 
  } 
)

step4：每个基因计算在两个细胞群的占比 pct.1和 pct.2，计算方式是：基因阳性表达细胞数/该群细胞总数；

pct.1 <- round(
  x = rowSums(x = object[features, cells.1, drop = FALSE] > thresh.min) /
    length(x = cells.1),
  digits = 3
)
pct.2 <- round(
  x = rowSums(x = object[features, cells.2, drop = FALSE] > thresh.min) /
    length(x = cells.2),
  digits = 3
)

每个基因的avg_log2FC的计算：
cells.1群细胞的每个基因表达量取expm1，求平均+1，再取log2；
cells.2群细胞的每个基因表达量取expm1，求平均+1，再取log2；
每个基因的avg_log2FC计算：fc <- (data.1 - data.2)

data.1 <- mean.fxn(object[features, cells.1, drop = FALSE])
data.2 <- mean.fxn(object[features, cells.2, drop = FALSE])
mean.fxn <- mean.fxn %||% switch(
  EXPR = slot,
  'data' = function(x) {
    return(log(x = rowMeans(x = expm1(x = x)) + pseudocount.use, base = base))
  },
  'scale.data' = rowMeans,
  function(x) {
    return(log(x = rowMeans(x = x) + pseudocount.use, base = base))
  }
)
fc <- (data.1 - data.2)
image.png

step5：我们进入到FindMarkers.default()这个函数，先判断cluster0和cluster1的细胞数是不是满足最小细胞数3，min.cells.group默认是3个，cluster0和cluster1最少要3个以上细胞；显然，符合要求。根据设置的最小百分比min.pct 和 log2FC阈值和仅保留阳性基因，对所有的36601基因进行过滤，仅保留789个基因；
step6：调用PerformDE()函数，执行wilcox检验和P值校正；
p_val的计算：test.use = "wilcox"，采用WilcoxDETest检测计算出p_val；

de.results <- switch(
  EXPR = test.use,
  'wilcox' = WilcoxDETest(
    data.use = data.use,
    cells.1 = cells.1,
    cells.2 = cells.2,
    verbose = verbose
  )
)

p_val_adj的计算：对p_val进行Bonferroni校正

de.results$p_val_adj = p.adjust(
      p = de.results$p_val,
      method = "bonferroni",
      n = nrow(x = object)
    )
image.png

汇总FindMarkers函数的主要步骤如下：


image.png

2.2 在SCT assay上做差异表达分析

经过SCTransform预处理后进行差异基因分析，看下哪些步骤和上面的不同：

DefaultAssay(seurat_obj) <- "SCT"  
DEG_wilcox <- FindMarkers(seurat_obj, ident.1 = 0, ident.2= 1,  test.use = "wilcox", min.pct = 0.5, logfc.threshold = 0.5, only.pos = T)
head(DEG_wilcox)
# p_val avg_log2FC pct.1 pct.2 p_val_adj
# RBP7         0   1.475970 0.709 0.004         0
# PGD          0   1.464552 0.804 0.045         0
# AGTRAP       0   1.296850 0.803 0.135         0
# TNFRSF1B     0   1.473898 0.830 0.120         0
# EFHD2        0   1.821583 0.915 0.086         0
# CDA          0   1.373028 0.745 0.001         0
dim(DEG_wilcox)
# [1] 686   5
image.png

阅读完源码，RNA/SCT assay做差异表达分析只有输入的归一化矩阵不同，fold change计算方法完全一致。

小心得：

  1. 求基因在细胞群的平均表达量
    通过阅读源码，seurat在计算基因在细胞群的平均表达量不是对细胞群的每个细胞归一化的表达值求均值；而是对将细胞群的细胞归一化的基因表达量进行expm1还原，再求均值。这也是合理的，AverageExpression()函数在求基因的平均表达也是如此处理的。
    这个怎么理解呢？
    可以把细胞群看成一个大细胞，包含每个细胞校正后的counts，我们求单个细胞的基因表达量=基因的总counts/细胞内总counts，另外，每个细胞考虑到测序深度的影响会乘以scale_factor进行校正；
    计算步骤：
  1. avg_log2FC的计算梳理
    前面我们说了基因的平均表达量计算：
    cells.1群细胞的每个基因表达量取expm1，求平均；
    cells.2群细胞的每个基因表达量取expm1，求平均；
    但是我们是计算基因的avg_log2FC；所以要取log2；
    那为什么会出现“均值+1”，这个怎么理解，计算方式如下：

每个基因的avg_log2FC的计算：
cells.1群细胞的每个基因表达量取expm1，求平均+1，再取log2；
cells.2群细胞的每个基因表达量取expm1，求平均+1，再取log2；
每个基因的avg_log2FC计算：fc <- (data.1 - data.2)

答：这是避免log2(基因的平均表达值)为负，
如果基因的平均表达值=0，log2(基因的平均表达值)会负；
但是log2(基因的平均表达值+1)=0，避免计算结果为负。

细胞群的平均表达量的方式：total sum scaling (TSS) normalization
The expm1 does un-log the data, but the normalization persists (this would be lost in the counts slot)
expm1() transformed in order to recover normalized values not in log scale
在网上也请教过一个大佬，以下是他的答复，非常感谢迷茫时有人点灯

你要取平均，很自然的应该用原始数据啊（总共三个苹果，五个小孩，每个小孩0.6个苹果）。log以后的数据取平均，你说这是个什么东西表达了什么含义呢，就很难让人理解的。（不过你非要log以后取平均，我猜测也能得到相似的结果，但是很难让人理解这东西到底是个啥……）

上一篇下一篇

猜你喜欢

热点阅读