先日WGCNA（Weighted Gene Coexpression Network Analysis）を試しに使ってみました（この記事 → 論文のメモ: 生態毒性研究へのWGCNAの適用）。しかし、遺伝子発現のデータは別に直線関係ばかりではないので（最近読んだ論文のメモ： omicsデータの用量応答モデル）、ピアソンの相関係数に基づいて計算することに違和感を覚えました。
そこで、非線形の関係も扱える相互情報量（mutual information）に基づくネットワーク分析を勉強してみました。
参考にした情報は以下の論文＋Wikipediaの"相互情報量"と"情報量"。
ARACNEをRで実装した論文。
Meyer PE, Lafitte F, Bontempi G, 2008, minet: AR/Bioconductor package for inferring large transcriptional networks using mutual information, BMC Bioinfor 9(1): 461.
相互情報量とピアソン相関係数のどちらが適切かは時と場合によるという話。
法隆大輔, 林武司, 2013, 遺伝子発現プロファイル類似度としてのピアソン相関係数と相互情報量の比較, 計量生物学 33(2):125-143.
まず相互情報量を求める
相互情報量の定義はwiki参照。数式で書くとこんな感じらしい。
(wikiより)連続値の場合の周辺確率と同時確率の求め方は、色々提案されているようですが（たぶん）、一番分かりやすいのはデータを分割して離散化する方法でしょう。
下のデータで遊んでみる。G1とG2が遺伝子発現量を、Sが試料番号を示してます。
> data
S1 S2 S3 S4 S5 S6 S7 S8 S9
G1 62 41 5 30 15 52 45 1 39
G2 167 42 136 50 98 114 90 167 38
ピアソンの相関係数は-0.20です。
一方、相互情報量の場合。
まずデータを分割して離散化します。区切り方は、区切りの幅を同一にする"equalwidth"か区切り内のデータ数を同一にする"equalfreq"か、など選択します。
> d <- discretize ( t(data) , nbins=3 ) #デフォルトだとnbinsはサンプルサイズの平方根
> d
G1 G2
1 3 3
2 2 1
3 1 3
4 2 1
5 1 2
6 3 2
7 3 2
8 1 3
9 2 1
## こういう感じに分割される
# このデータの場合equalfreqもequalwidthも同じ
上の図をもとに集計表を作成し、X（=G1）とY（=G2）のエントロピーを求めます。
- H(X) = - { 1/3*log(1/3) + 1/3*log(1/3) + 1/3*log(1/3) } = 1.584963
- H(Y) = - { 1/3*log(1/3) + 1/3*log(1/3) + 1/3*log(1/3) } = 1.584963
- H(X,Y) = - { 1/3*log(1/3) + 2*1/9*log(1/9) + 2*2/9*log(2/9) } = 2.19716
- I (X; Y) = 1.584963 + 1.584963 - 2.19716 = 0.972766
一応、Rのinfotheoライブラリーのmutinformation関数で検算しておきます。求め方は"empirical (=maximum likelihood)" か"shrink entropy"か、"schurmann-grassberger"など選べます。
empirical methodで求めた相互情報量は0.973で、上記の手計算と一致していますね。
> mutinformation(d, method="emp") ## 底はe
G1 G2
G1 1.099 0.674
G2 0.674 1.099> natstobits ( mutinformation(d, method="emp") ) ## 底2のビットに変換
G1 G2
G1 1.5849625 0.9727653
G2 0.9727653 1.5849625
ちなみにRのminetパッケージ（infotheoを利用している）だとbuild.mim関数で一気に求められます。
> build.mim( t(data), estimator="mi.empirical", disc="equalfreq")
G1 G2
G1 0.0000000 0.6742695
G2 0.6742695 0.0000000
相互情報量をもとにネットワーク分析をしてみる
ぶっちゃけRのminetライブラリにぶち込むだけですが、とりあえずやってみます。
次はデータを拡張して、下のものを使用します。G1~G10が遺伝子で、各サンプルS1~S6における発現量のデータです。サンプル6個で相互情報量を求めるのは厳しいかも…。
ピアソン相関係数の行列はこんな感じ。
相互情報量を、build.mim関数で求めたのが下。
> mim.eqwi <- build.mim( t(data) , estimator ="mi.empirical", disc="equalwidth")
このように類似度行列が得られたら、次は図示するための準備。WGCNAの場合と同様に、閾値で区切って0/1データに変換するなどの処理を行います。とりあえずminetライブラリに実装されているARACNEを用います。
> net.eqwi <- aracne(mim.eqwi)
#赤字の場所がARACNEによってカットされた部分。indirectな関係を除外できるということですが…？
最後に、相互情報量と相関係数の結果を比較。snaライブラリーで図示。
> mim.cor <- build.mim( t(data) , estimator ="pearson") # 相関係数使用, 離散化はされない
> net.cor <- aracne(mim.cor)
> library( sna )
> par (mfrow=c(2,1) )
> gplot(net.eqwi, usearrows=FALSE, displaylabels=TRUE)
> gplot(net.cor, usearrows=FALSE, displaylabels=TRUE)
