RNA-seq解析の発現変動遺伝子を検出する方法として、RパッケージのedgeRやDEseqが有名です。その際に、条件間でのFold Change、つまり発現差を計算するのですが、リードカウントが0の場合は0で割ることはできないのでFold Changeを算出することができません。
今回はedgeR内でどのような処理を行い、Fold Changeを算出可能なのか、備忘録としてまとめたいと思います。
結論だけ知りたい方は一番下のまとめをご覧ください。
シミュレーションデータについて
シミュレーションデータはRのTCCパッケージで作成したものを用いました。
データの作成方法については以下の記事から
edgeRでの遺伝子発現変動比較解析の流れ
まずは、edgeRで発現変動比較解析の流れを説明します。
edgeRについては、以下の文献を参照してください。
DOI: 10.1093/bioinformatics/btp616
今回はbiological replicateを１として2群間比較を行います。
以下にログを示します。
> #パッケージの読み込み
> library(edgeR)
> #パッケージバージョンの確認
> packageVersion("edgeR")
[1] ‘3.32.1’
> #カウントデータの読み込み
> count <- read.csv("TCC_simulation_data.csv", header = T, row.names = 1)
> #カウントデータの抽出
> count <- count[, c(3, 6)]
> head(count)
G1_rep3 G2_rep3
gene_1 32 10
gene_2 28 0
gene_3 158 74
gene_4 68 10
gene_5 235 115
gene_6 116 47
> #行列型への変換
> count <- as.matrix(count)
> dim(count)
[1] 10000 2
> #実験群の識別ラベルを作成し、group変数に代入
> group <- factor(c("condition1", "condition2"))
> #カウントデータと実験群の識別ラベルを一つの変数に代入
> d <- DGEList(counts = count, group = group)
> #TMM正規化
> d <- calcNormFactors(d)
> d
An object of class "DGEList"
$counts
G1_rep3 G2_rep3
gene_1 32 10
gene_2 28 0
gene_3 158 74
gene_4 68 10
gene_5 235 115
9995 more rows ...
$samples
group lib.size norm.factors
G1_rep3 condition1 1507360 0.851503
G2_rep3 condition2 982924 1.174394
> #edgeRによる2群間比較
> d <- estimateGLMCommonDisp(d, method = "deviance", robust = T, subset = NULL)
> #exact test
> result <- exactTest(d)
> result
An object of class "DGEExact"
$table
logFC logCPM PValue
gene_1 -1.5139769 4.207542 0.39146424
gene_2 -7.7396297 3.657546 0.01129089
gene_3 -0.9401836 6.573969 0.53604662
gene_4 -2.5983076 5.023355 0.12333235
gene_5 -0.8773113 7.160117 0.55946082
9995 more rows ...
$comparison
[1] "condition1" "condition2"
$genes
NULL
長くなりましたが、解析はこのような流れになります。
ここで注目していただきたいのは、G2_rep3のgene_2です。リードカウントは0になっていますが、excat testの結果の後の「result」では、logFC (Fold ChangeをLog表記したものです)が-7.7396・・・となっています。
単純に考えると不思議な結果です。
G1_rep3のgene_2のリードカウントは28、G2_rep3のgene_2は0ですから、0/28となり、計算できないはずです。
なぜ、このような現象が生じるのでしょうか。
TMM正規化で処理がされている？
TMM正規化とは、非発現変動遺伝子はMAプロット上でM値が0となると仮定し、そのようになるように係数を掛ける手法です(詳しくは以下のサイトを参照)。
つまり、TMM正規化ではリードカウントに係数を掛ける処理しかしないので、カウントが0であれば0のままです。
TMM正規化をするだけではFold Changeは計算できないままのようです。
edgeRによる2群間比較の際に処理されている？
edgeRで使用した関数は主に以下の二つです。
estimateGLMCommonDisp(d, method = "deviance", robust = T, subset = NULL)
exactTest(d)
このどちらかで何らかの処理がされていると考えられます。
まずはestimateGLMCommonDisp関数から調べました。
以下をRstudioのコンソールに記入します。
?estimateGLMCommonDisp
そうすると、help画面が起動し、estimateGLMCommonDisp関数のArgumentsなどを見ることができます。しかし、helpを見ても特に関連しそうな項目は見つかりません。
次に、exactTest関数について調べました。
?exactTest
すると、関連しそうな項目がありました(以下はRのhelp画面からの引用です)。
Usage
exactTest(object, pair=1:2, dispersion="auto", rejection.region="doubletail",
big.count=900, prior.count=0.125)
Arguments
prior.count
average prior count used to shrink log-fold-changes. Larger values produce more shrinkage.
Argumentsのprior.countが原因だったようです。
これにより、prior.countで指定した数字をTMM正規化後の数値に足すことで、リードカウントが0であるものを無くし、Fold Changeが正常に計算できるようになっています。
デフォルトは0.125のようです。
試しにprior.countを0に指定してみます。
> #exact test
> result <- exactTest(d, prior.count=0)
> result
An object of class "DGEExact"
$table
logFC logCPM PValue
gene_1 -1.5250327 4.207542 0.39146424
gene_2 -Inf 3.657546 0.01129089
gene_3 -0.9412882 6.573969 0.53604662
gene_4 -2.6124956 5.023355 0.12333235
gene_5 -0.8779877 7.160117 0.55946082
9995 more rows ...
$comparison
[1] "condition1" "condition2"
$genes
NULL
gene_2のlogFCが-Infとなっています。前述したようにG1_rep3のgene_2のリードカウントは28、G2_rep3のgene_2は0ですから、0/28となり、-Infとなるはずなので、予想された結果が得られたということになります。
デフォルト値について、usageに書いてある情報しか見つからなかったので、一応デフォルトが0.125なのかどうかも調べてみましょう。
> #exact test
> result <- exactTest(d, prior.count=0.125)
> result
An object of class "DGEExact"
$table
logFC logCPM PValue
gene_1 -1.5139769 4.207542 0.39146424
gene_2 -7.7396297 3.657546 0.01129089
gene_3 -0.9401836 6.573969 0.53604662
gene_4 -2.5983076 5.023355 0.12333235
gene_5 -0.8773113 7.160117 0.55946082
9995 more rows ...
$comparison
[1] "condition1" "condition2"
$genes
NULL
少しわかりにくいので以下に、デフォルト設定で行ったもの、0.125で指定したものを上下に並べてみました。
デフォルト設定
logFC logCPM PValue
gene_1 -1.5139769 4.207542 0.39146424
gene_2 -7.7396297 3.657546 0.01129089
gene_3 -0.9401836 6.573969 0.53604662
gene_4 -2.5983076 5.023355 0.12333235
gene_5 -0.8773113 7.160117 0.55946082
prior.countを0.125に設定
logFC logCPM PValue
gene_1 -1.5139769 4.207542 0.39146424
gene_2 -7.7396297 3.657546 0.01129089
gene_3 -0.9401836 6.573969 0.53604662
gene_4 -2.5983076 5.023355 0.12333235
gene_5 -0.8773113 7.160117 0.55946082
全遺伝子については大変なので比較しませんが、最初の5遺伝子を比べると、logFCを含めて全ての値がデフォルト設定とprior.countを0.125にしたものが同じ値となっており、デフォルトが0.125だったのは本当のようです。
まとめ
長くなりましたが、備忘録として、RパッケージのedgeRでリードカウントが0の場合になぜFold Changeが計算できるのかを調べてみました。
結論としては、edgeRのexactTest関数内でprior.countというArgumantにより、TMM正規化後の数値に特定の数値(デフォルトは0.125)を足すことでリードカウントが0にならないようにし、Fold Changeの算出(つまりリードカウント同士の割り算)ができるようになっていた、ということでした。