RでWGCNA解析～Soft Thresholdの選択をR上で行う～

　前回の続きです。スケールフリー性を持つようなネットワーク構築のためにSoft Thresholdの選択が大事であると言いました。ここではサンプルデータを用いて実際にSoft Thresholdの選択とそのR^2値を求めて妥当性の評価を行っていきます。

　今回はWGCNAチュートリアルにあったサンプルデータを使います。以下からダウンロードしましょう。

LiverFemale3600.csv

1.26 MB

ファイルダウンロードについて

ダウンロード

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上で検証するパートになります。

↓前回

↓次回の記事