DEGESに基づく正規化

DEGES とは differentially expressed gene elimination strategy のことで、RNA-seq のリードカウントデータを正規化する戦略手法です。TMM 正規化法はハウスキーピング遺伝子の log-fold-change をゼロに近づけるという点に着目していると同様に、DEGESに基づく正規化手法もこれに着目しています。

TMM 正規化法と DEGES-based 正規化法の違い

TMM 正規化法では log-fold-change の値(M値)を計算して、上位から30%、下位から30%ずつ遺伝子を取り除いています。平均発現量(A値)についても、上位5%、下位5%の遺伝子を取り除いています。次に、残った遺伝子から正規化係数を計算しています。TMM 正規化法では上位と下位から一定量ずつのデータを削除するのは、二つのサンプル間で発現量が変動している遺伝子を取り除くためです。これによって、残ったデータはハウスキーピング遺伝子などの非発現変動遺伝子となり、正規化係数の計算に用いられます。

一方、DEGES に基づく正規化法では二つのサンプル間で発現量が変動している遺伝子を取り除くために、別のアプローチをとっています。具体的に次のような手順になっています。

  1. オリジナルのリードカウントデータを TMM 正規化します。
  2. TMM 正規化されたデータに対して検定を行い、FDR < 0.1 を満たす遺伝子を選び、削除します。
  3. 残ったデータを利用して TMM 正規化法で正規化係数を計算します。

TMM 正規化法では上位と下位から一定量ずつの遺伝子を取り除くことに対して、DEGES に基づく正規化手法では、 FDR < 0.1 となる遺伝子を取り除くことにしています。これにより、二つのサンプル間で発現量が変動している遺伝子をより正確に除去でき、データをより正確に正規化できます。

TMM 正規化と DEGES-based 正規化の比較

次に示す二つのシナリオでシミュレーションデータを生成して、TMM 正規化と DEGES-based 正規化をデータに対して施します。

  1. 全遺伝子数が10,000個、そのうち2,000個が発現変動遺伝子であり、発現変動遺伝子のうち1,000個がG1群で高発現し、残りの1,000個がG2群で高発現している状況
  2. 全遺伝子数が10,000個、そのうち2,000個が発現変動遺伝子であり、発現変動遺伝子のうち1,600個がG1群で高発現し、残りの400個がG2群で高発現している状況

正規化後のデータを M-A plot に示すと、次のようになります。灰色の点は非発現変動遺伝子、マゼンタ色の点は発現変動遺伝子を表します。また、オレンジ色の線は、非発現変動遺伝子の y 座標(log-fold-change)の median です。この線が y = 0 に近づくほどよいです。状況1では TMM 正規化法と DEGES-based 正規化法が同じ結果となりましたが、状況2の場合は DEGES-based 正規化法が TMM 正規化法よりも優れいていることが確認できます。

状況1状況2
正規化されていない状況1データ 正規化されていない状況2データ
TMM正規化された状況1データ TMM正規化された状況2データ
DEGES正規化された状況1データ DEGES正規化された状況2データ

DEGES-based 正規化の実行方法

DEGES-based 正規化法は R の TCC パッケージに実装されています。このパッケージを利用します。(サンプルデータ

library(TCC)

data <- read.table("counts33.txt", header = TRUE)
group <- c("A", "A", "A", "B", "B", "B")

# 準備
count <- data[, -1]
rownames(count) <- data[, 1]
tcc <- new("TCC", count, group)

# DEGESに基づく正規化で正規化係数を計算
tcc <- calcNormFactors(tcc)
tcc$norm.factors                 # 計算された正規化係数
##        A1        A2        A3        B1        B2        B3 
## 0.9820031 0.9509060 1.0308594 1.0242789 1.0061235 1.0058292 

## DEGESに基づく正規化後のデータを取得と表示
data.normalized <- getNormalizedData(tcc)
head(data.normalized)
##               A1          A2        A3        B1        B2        B3
## gene_1  4.018350   1.9928582  1.994290  0.000000  0.000000  0.000000
## gene_2  1.004588   0.9964291  0.000000  8.938687  3.006814  1.006516
## gene_3 33.151390  53.8071709 74.785866 16.884186 13.029526 12.078189
## gene_4  0.000000   0.0000000  3.988580  0.000000  1.002271  0.000000
## gene_5 24.110102 145.4786471 53.845823 22.843311 15.034069 22.143347
## gene_6  3.013763   1.9928582  5.982869  0.000000  0.000000  0.000000

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
  • Sun J, Nishiyama T, Shimizu K, Kadota K. TCC: an R package for comparing tag count data with robust normalization strategies. BMC Bioinformatics. 2013, 14:219. PubMed Abstract