いろいろ倉庫

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

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

・お題:他にも遺伝子発現解析の例があったので、そちらもなぞってみたい。

 

・次の例は以下。とても詳しく記載してあるので、きちんとしたことは元サイトをご確認いただきたい。

bioconductor.org

・今回は前回の続きで、遺伝子発現の分布をTMM normalizeしたあと。

・探索的に遺伝子発現の関係を眺めるために、MDS(Multi-Dimensional Scaling)プロットを作成する。

・GlimmaパッケージのglMDSPlot関数を使う。

> glMDSPlot(lcpm, labels=paste(group, lane, sep="_"), 
+           groups=x$samples[, c(2, 5)], launch = TRUE)

・第1主成分と第2主成分で展開してみると、それぞれの細胞集団でくっきり分かれた。第1主成分ではbasalとその他が明らかに分かれている。それぞれの細胞集団を一対一で比較した際に、basalサンプルとその他の細胞集団の比較は、LPとMLを比較した場合よりも発現に差のある遺伝子がたくさんありそうと推察できる。

・このあと、発現に差がある遺伝子は何なのか、解析していくことになるが、その前に、design matrixを作成する。今回は、3種類の細胞種(basal、LP、ML)とバッチ情報(L004、L006、L008)が9つのサンプルにどう割り振られているかを示すdesign matrixを作成する。

・design matrixに関しては、以下を参照。なるほど分からない。。

www.ncbi.nlm.nih.gov

・今回は細胞ごとに比較したいので、basal、LP、MLをグループとして振り、切片は無くて良いので~0になりそう。バッチ情報を因子として加味すると、Studies with multiple factorsの記載が最も近いような感じになりそう。

> design <- model.matrix(~0 + group + lane)

> design

・列名がgroupナントカになっており、少々気持ち悪い。groupは余分なので、gsub関数で置換して消す。

> colnames(design) <- gsub("group", "", colnames(design))
> design

・この後、特定の細胞の遺伝子発現を比較することになる。contr.matrixを作成する。

> contr.matrix <- makeContrasts(
+   BasalvsLP = Basal-LP, 
+   BasalvsML = Basal - ML, 
+   LPvsML = LP - ML, 
+   levels = colnames(design))
> contr.matrix

RNA-seqカウントデータでは、分散と平均は独立ではないので、voomで正規化係数を抽出し、補正する必要がある。この関数では、voom-plotと呼ばれるlogカウントと分散の関係を示すプロットも作図してくれる。低発現の遺伝子をうまいこと除去できていないと、低発現のところで分散が異様に小さくなるらしい。

> v <- voom(x, design, plot=TRUE)

・ぱっとみ緩やかな右肩下がりのグラフであり、低発現で分散の低下はみられていない。低発現の遺伝子を上手いこと除けているらしい。

・次に、線形モデルにフィットさせる。

> vfit <- lmFit(v, design)
> vfit <- contrasts.fit(vfit, contrasts = contr.matrix)
> efit <- eBayes(vfit)

・さらに、発現データの平均-分散関係を図示して分布を確認する。plotSAで残差-分散対平均対数発現の散布図が得られることになる。

> plotSA(efit, main = "Final model: Mean-variance trend")

・発現量に応じて残差の傾向が変動している、ということはなさそうなので、確かに補正されているっぽい。
・次に、有意に上がった/下がった遺伝子の数を調べてみたい。デフォルトで5%に設定されている調整済みp-valueカットオフを使用して有意かそうでないか判別する。

> efit$p.value %>% head()#efitからp.valueを取り出す。

> decideTests(efit)#上がった、下がった、変わらないを1,-1,0で示す。

> summary(decideTests(efit))#上がった、下がった、変わらない遺伝子の数を集計。

・basal細胞はLP細胞と比べて4,648遺伝子が発現低下、basal細胞はLP細胞に対して4,863遺伝子が発現上昇したことになっている。先のMDSプロットで見られた「basal細胞は他と比べて遺伝子発現プロファイルが結構違うっぽい」という直感に一致した発現変動遺伝子の数になっている。

・p値だけではなく、変動した倍率に閾値を設けて、発現変動遺伝子を抽出することもできる。

> tfit <- treat(vfit, lfc=1)#lfc=1でlog2 fold changeが1以上のものをとって来る。
> dt<- tfit %>% decideTests()

> dt %>% summary()

・lfcを設定しない場合と比べ、確かに数が減っている。

・せっかくなので、これらの関係を図示してみる。BasalvsLPとBasalvsMLのUpとDownをベン図上に示してみる。

> vennDiagram(dt[,c(1,2)],
+             include=c("up", "down"),
+             counts.col=c("red", "blue"),
+             circle.col = c("gray", "green"))

・また、結果を出力することもできる。
> write.fit(tfit, dt, file="results.txt")

・top DE遺伝子はtreatの場合はtopTreat関数、eBaysの場合はtopTable関数でとってこれる。

> basal.vs.lp <- topTreat(tfit, coef=1, n=Inf)#coefは列番号の指定。1列目はbasalとLPの比較。
> head(basal.vs.lp)

・遺伝子発現解析の結果を視覚化するのに、横軸に平均log-CPMを、縦軸にlog FoldChangeをとったMDプロットを表示する。GlimmaのglMDPlot関数を使うと、インタラクティブなプロットが出てくる。

> glMDPlot(tfit, coef=1, status=dt, main=colnames(tfit)[1],
+          side.main="ENTREZID", counts=lcpm, groups=group, launch=TRUE)

・左のMDプロットの各点は遺伝子であり、気になる点をクリックすると、右に個々のサンプルの遺伝子発現量(log-CPM)、プロットの下に種々の情報が出てくる。

・gplotsパッケージのheatmap.2関数を使って、basalとLPの比較からadjusted p-value上位100のDEGについてヒートマップを作成する。

> library(gplots)

> basal.vs.lp.topgenes <- basal.vs.lp$ENTREZID[1:100] #上位100遺伝子を抜き出す
> i <- which(v$genes$ENTREZID %in% basal.vs.lp.topgenes) #v$genes$ENTREZIDのうち、先の上位100遺伝子にあるものを選抜
> mycol <- colorpanel(1000, "blue", "white", "red") #ヒートマップ用のグラデーションのかかったカラーコードを作成。

・次に、heatmap.2関数でヒートマップを描く。引数の詳細は以下を参照。

www.rdocumentation.org

> heatmap.2(lcpm[i,], scale="row",
+           labRow=v$genes$SYMBOL[i], labCol=group, ColSideColors=col.cell,
+           col=mycol, trace="none", density.info="none", 
+           margin=c(8,6), lhei=c(2,10), dendrogram="column")

・細胞種ごとにサンプルが正しくクラスタリングされており、MLとLPの遺伝子発現は、かなり似ているらしい(今回ヒートマップを描いた範囲は、だけれど)。

つづく。