Part1では正規化などのデータの前処理、Part2ではネットワークの構築とデンドログラムの作成を行ってきました！
Weighted Gene Coexpression Network Analysis (WGCNA)は遺伝子発現量の相関を利用して、互いによく似た発現変動をする遺伝子同士のクラスターを発見する手法です！ WGCNAは単純にクラスター[…]
前回のPart1ではrawデータのダウンロードから正規化を行い、サンプルと遺伝子フィルタリングまでを行いました。 [sitecard subtitle=関連記事 url= https://shiokoji11235.com[…]
今回はWGCNAで最も重要なパートであるモジュール検出を行います！
モジュール検出は自分で設定しなければならないパラメータが多くあり、解析者によって得られる結果が大きく異なります。ですので一つ一つのパラメータの意味を理解してデータに適したものを選択していくのがとても重要になります！
各パラメータについてひとつひとつ丁寧に説明していきますのでお付き合いください！
deepSplit: cutreeDynamicの引数。どれだけ細かくモジュールを分割するかを決める。
pamStage: cutreeDynamicの引数。モジュールに割り当てられなかった遺伝子を近傍のモジュールに再割当てを行うかどうかを決める。
module eigengene: モジュールに所属する遺伝子の発現マトリクスの第一主成分
準備
まずは今回使用するライブラリーをインストールします！
library(limma)
library(gplots)
library(tidyverse, warn.conflicts=F, quietly=T)
library(WGCNA)次にPart1で作成した正規化済み発現データとPart2で作成したネットワーク、デンドログラムをロードします！
lnames <- load('../Output/001_NormalizeRawData/SLE_GSE138458-01-dataInput.RData')
print(lnames)
# [1] "datExpr" "datTraits"
lnames <- load('../Output/002_NetworkConstruction/SLE_GSE138458-02-networkConstruction.RData')
print(lnames)
# [1] "adjacency" "dissTOM" "geneTree" ちゃんとロードできていますね！
前回と同じようにデンドログラムも描画してみます！
plot(geneTree, xlab="", sub="", main = "TOM-based dissimilarity", labels = FALSE, hang = 0.04)
モジュール検出の概要
ここからは前回作成したgeneTreeにもとづいて遺伝子のクラスターであるモジュールを切り出していきます！
近い枝同士を同じモジュールに切分けていくイメージを持っていただければこれからやることが理解しやすいと思います！
モジュール検出にはWGCNAパッケージのcutreeDynamic関数を使います！
モジュール検出用の関数としては他にもcutreeStaticなどの関数も用意してくれていますが、
基本的にはcutreeDynamic関数を使うことが推奨されています。
実際、WGCNAチュートリアルでもcutreeDynamicが使われています。なぜcutreeDynamicを使うべきなのでしょうか？
それは複雑な構造のgeneTreeに対しても自然なモジュールの割り当てを行うためです！
推奨されていない方の手法であるcutreeStatic（Dymanicじゃないよ！）はcutHeightを1つ指定して、その高さ以下でつながっている遺伝子同士を同じモジュールに割り当てます。
とても簡単なデンドログラムでの例を以下に示しています。
cutreeStaticは木の高さを指定して、それ以下の高さでつながっている遺伝子同士を同じモジュールとしてクラスタリングを行います。左側のCut Heght を６でクラスタリングしたものでは遺伝子１〜３がつながっていて同じクラスターで４が単独でクラスターとなっています。それに対して右側の左側のCut Heght を３でクラスタリングしたものでは遺伝子１と２が同じクラスターで３と４がそれぞれ単独でクラスターとなっています。
このような単純なデンドログラムではcutreeStaticはとてもうまく働きます！でも今回のようにかなり複雑なデンドログラムではどうでしょうか？
赤線Cut Height = 0.93を示しています。仮にこのCut Heightでクラスタリングを行うなら、右側の青い部分は比較的自然なクラスタリングができそうですが、左側のオレンジ色の部分はもう少しCut Heightを調節したほうがより自然なクラスタリングになりそうです。左側がCut Height = 0.95で右側はCut Height = 0.93でクラスタリングしたいな〜ていう願いを叶えてくれるのが、cutreeDynamic関数です！
cutreeDynamicが何をしてくれるかをざっくりと説明します！まずは枝の先（コア）から順に木をさかのぼっていきます。そして２つの枝がぶつかったときにshape criteriaを判定し、shape criateriaをみたしていればその２つの枝をマージします。それを木の最上部まで繰り返してマージされた遺伝子郡をモジュールとします。このshape criteriaがどんなものかというと、「最小遺伝子数」や「枝の先（コア）から２つの枝がぶつかった高さ（Gap）」などです。これについてはWGCNAワークショップのスライドが詳しいのでこちらをご参照ください。
パラメータの設定
さて！cutreeDynamicを使うことがわかったところで実際にこれを使ってモジュール検出をしてみましょう！
cutreeDynamicはクラスタリングする対象となるTreeとTreeを作成するのに使ったdissTOMをインプットとして渡すと各遺伝子のモジュールへの割り当てを行ってくれる関数です！
この際にいくつかの重要なパラメータがあります！
- deepSplit
- pamStage
順に説明していきます！
deepSplit
deep splitはこの３つのなかでも最も重要であると考えていいと思います！
deep splitは「どれくらい細かくモジュールを分けるか」を決めるパラメータで、0~4の整数をとります！
0が最も大ざっぱにモジュール分割を行い、４に近づくほどに細かくモジュールを分割します。つまりアウトプットされるモジュールの数は０が一番少なく、４が一番多くなります。
pamStage
pam stageは「どのモジュールにも割り当てられなかった遺伝子を最も近いモジュールに再分配する」かどうかを決めるパラメータでTRUEまたはFALSEをとります！
Partitioning Around Medoids (PAM)という教師なし機械学習の一種のアルゴリズムに基づいた処理なのでこの名前になっています！
pamStageについてはまずFALSEで試してみて、どのモジュールにも割り当てられない”grey”モジュールの遺伝子が多かったらpamStageをTRUEに変更するのをオススメします！なぜならpamはモジュールに割り当てるのが難しい遺伝子を強制的に近傍のモジュールに割り当ててしまうため生物学的な解釈が難しくなってしまう傾向があるからです。
pamStageに付随するパラメータとしてpamRespectsDendroというものもあります！こちらは「pam stageにおいてモジュールへの再割り当てでデンドログラムを参考にするかどうか」を決めます！これもTRUEまたはFALSEをとります！
pamRespectsDendroがTRUEなら遺伝子の再割り当ては、「その遺伝子とデンドログラム上で同じ枝に存在しているモジュールだけを対象」にして行われます！
pamRespectsDendroもモジュール割り当てを大きく変える重要なパラメータで、WGCNA Tutorialでも使用（pamRespectsDendro=FALSEが選択）されています！pamStageをTRUEとすると決めたならpamRespectsDendroについてもTRUEとFALSEどちらも試してみるのがいいと思います！
その他のパラメータ
モジュール割り当てに非常に大きく関わり、最優先にに検討すべきパラメータは上の２つですが、他にも多くのパラメータがあります。
上のcutreeDynamicは何をしてるの？の項で少しだけ触れたshape criteriaに関わるものとしてmaxCoreScatter、minGap、minClusterSizeなどがあります。これらについてはデフォルトの値で十分に良い結果を出してくれます！ですのであまり検討する必要はないと思いますが、もしdeepsplitとpamStageを変えても最適なモジュールが得られなかった際にはここを変えてみてもいいかもしれません。
モジュール検出の実行
非常に重要なパートだったため説明がとても長くなってしまいましたが、ここからコードを実行していきます！検討の結果、今回は以下のパラメータを採用しました！
pamStage = TRUE
pamRespectsDendro = TRUE
minModuleSize = 30
deepSplit <- 4
minModuleSize <- 30
dynamicMods <- cutreeDynamic(dendro = geneTree,
distM = dissTOM,
deepSplit = deepSplit,
pamStage = TRUE,
pamRespectsDendro = TRUE,
minClusterSize = minModuleSize)
print(table(dynamicMods))cutreeDynamicはモジュールを数字のラベルで返してくれます！どのモジュールにも割り当てられなかったもの（”grey”モジュール）は0のラベルがつけられます。また、最も遺伝子数の多かったモジュールが1でそこから遺伝子数の順番に2,3,4…と割り当てられていきます！
数字のままではモジュールを区別するのに苦労するので数字ラベルを色に変換します！
dynamicColors <- labels2colors(dynamicMods)
table(dynamicColors)色に変換できたのでデンドログラムと一緒にモジュールの割り当てをカラーバーで見ていきましょう！
plotDendroAndColors(dendro = geneTree,
colors = dynamicColors,
groupLabels = "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors")
似たモジュールのマージ
ここまででDynamicTreecutでモジュール検出ができました！さて、ここからはモジュールの微調整を行っていきます。
具体的に何をしていくかと言うと「module eigengeneが強く相関しているモジュール同士をくっつけて一つのモジュールに統合」していきます！モジュールのマージは必須ではないので、マージを一度試してみて、マージしていない結果のほうが生物学的に理解しやすい結果が得られているならマージをしない結果を採用するのがいいと思います！
いくつか理由はありますが、一言でいうと「module eigengeneが強く相関しているなら、それらのモジュールの中の遺伝子同士の発現プロファイルも相関している可能性が高いから」です！
まずmodule eigengeneとはそのモジュールの遺伝子発現マトリクスに主成分分析（PCA）を行ったときに得られる第一主成分です！主成分分析の結果なのでこれは各遺伝子の発現量の線形結合で以下のように表せます！qモジュールのeigengeneを$E^{(q)}$として各遺伝子の発現量を$Gene_i$、その遺伝子発現量に対する係数を$w_{i}$とすると
$$E^{(q)}=w_{1}Gene_1 + w_{2}Gene_2 + w_{3}Gene_3 … $$
のように表せます！もっとシンプルに説明するとeigengeneは「遺伝子発現量の重み付き平均」と考えられます！例えばこのeigengeneがRedモジュールとBlueモジュールで相関するならば、Redモジュールの中の遺伝子の発現量とBlueモジュールの中の遺伝子の発現量も相関しそうですよね？WGCNAの目的は発現プロファイルが似ている遺伝子同士をクラスター化して生物学的な解釈をしやすくすることです！eigengeneが相関する（≒遺伝子発現プロファイルが相関している）モジュール同士をくっつけて一つにするのはこの目的に沿った手法だと考えられると思います！
もう少しだけ補足するとmodule eigengene間で相関ネットワークを考えたときに似たモジュール同士でクラスターができるんですがこれをメタモジュール（meta-module）と呼びます！このメタモジュールを一つに統合する過程がこのモジュールのマージというステップです！
メタモジュールは病態スコアなどのTrait情報とモジュールとの関係性を見たいときにすごく有用です！module eigengeneにTrait情報も含めて相関ネットワークを描き、Trait情報と強く相関するmodule eigengeneを持つモジュール達をメタモジュールとして統合しパスウェイ解析などを行うことで、そのTraitと関わる機能を調べたりもできます！
実際にモジュール同士をマージするかどうかは「Gene ontologyやパスウェイ解析の結果を見て」決めます。
例えばRedモジュールとBlueモジュールがどちらも似た機能を持っているならマージした方が解釈しやすいと思います！逆に、たとえmodule eigengenが強く相関したとしても全く異なる機能を持っているモジュール同士であればマージしないほうがいいのではないかと考えられます。
このあたりに絶対的な正解はなく解析者の考え方や持っている仮設に依存すると思うのであまり深く考えすぎる必要ないと思います！
モジュールのマージの具体的な手順は以下です！
- moduleEigengenes関数で各モジュールのeigengeneを計算
- module eigengen同士の相関を計算
- 相関をもとにモジュール同士の距離を計算
- 階層的クラスタリングを行う
- モジュールを統合する基準を決めてマージする
以下のコードで上記を実行していきます！
MEList <- moduleEigengenes(datExpr, colors = dynamicColors)
MEs <- MEList$eigengenes
MEDiss <- 1-cor(MEs)
METree <- hclust(as.dist(MEDiss), method = "average")
plot(METree, main = "Clustering of module eigengenes",
xlab = "", sub = "")
MEDissThres <- 0.15
abline(h=MEDissThres, col = "red")
今回はmodule eigengen間の距離0.15を基準にモジュールのマージを行ってみました！
マージ前のものと比較してみます！
merge <- mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
mergedColors <- merge$colors
mergedMEs <- merge$newMEs
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
結果は上のようになりました！今回はマージした後もGene ontologyやパスウェイ解析の結果を見ても大きな変化がなかったのでマージ前のモジュールを採用することにしました！
moduleColors <- dynamicColors
colorOrder <- c("grey", standardColors(50))
moduleLabels <- match(moduleColors, colorOrder)-
save(MEs, moduleLabels, moduleColors, geneTree, file = paste0("../Output/003_ModuleDetection/SLE_GSE138458-03-ModuleDetectionn", "_", "signed", "_",deepSplit,".RData"))今回の結果を保存してモジュール検出は終了です！
今回もすごく長くなってしまいました笑 ここまで読んでいただいて本当にありがとうございました！
次回は得られたモジュールとTrait情報（今回は病態スコア）と関連付けて解析していきます！