ex402: Differential expression analysis with edgeR

arab2データの遺伝子発現の2群間比較を、edgeRで行う。

arab2.txt は、ex301, ex302でもあつかった、シロイヌナズナRNA-seqのデータである。細菌の感染と非感染でトランスクリプトームを比較している。実験デザインに関する記述を原著論文から引用する。

Bacteria were grown in King’s B media and infiltrated into plants as previously described. Briefly, we used a syringe lacking a needle to infiltrate the abaxial side of leaves of six-week old Arabidopsis plants. Plants were infected with either the ∆hrcC mutant of Pseudomonas syringae pv. tomato DC3000 (PtoDC3000) or mock inoculated with 10 mM MgCl2 7 hpi. Each treatment was done as biological triplicates with each pair of replicates done at separate times and derived from independently grown plants and bacteria. -- Cumbie et al., 2011 PLOS ONE

つまり、6 libraries (2 groups x 3 biological replicates) のデータ。すでにマッピング済みで遺伝子毎のリードカウントがタブ区切りテキストとして提供されている。m1, m2, m3 が mock (control), h1, h2, h3 が　Pseudomonas感染。

Import data

> x <- read.delim("arab2.txt", row.names=1)

ここまでは、ex301, ex302と同じ。

Import edgeR library

> library(edgeR)

2グループ、各3繰り返し実験、という実験デザインを定義する。（m1, m2, m3 が mock (control), h1, h2, h3 が　Pseudomonas感染。）

> grp <- c("M", "M", "M", "H", "H", "H")
> grp
[1] "M" "M" "M" "H" "H" "H"

edgeRのDGEList関数でカウントデータを読み込む。

> y <- DGEList(counts=x, group=grp)

DGEListオブジェクトはcountsとsamplesという重要なコンポーネントを持つ。countsはDGEListを作成するときにcountsで指定したデータ（matrix）。samplesはサンプル情報。ここで確認してみよう。

> y$samples
   group lib.size norm.factors
m1     M  1902032            1
m2     M  1934029            1
m3     M  3259705            1
h1     H  2129854            1
h2     H  1295304            1
h3     H  3526579            1

Filtering genes by expression level (optional)

発現量が全体的にゼロに近い遺伝子はDEG解析において情報量がないので除去するのが望ましい。ただし、このフィルタリング処理はやや発展的な内容なので、基礎を理解する過程ではスキップしてかまわない。（フィルタリング処理をスキップしても結果は大きくは変わらない。余裕がある方はフィルタリングの有無で結果を比較しても良いでしょう）

ここではCPM(count-per-million)が0.5以上のサンプルが３つ以上、をフィルタリングの条件とする。

> keep <- rowSums(cpm(y) > 0.5) >= 3
> table(keep)
keep
FALSE  TRUE 
 6663 19558 

y <- y[keep, , keep.lib.sizes=FALSE]

参考）edgeRには発現量が低い遺伝子を判定するコマンド、filterByExprがあり、これを使っても良い。使い方はマニュアル参照。ただデフォルトのパラメータのままだと厳しい印象。

参考）data.frame(y, keep)　を実行すると、どの遺伝子がkeepされたかわかる。

Normalization

TMM法で、ノーマライズする。normLibSizesを使う。(注：edgeRのバージョンが３系の場合は、calcNormFactors)

> y <- normLibSizes(y, method="TMM")

計算結果の確認

> y$samples
   group lib.size norm.factors
m1     M  1901107    1.0247736
m2     M  1917419    1.0537590
m3     M  3257400    0.8958933
h1     H  2126850    1.0297446
h2     H  1280838    1.1293858
h3     H  3524342    0.8887973

(Option) Normalize効果の可視化

> par(mfrow=c(1,2))
> boxplot(log10(x + 1))    # before normalization
> boxplot(log10(cpm(y) + 1)) # after normalization

注）par(mfrow=c(1,2)) の効果を無効化するには、 dev.off() 。もしくは、par(mfrow=c(1,1))でもよい。

DE testing

Estimate dispersion

負の二項分布の確立分布モデルに当てはめる。

> y <- estimateDisp(y)

(参考)単純な２群間比較の場合(正確にいうと１要因；edgeRではclassic modeと呼ばれる)、edgeRではdispersionの推定に、quantile-adjusted conditional maximum likelihood (qCML) method が内部的に使われる。

（参考）内部的には、全体のdispersionを推定するのと、遺伝子ごとにdispersionを推定する２ステップで処理される。この２つのステップを明示的に分けて実施することもできる；estimateCommonDisp()を先に計算してそのあとestimateTagwiseDisp()。

推定されたdispersionを確認してみよう

> y$common.dispersion
[1] 0.2647946

> summary(y$tagwise.dispersion)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.07852 0.13140 0.22346 0.43488 0.55465 2.69212 

(注：edgeRのバージョンで値は少し異なる。上記は、edgeR ver. 4.0.16 での結果)

DE testing and multiple comparison correction

上で推定した確立分布に基づいて正確確率検定を行う。

> et <- exactTest(y)

上位のDEG遺伝子(deferentially expressed genes)を確認する。topTagsコマンドを使う。このコマンドはmultiple comparisonの補正も同時に計算してくれる。補正の方法はadjust.methodオプションで指定する。デフォルトのBH (=FDR)で通常は問題ない。

> topTags(et, n=10, adjust.method="BH")

Comparison of groups:  M-H 
              logFC    logCPM       PValue          FDR
AT2G19190 -4.539137  7.388103 1.001454e-24 1.958644e-20
AT2G44370 -5.428882  5.211831 2.061797e-22 2.016231e-18
AT3G55150 -5.774035  4.916743 9.049096e-22 5.899407e-18
AT1G51820 -4.379040  6.385171 3.111640e-21 1.521436e-17
AT4G12500 -4.389112 10.442749 4.121689e-21 1.612240e-17
AT3G46280 -4.911817  8.119364 7.044169e-21 2.251923e-17
AT2G39380 -4.941156  5.785228 8.059852e-21 2.251923e-17
AT5G64120 -3.747018  9.703124 9.831508e-21 2.403558e-17
AT2G39530 -4.344088  6.721770 3.069247e-20 6.669814e-17
AT4G12490 -3.869536 10.203920 5.468145e-20 1.069460e-16

(注：edgeRのバージョンで微妙に結果は変わる。上記は、edgeR ver. 4.0.16 での結果)

Dump the result into a text file

exactTestの全遺伝子の結果をFDRの値付きで、タブ区切りテキスト形式で出力する。

> et.sorted <- topTags(et, n=nrow(et$table), adjust.method="BH")
> write.table(et.sorted$table, "et.sorted.txt", sep="\t", quote=F)

適切にファイルにデータが出力されたか確認する。MacやWindowsでは作業していたフォルダを開いて、et.sorted.txtというファイルができているかをチェックする。エディタや表計算ソフトウェアでオープンして中身を確認することもできる。UNIXコマンドラインからは、less などのコマンドで中身を確認することができる。

Inspect DE result

DEG遺伝子の数を確認しよう

example: fold-change > 10を抽出・カウント

## 準備
> et.sorted <- topTags(et, n=nrow(et$table))
> detab <- et.sorted$table
## get fold-change > 10
> detab[detab$logFC > log2(10),]
# => fold-change > 10のものが出力される。数だけ知りたい時は、
> nrow(detab[detab$logFC > log2(10),])
[1] 30

example: FC > 5 AND FDR < 0.01

> detab[(detab$logFC > log2(5)  & detab$FDR < 0.01), ]

edgeRに組み込まれている''decideTest''関数も便利。(edgeR 3系の場合は、decideTestsDGE)。上と同じ条件を満たす遺伝子を抽出するには、

> summary(decideTests(et, p.value=0.01, lfc=log2(5)))
         M-H
Down     130
NotSig 19417
Up        11

dumpしたタブ区切りテキストをMS Excelで読み込んで、フィルタ機能やソート機能を使ってデータを探索するのも良いだろう。

MA plot　を描いてみよう。

ex302 で、MA plotを描く学習をしたが、ここでは、edgeRに含まれているのplotSmear関数でMA plotを描いてみる。

plotSmear(y)

有意に発現差のある遺伝子を赤色でハイライトすることもできる。

>  de.names <- row.names(et[decideTestsDGE(et, p.value=0.05, lfc=log2(2)) !=0, ])
> plotSmear(y, de.tags=de.names)

Links

• http://www.bioconductor.org/packages/release/bioc/html/edgeR.html | edgeR
• http://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf | edgeR User's Guide

Notes

最近edgeRのバージョンが４に上がって、コマンドが変更になったものがいくつかある。本文には３系と４系の違いがわかるように記述した。

Revision history

2024-7-30

• edgeR v4系とv3系の両方に対応するように記述をアップデート。
• filter by expression level をオプション扱いに変更。

2024-2-25

• Major update.
• R v4.3.2; edgeR v4.0.16 on local Mac で動作確認
• 最新のedgeRのマニュアルに沿って改訂。
• filter by expression level を追加。

2023-2-24

• Minor updates.
• R ver.: 4.0.3; edgeR ver.: 3.32.1 on bias5 + X11 で動作確認

2022-8-27

• R ver.: 4.0.3; edgeR ver.: 3.32.1 on bias5 で動作確認

2021-0-11

• Minor updates.

2021-3-9

• Minor updates. R ver.: 4.0.4; edgeR ver.: 3.32.1