いろいろ倉庫

KNIME、EXCEL、R、Pythonなどの備忘録

【R】遺伝子発現解析を学んでみたい⑤

・お題:前回の続き。Differential expression analysis(DEG)という解析を少し齧ってみたい。私はWet専門で学んできた人間だから、この辺の解析方法にはとんと疎い。以下のHPを順になぞって学んでいきたいと思う。

combine-australia.github.io

・今回は前回の続きであり、DEGを導き出したりベン図を描いたりしてみたところ。参考サイトでいうと、Plots after testing for DEの直前まで終わっている感じ。

・解析がおかしくないことを確認するために、いろいろ図示してみる。ゲノム全体のプロットをMAplots (または MDplots)とvolcano plotsで確認してみる。limmaのfit.contをグラフを描くための関数に放り込む。

・とりあえず、fit.contが何だったか思い出す。

> View(fit.cont)

・次に、グラフを描く。MDplotsはMean-Difference Plotの略で、対数強度比(defference)と対数強度平均(mean)を各遺伝子ごとに算出して、プロットしたもの。横軸に発現量、縦軸に変動の幅をとったようなイメージ。volcano plotsは横軸に対数強度比を、縦軸にp値の対数変換値にマイナスをかけた値をとったグラフ。

> par(mfrow=c(1,2))#グラフを1行2列に配置する。
> plotMD(fit.cont, coef=1, status=summa.fit[,"B.PregVsLac"], values = c(-1, 1),  hl.col=c("blue","red"))#MDプロットを描く。
> volcanoplot(fit.cont, coef=1, highlight=100, names=fit.cont$genes$SYMBOL,  main="B.PregVsLac")#volcanoプロットを描く。

・次に、stripchart関数を使って、個々のサンプルの発現レベルを見てみる。voomオブジェクト(v$E)の正規化されたlog expression値を引数に使う。

> v$E %>% head()#v$Eが何だったか確認してみる。

・indexが24117になっているやつを取り出して、グループごとに発現量を図示する。

> nice.col <- brewer.pal(6,name="Set2")
> stripchart( v$E["24117",]~group, vertical=TRUE, las=2, cex.axis=0.8, pch=16,cex=1.3, col=nice.col, method="jitter", ylab="Normalised log2 expression", main="Wif1")

・発現量は群内で綺麗に揃っている。GlimmaパッケージのglXYPlot関数を使うと、インタラクティブ版が作れる。
> group2 <- group
> levels(group2) <- c("basal.lactate", "basal.preg", "basal.virgin", "lum.lactate", "lum.preg", "lum.virgin")
> glXYPlot( x=fit.cont$coefficients[,1], y=fit.cont$lods[,1], #volcano plotのx軸データとy軸データをfit.contから取得。
+          xlab="logFC", ylab="B", main="B.PregVsLac",#volcano plotの軸とタイトルを指定。
+          counts=v$E, groups=group2, status=summa.fit[,1], #発現量、群、色基準を指定。
+          anno=fit.cont$genes, side.main="ENTREZID", folder="volcano")#フロート情報、発現量プロットのタイトル基準、htmlを保存するフォルダ名を指定。

・ボルケーノプロットにカーソルを合わせると、情報がフロートするうえ、自動でstripプロットが作成される。

・気になる遺伝子をリストから検索することもでき、気になる遺伝子がボルケーノプロットのどこにあって、各群の遺伝子発現がどうだったかが出てくる。

・ここまでで、遺伝子発現データセットからたくさんのDEGを見つけてきて、図示したり発現を確認したりすることができた。ここからさらに、p値の閾値だけでなく、発現の変動倍率の閾値でカットオフしたりする。単純にp値とlogFCだけを基準にDEGを見直すとエラーが増えるので、limmaのtreat関数で解析する。treat関数は、fit.contオブジェクトと、ユーザー指定のlog fold changeカットオフを渡すことで、moderated t統計とp値を再計算してくれる。
・とりあえず、やってみる。

> fit.treat <- treat(fit.cont,lfc=1)#fit.contオブジェクトと、log FCの基準を渡す。lfc=1は、log2換算で絶対値1が閾値になるので、0.5倍と2倍が閾値になる。
> res.treat <- decideTests(fit.treat)#上がった下がった変わらないを1, -1, 0に振りなおす。
> summary(res.treat)#上がった、下がった、変わらない、を集計。

・MAプロットを描いてみる。

> plotMD(fit.treat,coef = 1, status = res.treat[,"B.PregVsLac"], values = c(-1,1),  hl.col = c("blue","red"))
> abline( h = 0, col = "grey")

・ちなみに、glMDPlot関数を使えば、MDプロットもインタラクティブ化できる。

> glMDPlot(fit.treat, coef=1, counts = v$E, groups = group2,
+          status=res.treat,  side.main = "ENTREZID", main = "B.PregVsLac", 
+          folder="md")

 

・これでどうも発現に差があるらしい遺伝子をとってこれることになった。ただ、この遺伝子のリストが膨大な量になった場合、生物学的意義を理解することは困難になる。そこで、これら発現変動が見られた遺伝子が、どのパスウェイ/ネットワークに関与しているかを解析していくことになる。

・生物学的パスウェイの濃縮の検定方法は大きく分けて、self-containedとcompetitive の2種類がある。self-containedは「パスウェイ内の遺伝子は全体として発現変動しているか?」を評価する手法で、ROASTなどがある。competitive は「他のすべての遺伝子と比較して、発現変動した遺伝子がセット内にたくさん出現する傾向があるか」を評価する手法で、goanaやcameraなどがあるらしい。

・まず、limmaのgoana関数を使ってエンリッチメント解析をやってみる。例えば、basal細胞に関して、pregnantとlactateの群を比較し、後者でたくさん出てくるGOタームを調べる。goana関数にfit.contオブジェクトと気になる係数と生物種の略称入力する。最も濃縮されたGOタームのトップセットは、topGO関数で確認できる。

> go <- goana(fit.cont, coef="B.PregVsLac",species = "Mm")#Mmはマウスのこと
> topGO(go, n=10)

・N列の数字は、各GOタームでアノテーションされた遺伝子の総数。Up列とDown列は、GOタームの遺伝子と重複しているDEGの数。P.Up列とP.Down列は増えた又は減った遺伝子のセットで、GOタームの過剰出現のp値を示す。

・更に、遺伝子の長さをcovariateに与えることにより、補正できるらしい(割愛)。この後、バーコードプロットを描いたりする流れになるが、詰まった。他のサイトなども参考にしながら、そのうち続きもやってみたい。

 

おわり。