WGCNA

解析の流れ

重み付き相関行列

パッケージとデータの読み込み

RでWGCNAパッケージのインストールを行います。インストールはBiocManagerを使って行うので、BiocManagerがインストールされていない場合はこちらのインストールします。

## BiocManagerがインストールされていない場合
#install.packages("BiocManager")

### WGCNAのインストール
BiocManager::install("WGCNA")

インストールしたパッケージを読み込みます。

library(WGCNA)

WGCNAパッケージの他に使用するパッケージも読み込みます。それぞれ、インストールしていない場合は、install.packages("")を使ってCRANからインストールします。

library(igraph) #グラフの計算
library(visNetwork) #グラフの可視化
library(mixOmics) #TGCAデータを取得するため
library(colorBlindness) #カラーパレットの取得
library(tidyverse) #ggplot,dplyrなどをまとめて読み込み

次に、mixOmicsパッケージに含まれているマルチオミクスデータを読み込みます。TCGAは乳がん患者をサンプルとしたマルチオミクスデータで、mixOmicsにはTCGAの抜粋データが入っています。読み込んだデータは、150人の乳がん患者をサンプルに、miRNA（184）、mRNA（２００）、タンパク質（142）の３つのオミクスデータを含む、150×（184+200+142）のデータ行列です。サンプルの患者は、Basal、Her2、LumAという３種類の乳がんのサブタイプを患った患者に分類されます。乳がんはサブタイプごとに特徴が異なり、サブタイプを詳しく調べることで治療法の発見につながる可能性があります。

data('breast.TCGA')
data <- cbind(breast.TCGA$data.train$mrna,
             breast.TCGA$data.train$mirna, 
             breast.TCGA$data.train$protein) %>% #パイプ演算子
        data.frame()

soft-Thresholding　powerの設定

まずは、相関行列の重み付けに用いるsoft-Thresholding　powerの値を、ネットワークのスケールフリー性を基準に決定します。スケールフリー性とは、ごく少数のノードがたくさんのエッジを持っているようなネットワークの性質を言い、次数分布がべき乗で表現されます。スケールフリーなネットワークでは、次数とその頻度（確率）の両対数が直線上に分布するので、ここでは回帰直線の決定係数（R^2)が、設定した閾値よりも高くなるような"最小"の値（大きくしすぎるとネットワークのコネクティビティが失われてしまうため）をsoft-Thresholding　powerに採用します。

WGCNAのpickSoftThreshold関数を用いると、上記の観点からsoft-Thresholding　powerの値を選択します。引数には、サンプル×変数の行列データ(data)、soft-Thresholding　powerの候補となる値を入れたベクトル（powerVector）、ネットワークの種類（networkType）、フィッテイングの閾値(RsquaredCut)を渡します。ここでは、無向グラフで、R^2が0.8以上になるようにsoft-Thresholding　powerを設定します。

fit.index <- 0.8 #0.8以上推奨
powers  =  c(1:10)

sft <- pickSoftThreshold(data, powerVector=powers, networkType = "unsigned", RsquaredCut = fit.index, verbose=5)

soft-Thresholding　powerの値とR^2の値をプロットして確認します。今回のデータでは、５以上でR^2が0.8を超えました。soft-Thresholding　powerの増加に伴い、ネットワークの連結性が急速に減少することがわかります。

par(mfrow		=		c(1,2))
cex1		=		0.8

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=cex1)

abline(h=0.8, col = "red")


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)

隣接行列の作成

選択したsoft-Thresholding　powerを使って、重み付きの隣接行列を計算します。WGCNAのadjacency関数を使います。選択したsoft-Thresholding　powerは格納した変数のsft＄powerEstimateに入っています。行名・列名がリセットされるので入れ直します。

Adj <- adjacency(data, power = sft$powerEstimate)

id <- colnames(data)
colnames(Adj) <- id
rownames(Adj) <- id

ネットワークの可視化のため、値の小さなエッジを枝刈りします。ここでは閾値を0.1に設定して枝刈りを行いました。

thre <- 0.1
subAdj <- (Adj>thre)*Adj

次数分布のヒストグラムと両対数グラフを出して確認します。R^2の値が0.97で、スケールフリー性の高いネットワークが構築されたことが確認できます。

k <- subAdj %>%
  apply(2, sum, na.rm=T) %>%
  as.vector()

par(mfrow=c(1,2))
hist(k)
scaleFreePlot(k, main="Check scale free topology\n")

モジュールの計算

次に、構築したネットワーク内のモジュールを、igraphパッケージを使って計算します。

まずは、igraphのgraph.adjacencyで、隣接行列をネットワーク作成用のデータフォーマットに変換します。引数で、無向（mode = "undirected"）、重みあり（weighted = TRUE）、自分自身に向かうエッジはなし（diag = FALSE）を設定します。

subnet <- graph.adjacency(subAdj, mode = "undirected", weighted = TRUE, diag = FALSE)

次に、モジュラリティ最適化の手法を用いてモジュールを計算します。モジュラリティの計算には、Fast Greedy Algorithmを用いて計算します。igraphのfastgreedy.community関数にsubnetを渡すと、自動で計算が行われます。モジュールの結果を変数moに格納し、10以上のノードを含むモジュールをmo.sに格納しました。

m <- fastgreedy.community(subnet)

mo <- m$membership %>%
  table() %>%
  data.frame()

colnames(mo) <- c("module.ID", "No.of.node")

mo <- mo[mo$No.of.node >=1, ]
mo.s <- mo[mo$No.of.node >= 10, ]
#mo.s
#sum(mo.s$No.of.node)

ネットワークのトポロジーの計算

構築されたネットワークの各変数の、①degree（次数）、②strength（強度）、③betweenness（媒介中心性）、④page rank、⑤モジュール番号をまとめます。

degree <- degree(subnet)
strength <- strength(subnet)
between <- betweenness(subnet, normalized = TRUE)
page.rank <- page_rank(subnet)$vector
module <- m$membership
result_adj <- cbind(degree, strength, between, page.rank, module) %>%
  data.frame()

ネットワークの可視化

# visualization by igraph
V(subnet)$color <- "#999"
#V(subnet)

for (i in c(1:length(mo.s$module.ID))) {
  V(subnet)[m$membership == mo.s$module.ID[i]]$color <- colorBlindness::safeColors[i]
  V(subnet)[m$membership == mo.s$module.ID[i]]$group <- paste("Module.", i, sep = "")
}

n <- induced_subgraph(subnet, m$membership %in% mo.s$module.ID) 

set.seed(108)
# visualization by visNetwork
netv <- toVisNetworkData(n)
# construct interactive network
visNetwork(netv$nodes, netv$edges) %>%
  visIgraphLayout(layout = "layout_with_fr") %>%
  visNodes(shadow = list(enabled = TRUE ,size=20)) %>%
  visOptions(highlightNearest = list(enabled = T, hover = T),
             selectedBy = "group",
             nodesIdSelection =TRUE)

Topological Overlap Matrix

これまでは重み付きの相関行列を基にネットワークを構築しましたが、次はTopological Overlap Matrixを用いて同様にネットワークを構築します。Topological Overlap Matrixは、WGCNAのTOMsimilarity関数に隣接行列を渡して計算します。同様のp×p行列が出力されます。

TOM <- TOMsimilarity(Adj)

colnames(TOM) <- id
rownames(TOM) <- id

以下、同じようにネットワークを構築して可視化します。

thre <- 0.07 # エッジの枝狩り
subTOM <- (TOM>thre)*TOM 
subnet <- graph.adjacency(subTOM, mode = "undirected", weighted = TRUE, diag = FALSE)

k <- subTOM %>%
  apply(2, sum, na.rm=T) %>%
  as.vector()

par(mfrow=c(1,2))
hist(k)
scaleFreePlot(k, main="Check scale free topology\n")

m <- fastgreedy.community(subnet)
mo <- m$membership %>%
  table() %>%
  data.frame()
colnames(mo) <- c("module.ID", "No.of.node")
mo <- mo[mo$No.of.node >=1, ]
mo.s <- mo[mo$No.of.node >= 10, ]
#mo.s
#sum(mo.s$No.of.node)

# network property
degree <- degree(subnet)
strength <- strength(subnet)
between <- betweenness(subnet, normalized = TRUE)
page.rank <- page_rank(subnet)$vector
module <- m$membership
result <- cbind(degree, strength, between, page.rank, module) %>%
  data.frame()

# visualization by igraph
V(subnet)$color <- "#999"
#V(subnet)

for (i in c(1:length(mo.s$module.ID))) {
  V(subnet)[m$membership == mo.s$module.ID[i]]$color <- colorBlindness::safeColors[i]
  V(subnet)[m$membership == mo.s$module.ID[i]]$group <- paste("Module.", i, sep = "")
}

n <- induced_subgraph(subnet, m$membership %in% mo.s$module.ID) 

set.seed(108)
# visualization by visNetwork
netv <- toVisNetworkData(n)
# construct interactive network
visNetwork(netv$nodes, netv$edges) %>%
  visIgraphLayout(layout = "layout_with_fr") %>%
  visNodes(shadow = list(enabled = TRUE ,size=20)) %>%
  visOptions(highlightNearest = list(enabled = T, hover = T),
             selectedBy = "group",
             nodesIdSelection =TRUE)