理论-edgeR-TMM-normalization

最近在看差异分析当中原始read counts是如何被校正的，自然就不会放过差异分析的经典之一 —— edgeR.

edgeR使用的校正方法称为trimmed mean of M values (TMM)，其前提假设为样本对照组和处理组间绝大多数基因表达不发生差异。

如何界定绝大多数基因这一点我个人还没有看到一个量化的指标，是50%还是80%才算绝大多数。

edgeR的TMM校正方法其实在其分析流程中就是一句命令而已

1	y <- calcNormFactors(y)

但由于我十分好奇其背后计算的原理和计算过程，我便搜索一番，才发现无论是中文还是英文的帖子对于TMM的具体运算步骤及代码都没有很好的整理。因此，本文记录我通过阅读TMM提出的原文和edgeR源码所了解到的TMM校正。

TMM校正示例

我们通过airway数据展示实际TMM校正的过程

# TMM normalization
data("airway", package='airway')

set.seed(123)
counts1 <- as.data.frame(assay(airway)[sample(1:nrow(airway), 5000),1:4])

Step 1: 选择参考样本

TMM normalization首先需要选择一个参考样本，以它为基准进行校正。

默认下，参考样本的选择是通过比较每个样本的CPM (counts per million)的上四分位数与所有样本CPM的平均上四分位数之间的差值，找出差值最小的样本作为参考样本。

# library size of each samples
lib.size <- colSums(counts1)

# library size normalized read count for each gene in each sample
counts1.cpm <- t(t(counts1)/lib.size) * 1e6

# reference library for TMM normalization
f75 <- apply(counts1.cpm, 2, function(x) quantile(x, probs = 0.75))
# select the one with upper quartile closest to mean upper quartile
refColumn <- which.min(abs(f75 - mean(f75)))
> refColumn
SRR1039512 
         3

在这个例子中refColumn是SRR1039512

在edgeR::calcNormFactors()中，我们也可以通过refColumn=参数指定参考样本

Step 2: calculation of sample-reference pairwise M and A

接着，在参考样本和非参考样本间两两计算校正因子（normalization factors）。

我们首先需要计算参考样本和非参考样本间的Fold change (M)和平均表达量 (A)

$M = log2\frac{non-reference\ sample\ count}{reference\ sample\ count}$

$A=\frac{log2(non-reference\ sample\ count)\ +\ log2(reference\ sample\ count)}{2}$

注意这里使用的count都是校正过文库大小的

以下代码用到的变量名含义：

obs: 非参考样本的原始count

ref: 参考样本的原始count

libsize.obs: 非参考样本的原始文库大小

libsize.ref: 参考样本的原始文库大小

# 第一列是非参考样本
obs <- as.numeric(counts1[, 1])
ref <- as.numeric(counts1[, refColumn])
libsize.obs <- lib.size[1]
libsize.ref <- lib.size[refColumn]

分别计算M value和A value，为了对M value加权，我们还需要通过delta method估计渐近方差

# M value: log ratio of expression, accounting for library size
M <- log2((obs/libsize.obs)/(ref/libsize.ref))
# A value:absolute expression
A <- (log2(obs/libsize.obs) + log2(ref/libsize.ref))/2
# estimated asymptotic variance
v <- (libsize.obs - obs)/libsize.obs/obs + (libsize.ref - ref)/libsize.ref/ref

保留M和A值均为有限值的基因，并过滤极低表达量的基因

Acutoff = -1e10
#	remove infinite values, cutoff based on A
fin <- is.finite(M) & is.finite(A) & (A > Acutoff)
M <- M[fin]
A <- A[fin]
v <- v[fin]

Step 3: trimmed mean of M values

接着，我们对M和A值进行双重截值，截掉M值排在前30%和后30%，A值排在前5%和后5%的基因，计算中间这部分基因M值的加权平均值

logratioTrim <- 0.3
sumTrim <- 0.05
# Double trim the upper and lower percentages of the data
# trim M values by 30% and A values by 5%
  
  n <- length(M)
  loM <- floor(n * logratioTrim) + 1
  hiM <- n + 1 - loM
  loA <- floor(n * sumTrim) + 1
  hiA <- n + 1 - loA
  
  keep <- (rank(M)>=loM & rank(M)<=hiM) & (rank(A)>=loA & rank(A)<=hiA)
  
  # Weighted mean of M after trimming 
  f <- sum(M[keep]/v[keep], na.rm=TRUE) / sum(1/v[keep], na.rm=TRUE)
  f <- 2^f
  #	Factors should multiple to one
  f <- f/exp(mean(log(f)))

  #	Output
  names(f) <- colnames(counts1)

上述步骤是计算第一个样本与参考样本的校正因子，下面我们通过一个循环计算所有样本的校正因子.

以下循环在选取refColumn后开始

nsamples <- ncol(counts1)

logratioTrim <- 0.3
sumTrim <- 0.05
Acutoff = -1e10

f <- rep_len(NA_real_, nsamples)
for (i in 1:nsamples) {
  obs <- as.numeric(counts1[, i])
  ref <- as.numeric(counts1[, refColumn])
  libsize.obs <- lib.size[i]
  libsize.ref <- lib.size[refColumn]
  
  # M value: log ratio of expression, accounting for library size
  M <- log2((obs/libsize.obs)/(ref/libsize.ref))
  # A value:absolute expression
  A <- (log2(obs/libsize.obs) + log2(ref/libsize.ref))/2
  # estimated asymptotic variance
  v <- (libsize.obs - obs)/libsize.obs/obs + (libsize.ref - ref)/libsize.ref/ref
  #	remove infinite values, cutoff based on A
  fin <- is.finite(M) & is.finite(A) & (A > Acutoff)
  
  M <- M[fin]
  A <- A[fin]
  v <- v[fin]
  
  # Double trim the upper and lower percentages of the data
  # trim M values by 30% and A values by 5%
  
  n <- length(M)
  loM <- floor(n * logratioTrim) + 1
  hiM <- n + 1 - loM
  loA <- floor(n * sumTrim) + 1
  hiA <- n + 1 - loA
  
  keep <- (rank(M)>=loM & rank(M)<=hiM) & (rank(A)>=loA & rank(A)<=hiA)
  
  # Weighted mean of M after trimming 
  f[i] <- sum(M[keep]/v[keep], na.rm=TRUE) / sum(1/v[keep], na.rm=TRUE)
  f[i] <- 2^f[i]
}

#	Factors should multiple to one
f <- f/exp(mean(log(f)))

#	Output
names(f) <- colnames(counts1)

f
#SRR1039508 SRR1039509 SRR1039512 SRR1039513 
# 1.0225375  1.0111616  0.9920515  0.9749133

我们与edgeR::calcNormFactors()比较一下

library(edgeR)
y <- DGEList(counts = counts1)
y <- calcNormFactors(y)
nf <- y$samples$norm.factors

f == nf
#SRR1039508 SRR1039509 SRR1039512 SRR1039513 
#      TRUE       TRUE       TRUE       TRUE

Step 4: TMM normalized counts

最后，我们获取TMM校正后的read counts。实际上，上述计算的校正因子是对文库大小的校正，edgeR再利用校正后的文库大小对read counts进行校正。

${TMM\ Normalized\ Counts = \frac{Raw \ Count\ *\ 10^6}{Library\ size\ *\ NormFactor}}$

# TMM normalized counts
counts1.tmm <- t(t(counts1) / (lib.size * f)) * 1e6
counts1.tmm <- round(counts1.tmm, 4)

head(counts1.tmm)

	SRR1039508	SRR1039509	SRR1039512	SRR1039513
ENSG00000260166	0.0000	0.0000	0.000	0.000
ENSG00000266931	0.0000	0.0000	0.000	0.000
ENSG00000104774	1375.7505	1510.2195	1431.252	1301.323
ENSG00000267583	0.0000	0.7924	0.000	0.000
ENSG00000227581	0.7095	0.0000	0.000	0.000
ENSG00000227317	0.0000	0.0000	0.000	0.000

使用edgeR::calcNormFactors()计算校正因子后，我们可以使用edgeR::cpm()提取校正后的counts矩阵

# by edgeR
counts.tmm.edger <- round(cpm(y), 4)

identical(counts1.tmm, counts.tmm.edger)
# [1] TRUE

以上就是TMM校正的计算过程。

完。

Ref:

A scaling normalization method for differential expression analysis of RNA-seq data: https://doi.org/10.1186/gb-2010-11-3-r25

calcNormFactors.R: https://rdrr.io/bioc/edgeR/src/R/calcNormFactors.R

Gene expression units explained: RPM, RPKM, FPKM, TPM, DESeq, TMM, SCnorm, GeTMM, and ComBat-Seq: https://www.reneshbedre.com/blog/expression_units.html

edgeR提供的TMM归一化算法详解: https://cloud.tencent.com/developer/article/1625225