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の算出(つまりリードカウント同士の割り算)ができるようになっていた、ということでした。