RでWGCNA解析～Soft Thresholdの選択をR上で行う～
前回の続きです。スケールフリー性を持つようなネットワーク構築のためにSoft Thresholdの選択が大事であると言いました。ここではサンプルデータを用いて実際にSoft Thresholdの選択とそのR^2値を求めて妥当性の評価を行っていきます。
今回はWGCNAチュートリアルにあったサンプルデータを使います。以下からダウンロードしましょう。
Soft Threshold(sft)の選択
まずはデータの読み込みとライブラリの準備です。データ自体がR1C1参照方式になっていることも確認してください。
library(WGCNA);
options(stringsAsFactors = FALSE);
#データを行列として読み込み
Data1 = read.csv("LiverFemale3600.csv",header=T,row.names=1)
Data2<- t(as.matrix(Data1))
#ここの転置しないと組織ごとにモジュールを検出してしまいます次にこれを使ってsftを選択するためのグラフをプロットします。データによって変えるべき場所は説明するので、それ以外はそのままコピペで大丈夫です。
#R^2値の計算
powers = c(1:30)
sft <- pickSoftThreshold(Data2, powerVector = powers, verbose = 5)
#計算した値のプロット
cex1 = 0.9;
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
abline(h=0.90,col="red")ここでの"power"はsftそのものになります。前回の記事でwgcnaでは類似度を冪乗することでノード同士のつながりをより強調してスケールフリー性を持たせるといいましたが、ここではとりあえず1乗から30乗まで試しています。powers=c(1:30)がそれに該当します。対象データによって値を変えましょう。
あとはそれをグラフとしてプロットします。縦軸をR^2値、横軸をsftとしています。下に結果を示し、そのグラフを読み取っていきます。
基本的にsftを上げていくと、R^2値は頭打ちになります。今回の場合9~10の辺りからほぼ横ばいになっているのが分かります。前回の記事でも言いましたが、スケールフリー性を持たせることのできるsftの中で最小のものを選び取りたいので、今回は9か10を選択することになります。どちらを選択するかは恣意的になるので、そこを理解したうえで解析を進めましょう。今回は10を選択します。
コードの中の"cex"は文字の大きさを指定するもので、今回は0.9にしています。今回は30乗までしか計算していませんが、場合によっては100乗まで見ることもあるので、文字が大きいとグラフが見づらくなります。適宜変更しましょう。これはあってもなくても良いです。
"abline"はsftの選択をしやすくするためにラインを引くためのコードです。今回は高さ0.90に引きましたが、これも適宜高さを変えましょう。目視できる場合は引かなくても良いです。
各ノードがどれくらいの次数を見てみる
powerによって各ノードの次数(ほかのノードとの類似度の総和)の平均値を見てみます。スケールフリーなデータの場合、ある特定のノードだけがおおきなハブになるので、あるsftを境に小さくなるはずです。
#スケールフリー性を持つかどうかを見る
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
おそよ10近辺でほぼ横ばいになっていることが分かりました。ある程度スケールフリー性が確認できたので、次回はこの冪乗した類似度を用いてTOMの計算、それからもう一度スケールフリー性が担保されているかをR上で検証するパートになります。
↓前回
↓次回の記事