前回のPart1ではrawデータのダウンロードから正規化を行い、サンプルと遺伝子フィルタリングまでを行いました。
Weighted Gene Coexpression Network Analysis (WGCNA)は遺伝子発現量の相関を利用して、互いによく似た発現変動をする遺伝子同士のクラスターを発見する手法です！ WGCNAは単純にクラスター[…]
今回はいよいよメインのWGCNAによるネットワーク構築をしていきます！
このPartのコードは以下のチュートリアルを参考にしています。
WGCNAのネットワーク構築の流れ
WGCNAによるネットワーク構築の流れは以下です。
- Soft Thresholdの選択
- Scale Free Topologyの確認
- 階層的クラスタリングによるモジュール検出
最初は単語一つ一つが何を意味しているのかを捉えるのが難しいと思いますが、適宜説明していきますので今はそういうものがあるんだくらいで読み進めてください！
まず必要なライブラリーと前回作成した前処理済みの発現データ、Traitデータをロードします。
library(limma)
library(gplots)
library(tidyverse)
library(WGCNA)
lnames <- load("../Output/001_NormalizeRawData/SLE_GSE138458-01-dataInput.RData")
lnames
WGCNAでは遺伝子の発現プロファイルが互いにどれくらい似ているのか？を指標にして遺伝子同士をつなぎ合わせていきネットワークを構成します。（遺伝子を頂点として、似ている遺伝子同士は辺でつないだネットワーク）
その発現プロファイルの類似度にはその遺伝子同士の相関を用います。
ここで問題になるのがどれくらい相関が強ければその遺伝子の間に辺を持たせるべきなのかという基準です。
シンプルな考え方だとある値以上の相関係数（例えば0.8以上）があればその2つの遺伝子はネットワーク上でつながっているいるとするというものがあります。（Hard thresholdingといいます）
$$
\begin{eqnarray}
a_{i,j}=\left\{ \begin{array}{ll}
1 & (s_{ij}\geq0.8) \\
0 & (otherwise) \
\end{array} \right.
\end{eqnarray}
$$
しかしこの基準だと相関係数が0.80だとリンクを持っていて、0.79だとリンクを持っていないとすることになります。相関係数で0.01の違いをこれほどドラスティックにしてしまってよいのでしょうか？（実際はこのUnweighted Networkもよく使われますが）
この基準では明らかに情報をロスしてしまっています。また、共発現ネットワークにおいて繋がりの強度を連続的に表現できるもののほうが、生物学的に自然なネットワークを反映できているとも考えられます。
そこで遺伝子間のつながりを連続値で表現するために用いられるのがsoft thresholdです！
上記の議論はBMC Bioinformatics, 2008を参考にしています
Soft Thresholdの選択
WGCNAでは遺伝子間のつながりを以下のように計算します。
$$a_{i,j}=s_{i,j}^\beta$$
このように遺伝子間のつながりが連続値の重みをもつものをWeighted Networkと呼び、上述のような遺伝子間のつながりを0 or 1で表現するUnweighted Networkと対比されます。
Soft Thresholdとは$a_{i,j}=s_{i,j}^\beta$の$\beta$のことを指します。
なんのためにSoft Thresholdを使うかと言うと、弱いリンクをさらに弱くすることで強いリンクを際立たせるために使います！
例えば遺伝子1と2が相関係数0.8で遺伝子1と3が相関係数0.2だったとして、$Soft Threshold=5$を使ったとすると遺伝子間のリンク$a_{1,2}$ と$a_{1,3}$はそれぞれ
$$a_{1,2}=0.8^5=0.32$$
$$a_{1,3}=0.2^5=0.00032$$
となり$a_{1,2}$ のリンクは$a_{1,3}$のリンクの1000倍も強いことになり、もとの値より差が際立っています。
WGCNAパッケージではpickSoftThreshold関数で最適なsoft thresholdを探索できます。
また、上述のHard thresholdingを行いたい場合はpickHardThreshold関数を用意してくれています。
結論から言うと、「ネットワークがScale Free Topologyを持つことができるSoft Thresholdの中で最小のもの」です。
Scale Free Topologyってなんぞやって話ですよね(笑)
本当にざっくり説明するとScale Freeなネットワークとは「一部のノード（頂点のこと：今回だと遺伝子）にリンクが集中しているネットワーク」のことです。
現実の事例をあげるとインターネットでの各サイトの繋がり方はまさにScale Freeなネットワークで、GoogleやAmazonなど一部は莫大なリンクを持っていますが、個人のブログなどの大部分のサイトは少数のリンクしか持っていません。
transcription networkやprotein interaction networkなど生物学的なネットワークもScale Free性を持っているとされています(Journal of Cell Science, 2005)のでWGCNAではこれを再現するようにネットワークを構成したいというわけです！
Scale free なネットワークの定義についてもう少し厳密に話すと次数（リンクの数のこと）$k$を持つノードの割合が$k^{-\gamma}$に比例するようなネットワークのことを言います。式で書くと以下です。
$$P(k)\sim k^{-\gamma}$$
$P(k)$がネットワーク中のすべてのノードの中で$k$本のリンクを持つものの割合です。その割合が$k$が大きくなるほどに指数関数的に減少していくような分布を示すネットワークのことです。
つまりScale Freeなネットワークかどうかは次数$k$のヒストグラムを両対数グラフで描画したとき分布が$-\gamma$に比例して減少しているかどうかを確認すればいいことになります。
Soft Thresholdを大きくすると基本的にはScale Free Topologyを持つネットワークを得られるのですが、あまりに大きすぎると生物学的に真に関連のある遺伝子同士のつながりを見逃してしまうかもしれないので、「ネットワークがScale Free Topologyを持つことができるSoft Thresholdの中で最小のもの」がベストだと考えられます。
これを以下で確認していきます。
Scale free natworkについてはwikipediaを参考にしています。
Soft Thresholdの決定
WGCNAパッケージのpickSoftThreshold関数で最適なSoft Thresholdを探します！
powers = c(c(1:10), seq(from = 12, to=30, by=2))
sft = pickSoftThreshold(datExpr, powerVector = powers, networkType = "signed", verbose = 5)
pickSoftThreshold関数のアウトプットのsftの中身は上のようになっています。
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=0.9,col="red")
abline(h=0.90,col="red")
このプロットは各Soft ThresholdでのネットワークのScale Free性を表しています！
横軸はSoft Thresholdで縦軸は次数の分布の回帰直線の$R^2$値です！縦軸は「次数の分布の両対数が直線で近似できるか？」を指標としてScale Free性を確認しようとしています！
縦軸についてもう少し詳しく話すと、各Soft Thresholdで遺伝子間のつながりの強度を計算した際の以下のようなプロット（横軸がネットワーク中のノードの次数、縦軸がその次数を持つノードの割合の散布図を両対数取ったもの）の$R^2$値です。このプロットの描画方法は後述の遺伝子間の非類似度の計算で出てきます。
上で ”Scale Freeなネットワークかどうかは次数$k$のヒストグラムを両対数グラフで描画したとき分布が$-\gamma$に比例して減少しているかどうかを確認すればいいことになります” と書きましたが、それがまさにこのプロットになります。
このプロットで高い次数を持つノードの数が指数関数的に減少している（両対数なのでグラフの中では線形に減少している）ことを回帰直線を引いてその$R^2$値で確認しているわけです。
上のように$R^2$値で0.85とある程度線形に高次数のノードが減少していればそのネットワークはScale Free性を持っていると考えられます。
さて、以下のプロットに話を戻します。
WGCNAパッケージの作者らは$R^2$値で0.9をScale Free ネットワークの一つの基準として提案しています。上のプロットの赤い直線がそれです。
今回のデータで最初はこの基準を満たすSoft Threshold 22を選択してネットワークを構築してみましたが、ほとんどの遺伝子がどのモジュールにも割り当てられないという結果となってしまいました。
そこである程度$R^2$値 0.9の基準に近く、値の上昇が飽和し始めているSoft Threshold 10を用いて再びネットワーク構築を行ったところよい結果が得られたため今回はSoft Threshold 10を採用することにしました。
Soft Thresholdごとのネットワークの平均次数も確認しておきます。
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")
平均次数はSoft Threshold 6程度まで減少し、Soft Threshold 10でほぼ一定になるようです。
遺伝子間の非類似度の計算
さて、Soft Thresholdも決まったのでこれを用いて遺伝子間の”非”類似度を計算していきます。
まずAdjacency Matrixを計算し、さらにノイズなどの影響を小さくするためにAdjacency MatrixをTopological Overlap Matrix（TOM）に変換します。
このTOMはあとで詳しく説明しますが、ざっくり言うと「ネットワークの中で2つのノードがどれくらい似ているのか？」の指標です。値は０〜１の範囲で、大きいほど2つのノードが似ていることを示します。
WGCNAでは遺伝子間の”非”類似度を用いて階層的クラスタリングを行うため遺伝子の”類似度”の指標であるTOMを”非類似度”の指標であるdissimilarity TOMに変換します。これはTopological Overlapを指標にして、似ている遺伝子同士を近くにクラスタリングするためです。
softPower = 10
adjacency = adjacency(datExpr, power = softPower, type="signed")
TOM = TOMsimilarity(adjacency,TOMType = "signed")
dissTOM = 1-TOM
Scale Free性の確認
選択したsoft threshold 10を用いて非類似度の計算ができました。ここで作成したネットワークがスケールフリー性を持っているかどうかを確認しておきましょう！
次数のヒストグラムとscaleFreePlot関数で確認します！
これらは上でも述べたようにScale Free性を「次数分布の両対数が直線で近似できるか？」を指標として確認するためのものです！
まずネットワーク中のノードの次数の分布をヒストグラムで確認してみます。
k=as.vector(apply(adjacency, 2, sum, na.rm=T))
hist(k)次数$k$が増えるほど指数関数的に減少していそうですね。
今度はWGCNAパッケージのscaleFreePlot関数で次数$k$が$P(k)\sim k^{-\gamma}$の分布になっているかを確認してみます！
scaleFreePlot(k, main="Check scale free topology\n")
$R^2$値が0.85で回帰直線の傾きが$-1.7$なのでおおよそ$P(k)\sim k^{-\1.7}$の分布になっていて、ある程度Scale Free性があるネットワークが得られたと考えても良さそうです。
TOMについて
Topological Overlap Matrix（TOM）ですが、具体的に何を意味しているかと言うと、ネットワーク中の2つのノードを比較して、「リンクしている先がどれくらい重複しているか？」です。読んで字のごとくで、ノードのトポロジー（繋がり方）がどれくらいオーバーラップ（重複）しているのかということです。
ノード$i$とノード$j$についてのTOMの計算式は以下です。(WGCNAパッケージTOMsimilarityFromExpr関数ヘルプ参照)
$$TOM_{i,j}=\frac{a_{i,j}+\sum_{u}^n a_{i,u}a_{u,j}}{min(k_i,k_j)+1-a_{i,j}}$$
落ち着いてみれば数式の解釈も難しくはないです。理解を簡単にするためにノード間のリンクに重み考えない（つながってたら１, つながってなかったら０）Unweighted Networkで考えてみます。
まず分母は$min(k_i,k_j)$の部分で$i$と$j$で次数kの小さい方をとってきます。次に$i$と$j$の間にリンクがあれば（$a_{i,j}=1$)0を、$i$と$j$の間にリンクがなければ（$a_{i,j}=1$)1を足します。
最後の$+1-a_{i,j}$は$i$と$j$の間のリンクの有るとき$min(k_i,k_j)$の値が１大きくなってしまうのを調整するためのものなのであまり気にする必要はありません。
次に分子を見てみます。$\sum_{u}^n a_{i,u}a_{u,j}$はノード$u$に対して$i$も$j$も両方がリンクを持つときに+1され、どちらかが$u$にリンクを持たないときは+0されます。つまり共通するリンク先の数を計算しています。これに$i$と$j$同士のリンクの有無を加えています。
まとめると (重複しているリンク先の数) / (重複するリンクの理論上最大値) を計算してることがわかると思います。
上でも述べたように文字通りノードのトポロジー（繋がり方）がどれくらいオーバーラップ（重複）しているかを表していますね！
Unweighted Networkでなくリンクに重みをもたせたWighted Networkでも同様に考えられます！
signed unsignedについて
TOMsimilarity関数の最も大事な引数はTOMTypeです。 “signed”, “unsigned”などの選択肢があります。”unsigned”はネットワークのリンクの重みの記号（プラスかマイナスか）を考慮しないという意味です。つまり遺伝子の発現が正相関しようが逆相関しようが関係なくリンクには正の値を持たせます。逆に”signed”は正相関と逆相関を区別して、遺伝子同士が逆相関するならその遺伝子間のリンクはマイナスの重みをもたせます。
これは非常に重要な選択肢で得られる結果がまるで違うものになるので慎重に検討しましょう。
また、WGCNA FAQで議論されていますが、一般的にはSignedネットワークのほうが好ましいとWGCNAパッケージの作成者によって言及されています。しかし、今回用いたadjacency()とTOMsimilarity()のデフォルトはUnsignedです。Signedネットワークを構築したいときには今回のようにadjacency(datExpr, type=”signed”) TOMsimilarity(adjacency, TOMType=”signed”)というように引数で指定してあげてください！
階層的クラスタリング
最後に上で作成した非類似度を用いて階層的クラスタリングを行います。
geneTree = hclust(as.dist(dissTOM), method = "average")
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04)
上のようなデンドログラムを作ることができました！
ここまでの結果を保存しておきます！
save(adjacency, dissTOM, geneTree,
file = '../Output/002_NetworkConstruction/SLE_GSE138458-02-networkConstruction.RData')
長くなってしまったので今回の記事はここまでにして次回、このデンドログラムを用いてモジュール(遺伝子のクラスター)検出をやっていきます！
Part1では正規化などのデータの前処理、Part2ではネットワークの構築とデンドログラムの作成を行ってきました！ [sitecard subtitle=関連記事 url=https://shiokoji11235.com[…]
参考