RでWGCNA解析～クラスタリングとモジュール検出～

前回はWGCNAのメインである重み付きの相関を計算しました。今回はこの結果に基づいてクラスタリングとモジュール検出を行います。モジュール検出を行うことによって関連のある遺伝子を可視化することができます。

クラスタリング

例によって前回の記事から地続きでコードを書くので、セーブした人は読み込んでから始めてください。

#樹形図の作成。距離としてTOMを採用しています。
geneTree = hclust(as.dist(dissTOM), method = "average");

#描写
sizeGrWindow(12,9)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04);

以下が出力です。3600ある遺伝子が分類されています。

モジュール検出

ここからモジュール検出に移ります。ある程度の閾値を設けて、遺伝子を発現量に基づいてグループ分けする作業です。コードは以下の通り。

# モジュールサイズの最小値を決定します。一つのグループの最小の大きさで、小さいほどモジュールが細かく分かれる。
minModuleSize = 30;
# モジュールの特定
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
deepSplit = 2, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize);
table(dynamicMods)

　deepspitは0～4の整数値を入力することで、どの程度大雑把にモジュール分割を行うかを選択できます。0が最も大ざっぱにモジュール分割を行い、４に近づくほどに細かくモジュールを分割します。適当に変えて結果の変化を確かめてみると良いと思います。
　pam stageはどのモジュールにも割り当てられなかった遺伝子を最も近いモジュールに再分配するかどうかを決めるパラメータでTRUEまたはFALSEで指定します。特に、割り当てられなかったものの数が多いときに使用します。
　最後のtable(dynamicMods)で、各モジュールに割り当てられた遺伝子数が表示されます。16のモジュールに分かれているのが分かります。

次に検出したモジュールを図に示します。

#モジュール番号を色に変換
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
#樹形図の下にモジュールの色を追加して表示
sizeGrWindow(8,6)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors")

出力結果です。

モジュール検出ができました。これで基本的なWGCNAの手順は終了です。次回はどの遺伝子がどのモジュールに属しているかをCSVファイルで出力し、ヒートマップに落とし込むとこまでやります。

前回の記事↓

次回の記事↓