・お題：前回の続き。Differential expression analysis（DEG）という解析を少し齧ってみたい。私はWet専門で学んできた人間だから、この辺の解析方法にはとんと疎い。以下のHPを順になぞって学んでいきたいと思う。
・今回は前回の続きであり、ヒートマップなどを使って、どのサンプルとどのサンプルが似ているのかをざっくり見たあと。参考サイトでいうと、Normalisation for composition biasの直前まで終わっている感じ。
・この後、遺伝子発現が細胞の状態や種類によってどのように変化しているのか調べていく。1つのサンプル内での発現の高い低いではなく、サンプルを跨いで比較していくことになるので、サンプルのバイアスを補正してあげてから、比較したい。
・今回は、edgeRのcalcNormFactors関数を使って、TMM正規化してみる。TMM正規化に関して、詳細は以下を参照。「バイアスがかかっていなければ、ハウスキーピング遺伝子はほぼ発現変動しない筈だよね」という考え方をもとにした補正方法らしい。
・まず、元データを確認する。
> y
・次に、calcNormFactors関数を適用するとどうなるか見てみる。
> y <- calcNormFactors(y)
> y
・何が変わったかというと、y$samplesのnorm.factorsが1から1以外の値に変わっている。これが補正値になるらしい。補正前と補正後でMDプロットを並べてみた（左側が補正前、右側が補正後）。
> par(mfrow=c(1,2))
> plotMD(logcounts,column = 7, ylim=c(-8,8), yaxp=c(-8, 8, 8))#補正前
> abline(h=0,col="pink")#ゼロのところにピンクのラインを引く
> plotMD(y,column = 7, ylim=c(-8,8), yaxp=c(-8, 8, 8))#補正後
> abline(h=0,col="pink")
・補正後で、プロットが全体的に下に移動しており、平均がゼロに近くなっている。データそのものを弄ったのではなく、補正係数が弄られる形になるらしい。
・次に、limmaというパッケージのvoomという関数を使って、遺伝子発現解析を進める。まず、design matrixというやつを作成する。
・今回、サンプル構成の情報はgroupに入っている。
> group
[1] basal.virgin basal.virgin basal.pregnant basal.pregnant basal.lactate basal.lactate luminal.virgin luminal.virgin luminal.pregnant
[10] luminal.pregnant luminal.lactate luminal.lactate
Levels: basal.lactate basal.pregnant basal.virgin luminal.lactate luminal.pregnant luminal.virgin
・このgroupをもとにdesign matrixを作成する。design matrixというのは、各サンプルが、どの試験条件（2細胞種×3マウス状態＝6条件）に属するものかを示す行列らしい。
> design <- model.matrix(~ 0 + group)
> design
・design matrixでは、どのサンプルがどの試験条件なのか、が0と1で記されている。今回、replicateが2なので、例えばサンプル1とサンプル2に同じ試験条件が付与されている。designの列名が「group何とか」になっていて鬱陶しい。列名を少しシンプルにする。
> colnames(design) <- levels(group)#列名にもともとのgroup名を設定。
> design
・次に、voom変換する。引数plotにTRUE を与えると、平均分散トレンドを図示してくれる。データに変数があるっぽい遺伝子があるか、低いカウントを適切にフィルタリングしたか情報を得られるらしい。
・voomに関して正しく知りたい方は、以下の論文を参照いただきたい。私は理解が追い付いていない。。
genomebiology.biomedcentral.com
> v <- voom(y,design,plot = TRUE)
・ついでに、vが何か見てみる。
> v
・v$Eに補正後のlog2 countsが、v$weightsに線形モデルに適用される重みが格納されているっぽい。補正後のカウント数の分布が、補正前と比べてどうなっているのか、箱ひげ図を描いて調べる。
> par(mfrow=c(1,2))
> boxplot(logcounts, xlab="", ylab="Log2 counts per million",las=2,main="Unnormalised logCPM")
> abline(h=median(logcounts),col="blue")
> boxplot(v$E, xlab="", ylab="Log2 counts per million",las=2,main="Voom transformed logCPM")
> abline(h=median(v$E),col="blue")
・次に、発現の差異を調べてみる。limmaパッケージのlmFit関数を各遺伝子の線形モデルに当て嵌める。
> fit <- lmFit(v)
・発現の差異を見るのだから、どの試験条件とどの試験条件の差を見たいのか、指定する必要がある。今回は、細胞種をbasalに固定し、マウスの状態がpregnantとlactateの場合の遺伝子発現を比較したい。帰無仮説はそれぞれの試験条件で遺伝子発現に差がない、basal.pregnant - basal.lactate = 0 と定義される。グループ名はdesign matrixのカラム名と一致させる。
> cont.matrix <- makeContrasts( B.PregVsLac = basal.pregnant - basal.lactate, levels = design)
> cont.matrix
・contrasts matrixはdesign matrixのどれとどれを比較するかを示すmatrixらしい。ちなみに、今回は1対1の比較だけれど、複数の比較もできるらしい。
・次に、limmaのcontrasts.fit関数を用いて、contrasts行列をfitオブジェクトに適用して、興味のある比較統計量と推定パラメータを得る。
> fit.cont <- contrasts.fit(fit, cont.matrix)
・次に、分散に対して経験的ベイズ推計を適用し、moderated t-statisticsと関連するp値を推定する。何をやっているのか全く分からないが、とにかくp値が出るらしい。この辺の理屈は数学の基礎から学びなおしていこうと思うが、本稿でそれをやっていてはいつまでたっても終わらなので、いったん最後までいく。
> fit.cont <- eBayes(fit.cont)
・fit.contの中身はどうなっているかというと、以下のような感じ。
> View(fit.cont)
・ limmaのdecideTests関数で、DEGのサマリーが見れる。
> summa.fit <- decideTests(fit.cont)
> summa.fit
・ちなみに、DEGか否かの判断にはデフォルトの設定が存在する。decideTests(object, method="separate", adjust.method="BH", p.value=0.05, lfc=0)とのこと。詳細は以下を参照。
・DEGの数を数える。上がった遺伝子が2705個、変わってないのが10464個、下がったのが2635個になった。
> summary(summa.fit)
・ここまでで、一応DEGが出てきたことになる。先ほどの例では、細胞種をbasalに固定し、マウスの状態がpregnantとlactateの場合の遺伝子発現を比較した際のDEGを見てみた。せっかくなので、細胞の種類をluminalに変えてみてマウスの状態がpregnantとlactateの場合の遺伝子発現を比較した際のDEGを導き出し、どちらの細胞でもDEGになっている遺伝子、どちらかの細胞でしかDEGになっていない遺伝子、を求めてみたい。
・とりあえず、細胞の種類をluminalに変えてみてマウスの状態がpregnantとlactateの場合の遺伝子発現を比較した際のDEGを導き出す。
> cont.matrix2 <- makeContrasts(L.PregVsLac = luminal.pregnant - luminal.lactate,levels = design)
> fit.cont2 <- fit %>% contrasts.fit(cont.matrix2) %>% eBayes()
> summa.fit2 <- decideTests(fit.cont2)
> summary(summa.fit2)
・どれだけのDEGがそれぞれの比較で固有または共通なのか、ベン図で示してみる。limmaパッケージのvennDiagramが便利。以下を参考にさせて頂いた。
> vennDiagram(cbind( summa.fit, summa.fit2),
+ include=c("up", "down"),#upで上がったもの、downで下がったもの、bothでupまたはdownを示す。
+ counts.col=c("red", "blue"),#先の上がったものの色、下がったものの色を指定。
+ circle.col = c("red", "blue"))#円の色を指定。
・形式上、お絵かきすることはできた。原理がよく分かっていないのは致命的なので、徐々に勉強していきたい。
つづく。
