edgeR で発現変動遺伝子を検出する際に利用する RNA-Seq データの正規化法

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),    # add counts from technical replicates
  rowSums(count.liver)      # add counts from technical replicates
)
rownames(counts) <- x[, 1]
counts <- counts[rowSums(cpm(counts) > 0) >= 1, ]
group <- c("Kidney", "Liver")

 
d <- DGEList(counts = counts, group = group)
d <- calcNormFactors(d)
 
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. その他のライブラリー」で正規化係数を算出している。

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