ネットワーク解析wgcna

WGCNA解析 part 1

2025/11/02
#deg#gsva#rna-seq#wgcna#エンリッチメント解析#ネットワーク解析#モジュール検出#遺伝子優先順位付け

普段のRNA-seq等の遺伝子発現解析をだだのDEG解析から一歩進んだ解析としてWGCNA解析を紹介したいと思います。

コードは載せず、WGCNAの流れが分かるような内容にします。
コードありバージョンはいつか何らかの形で出します。

WGCNAって？

WGCNA (Weighted Gene Correlation Network Analysis)とは大量の遺伝子発現データから「一緒に上下する（共発現）」遺伝子群＝モジュールを見つけ、モジュールと表現型（群・連続値など）の関連を調べ、ハブ遺伝子（モジュールの中心）を仮説として挙げる手法です。因果ではなく相関ベースの探索です。

流れ（ざっくり）

  1. 前処理：正規化・低発現除去・バッチ補正、外れ値サンプル除去

  2. ネットワーク：相関を**softPower(β)**で重み付け（近いほど強い）→TOM（トポロジカル重み）→階層クラスタ→モジュール抽出

  3. 要約：各モジュールの第1主成分＝モジュール固有ベクトル（ME）

  4. 解釈：MEと形質の相関で“どのモジュールがどの表現型と動くか”を見る

  5. 優先度付け：kME（モジュール所属度）やkWithin（モジュール内連結度）、GS（Gene Significance：遺伝子と形質の相関）でハブ候補を選ぶ（高kMEかつ高GSが狙い目）

今回扱うデータに関して

今回使うデータセットは「GSE186294」にしました。

被験者は4週間にわたり、2日おきに体重1kgあたり50IUの組換えヒトエリスロポエチン（rHuEPO）注射を受けた。rHuEPO投与期間中（4週間）、毎日経口鉄剤（元素鉄100mg、硫酸第一鉄錠、Almus社製、英国バーンスタブル）が投与された。ここで、ベースライン時（Base1およびBase2）、rHuEPO投与期間中（初回投与後2日目および2週間後；EPO3およびEPO4）、ならびにrHuEPO投与終了後2週間（Post7）に採取した全血サンプルを、イルミナNextSeq 500を用いて解析した。

Base 1 → Base 2 → EPO3 → EPO4 → Post 7 という時系列のサンプルです。

WGCNA解析実行結果

1. softPower(β)の設定

softPower(β)とは？

まずは**softPower(β)**を設定します。

各遺伝子の相関cor(i, j)を計算した後、この相関に掛ける重****みのことです。
βを上げるほど、強い相関だけが残りやすく（スパース化）、弱い相関は0に近づきます。

unsigned:  aij=∣cor(i,j)∣βa_{ij} = |cor(i,j)|^{\beta}

signed:  aij=(1+cor(i,j)2)βa_{ij} = \left(\frac{1+cor(i,j)}{2}\right)^{\beta}（負相関は小さくなる）

この重みを使って TOM → 階層クラスタ → モジュール が決まります。

βはどうやって決める？

それではβを決めていきます。

WGCNAでは、遺伝子ノード次数分布がスケールフリー（べき乗則）に近くなるようにβを定める。

WGCNAのpickSoftThreshold()関数を実行し複数のβを試し、以下の2つを検討する：

  1. Scale-free fit（R²）：べき乗則にどれだけフィットするか

  2. Mean connectivity：平均連結度（βを上げると下がる）

基本ルール：

  • 最小のβでR² ≳ 0.8（小規模やノイジーなら 0.7 も妥協点）

  • かつ mean connectivity が極端に低くなり過ぎない地点を選ぶ

今回はβ=7 が選択されました。

系統樹の確認

TOMij=∑uaiuauj+aijmin⁡(ki,kj)+1−aij,  ki=∑uaiu\mathrm{TOM}_{ij}=\frac{\sum_u a_{iu}a_{uj} + a_{ij}}{\min(k_i,k_j)+1-a_{ij}},\ \ k_i=\sum_u a_{iu}

aija_{ij}は既出ですので省きます。

ki=∑uaiuk_{i}=\sum_u a_{iu}はノード(遺伝子)iの重み付き次数

TOMの分子の左側 ∑uaiu auj\sum_{u} a_{iu}\,a_{uj} は、
遺伝子iと遺伝子jがその他の遺伝子とどれぐらい共通で繋がっているかの指標です。

次は分母の説明ですが、
min⁡(ki,kj)+1−aij\min(k_i,k_j) + 1 - a_{ij}
これは理論上取りうる最大の重なりで割って 0 - 1に正規化しています。

dissTOMij=1−TOMij\mathrm{dissTOM}_{ij}=1-\mathrm{TOM}_{ij}
は距離に変換しています。TOMが近さだから、逆にすることで遠さ=距離に変換しているわけです。

これをクラスタリングや樹形図を作る際には距離を使用するためこの変換が必要です。

その結果が以下の樹形図です。

また各モジュールのサイズ（遺伝子数）は以下の通り

また、各モジュールのmodule eigengenes (ME)を算出し、モジュール間の類似性を比較しています。

MEとは各モジュールの遺伝子発現を主成分分析したときの第1主成分＝サンプル方向の代表発現（モジュール活動度）
1モジュール1サンプルにつき1つの値になる

着目モジュールの選定

次は、作成したモジュールのうち、どのモジュールに着目するかを決めていきます

どのように着目するモジュールを選んでいくかには、都度変わってくるかと思いますが、今回は各モジュールと各処置との関係から選んでいこうと思います。

各サンプルの各モジュールの「値」が処置(Base, EPO, Postなど)に伴って変動するどうかを検討すれば良いので、サンプルラベルと各サンプルのモジュールの「値」の相関を見ていきます。

この各サンプルのモジュールの「値」にMEを使います。

その結果が以下です。

まずは、EPO4と高い相関が見られるModule 3に着目していきたいと思います。

今回はココまで、次回からModule 3について詳しく解析していきたいと思います。

関連記事WGCNA解析 part 2 →