TMM正規化

TMM 正規化について

Trimmed mean of M values (TMM) 正規化は、RNA-seq のリードカウントデータを正規化する方法の一つである。細胞内で発現している遺伝子は、ハウスキーピング遺伝子などのような発現変動のない遺伝子(非発現変動遺伝子)が多いことに着目した正規化法である。

例えば 2 群間比較を考えたとき、非発現変動遺伝子の発現量は、2 つの実験群においてその対数比(M 値)はゼロになると期待される。しかし、何らかの技術的な原因などにより、その M 値の期待値がゼロとかけ離れる場合がある。TMM 正規化法は、非発現変動遺伝子の M 値の期待値をゼロに持っていくような係数を計算している。この係数のことを正規化係数(normalization factor)という。非発現変動遺伝子の M 値の期待値がゼロでないデータに、この正規化係数を作用させることで、その期待値が(理論上)ゼロになる。

例えば、以下の腎臓と肝臓 (Marioni et al.) の M-A plot でみられるように、正規化前ではハウスキーピング遺伝子M 値の mean が -0.845 であるのに対して、TMM 正規化後は -0.249 となった。プロットを眺めると、すべての点が y 軸方向に +0.596 だけシフトしていることがわかる。正規化係数はこのように y 軸(= M 値 = 対数比)におけるシフトの程度を意味しています。

正規化前のmarioniのデータのMAplot TMM正規化後のmarioniのデータのMAplot

TMM 正規化法

TMM 正規化(正規化係数の計算方法)は以下の手順で行われる。

  1. すべての遺伝子について M 値と A 値を計算する。
  2. すべての M 値について、上位 30% と下位 30% を除去する。
  3. すべての A 値について、上位 5% と下位 5% を除去する。
  4. 残ったデータを利用して正規化係数を計算する。

M 値と A 値に基いてトリムしたあとに残ったデータは、下図のマゼンダ色の点で表している。正規化係数はこれらのマゼンダ色から計算される。

TMM正規化法のトリーム領域

マゼンタ色の遺伝子の mean が -0.621 であるから、M-A プロットでいうと、これがゼロとなるように y 軸方向の +0.621 だけシフトさせるような係数を計算しています。

edgeR パッケージで正規化係数を求める例。利用するデータは Marioni et al. の肝臓と腎臓のデータ(SupplementalTable2.txt)。ダウンロードしたデータをそのまま R で読み込んで実行する。

library(edgeR)

# 解析用のデータの読み込みと整形
f <- read.table("humanhousekeeping.txt", header = TRUE)
x <- read.table("SupplementaryTable2.txt", header = TRUE)
count.kidney <- x[, grep("Kidney", colnames(x))]
count.liver  <- x[, grep("Liver", colnames(x))]
counts <- cbind(
  rowSums(count.kidney),    # technical replicates のためリードカウントをすべて足す
  rowSums(count.liver)      # technical replicates のためリードカウントをすべて足す
)
rownames(counts) <- x[, 1]
counts <- counts[rowSums(cpm(counts) > 0) >= 1, ]
group <- c("Kidney", "Liver")

 
# edgeRを利用して正規化係数を求めます
d <- DGEList(counts = counts, group = group)
d <- calcNormFactors(d)
 
# TMM正規化係数
d$samples$norm.factors
## [1] 1.2297972 0.8131422
d$samples$norm.factors / d$samples$norm.factors[1]
## [1] 1.0000000 0.6477503

この結果から、TMM 正規化係数をデータ全体に掛けることで、データ全体が y 軸方向に + 0.648 だけ移動することになります。これは、M-A プロット上のマゼンタ色の点の mean をゼロにシフトさせることを意味しています。(上で見た M-A プロットでみた数値と実際に異なっているのはなぜだろう???)

M値およびA値の求め方

遺伝子 g = 1, 2, 3, ..., G について、M 値及び A 値は次のようにして計算する。

\[ \begin{eqnarray} M_{g} &=& \log_{2}\frac{\left( \frac{Y_{gk}}{N_{k}}\right)}{\left(\frac{Y_{gl}}{N_{l}}\right)} \\ A_{g} &=& \frac{\log_{2}\left( \frac{Y_{gk}}{N_{k}}\times\frac{Y_{gl}}{N_{l}} \right)}{2} \end{eqnarray} \]

ここで、Ygkはライブラリー k の遺伝子 g のリードカウントデータを表す。また、Nk はライブラリー k のライブラリーサイズです(Nk = Σg=1GYgk)。

TMM 正規化係数の求め方

M 値および A 値に基いてトリムしたあとに残ったデータを利用して計算する。ライブラリー k の正規化係数はライブラリー r を対照群として次のように計算する。

\[ \log_{2}(TMM_{k}^{r}) = \frac{\sum_{g=1}^{G}w_{gk}^{r}M_{gk}^{r}}{\sum_{g=1}^{G}w_{gk}^{r}} \] \[ w_{gk}^{r} = \frac{N_{k} - Y_{gk}}{N_{k}Y_{gk}} + \frac{N_{r}-Y_{gr}}{N_{r}Y_{gr}} \]

TMM 正規化と多群間比較について

TMM 正規化はハウスキーピング遺伝子などの非発現変動遺伝子の log-fold-change をゼロにシフトさせる方法である。二つのライブラリーを比較して、そのシフトする程度となる TMM 正規化係数を計算する。二つのライブラリーというのは、「肝臓 vs. 腎臓」を意味しているだけではなく、「肝臓1 vs 肝臓2」でも構わない。

2群間の場合、肝臓サンプルが三つの biological replicate があり、腎臓も同様に 3 biological replicate あると仮定した時、TMM 正規化係数は次のように計算される。

  1. 6 biological replicate の中から、reference library を決める。(確か upper quartile の最も小さいライブラリーをリファレンスとしている)
  2. 「リファレンスライブラリー vs. 残りの 5 ライブラリー」で 5 つの正規化係数を算出します。この時点ではリファレンスの正規化係数は 1 で、他は 0.xx や 1.xx などようになる。
  3. 最後に、正規化係数の平均を揃えたりなどして、プログラムから出力されるような、すべての係数が小数値となるような正規化係数になる。

この操作からして、正規化係数を計算するとき、肝臓か腎臓かをそもそも配慮していない。あくまで「リファレンスライブラリー vs. その他のライブラリー」である。

このことから、「肝臓 vs. 腎臓 vs. 心臓」などの多群間比較のデータに対しても、結局、リファレンスライブラリーを決めて、「リファレンスライブラリー vs. その他のライブラリー」というように正規化係数を計算することになる。

つまり、二群間であろうか、多群間であろうか、Robinson らはそもそもそれを気にしていなく、あくまで「リファレンスライブラリー vs. その他のライブラリー」で正規化係数を算出している。

プロットの作成スクリプト

上の M-A プロットは次のコードで作成した。必要なデータは二つで、どちらもダウンロードしてそのまま R に読み込ませる。

library(TCC)
 
# TMM トリム領域を描く関数
.trimmedGenes <- function(x) {
  r <- x[, 1]
  o <- x[, 2]
  allzero <- (rowSums(x) == 0)
  r <- r[!allzero]
  o <- o[!allzero]
  nO <- sum(o)
  nR <- sum(r)
  logR <- log2((o / nO) / (r / nR))
  absE <- (log(o / nO) + log2(r / nR)) / 2
  fin <- is.finite(logR) & is.finite(absE) & (absE > -1e10)
  logR <- logR[fin]
  absE <- absE[fin]
  n <- sum(fin)
  loL <- floor(n * 0.3) + 1
  hiL <- n + 1 - loL
  loS <- floor(n * 0.05) + 1
  hiS <- n + 1 - loS
  keep <- (rank(logR) >= loL & rank(logR) <= hiL) &
          (rank(absE) >= loS & rank(absE) <= hiS)
  keep.gene <- names(keep)[keep]
  tag <- rep(0, length = nrow(x))
  names(tag) <- rownames(x)
  tag[keep.gene] <- 1
  return(tag)
}
 


# 解析用のデータの読み込みと整形
f <- read.table("humanhousekeeping.txt", header = TRUE)
x <- read.table("SupplementaryTable2.txt", header = TRUE)
count.kidney <- x[, grep("Kidney", colnames(x))]
count.liver  <- x[, grep("Liver", colnames(x))]
counts <- cbind(
  rowSums(count.kidney),    # technical replicates のためリードカウントをすべて足す
  rowSums(count.liver)      # technical replicates のためリードカウントをすべて足す
)
rownames(counts) <- x[, 1]
counts <- counts[rowSums(cpm(counts) > 0) >= 1, ]
group <- c("Kidney", "Liver")
tcc <- new("TCC", counts, group)

 
# housekeeping genes を緑色(2)で描くように設定する
col.tag <- rep(1, length = nrow(counts))
names(col.tag) <- rownames(counts)
inter <- intersect(names(col.tag), as.character(f[, 1]))
col.tag[inter] <- 2
 

# 正規化前の M-A プロットを描く
xlab <- expression(A == (log[2] * Kidney + log[2] * Liver)/2)
ylab <- expression(M == log[2] * Liver - log[2] *  Kidney)
png("marioni.org.png", 320, 340)
par(mar = c(4.5, 4.5, 1.5, 0.5))
xy <- plot(tcc, col = c("darkgrey", "green"), col.tag = col.tag,
           ylim = c(-8, 8), xlab = xlab, ylab = ylab, main = "un-normalized")
m <- mean(xy$m.value[col.tag == 2])
abline(h = m, col = "orange", lwd = 2)
text(13, m + 0.5, sprintf("%.3f", m), cex = 1.4, offset = 0, pos = 3)
dev.off()
 

# 正規化後の M-A プロットを描く
png("marioni.tmm.png", 320, 340)
par(mar = c(4.5, 4.5, 1.5, 0.5))
tcc <- calcNormFactors(tcc, norm.method = "tmm", iteration = F) # TMM normalization
xy <- plot(tcc, col = c("darkgrey", "green"), col.tag = col.tag,
           ylim = c(-8, 8), xlab = xlab, ylab = ylab, main = "TMM normalized")
m <- mean(xy$m.value[col.tag == 2])
abline(h = m, col = "orange", lwd = 2)
text(13, m + 0.5, sprintf("%.3f", m), cex = 1.4, offset = 0, pos = 3)
dev.off()
 

# TMM 正規化法トリム領域を描く
tcc <- new("TCC", counts, group)
trimmed.genes <- .trimmedGenes(tcc$count)
png("marioni.trimarea.png", 320, 340)
par(mar = c(4.5, 4.5, 1.5, 0.5))
xy <- plot(tcc, col = c("darkgrey", "magenta"), col.tag = trimmed.genes + 1,
           ylim = c(-8, 8), xlab = xlab, ylab = ylab, main = "")
m <- mean(xy$m.value[trimmed.genes == 1])
abline(h = m, col = "orange", lwd = 2)
text(13, m + 0.5, sprintf("%.3f", m), cex = 1.4, offset = 0, pos = 3)
dev.off()

References

  • Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010, 11(3):R25. PubMed Abstract
  • Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y. RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Res. 2008, 18(9):1509-17. PubMed Abstract
  • Eisenberg E, Levanon EY. Human housekeeping genes are compact. Trends Genet. 2003, 19(7):362-95. PubMed Abstract