Weighted Gene Coexpression Network Analysis (WGCNA)は遺伝子発現量の相関を利用して、互いによく似た発現変動をする遺伝子同士のクラスターを発見する手法です！
WGCNAは単純にクラスターを見つけるだけではありません！
それぞれのクラスターの中で中心的な役割を果たす遺伝子(hub gene)を定量的に見つけたり、メタデータ（体重や病気の進行具合のスコアなど）と遺伝子クラスターとの相関を解析することができるめちゃくちゃ強力な手法です！
実際に、WGCNAでバイマーカーや創薬標的の探索を行なった論文がいくつも報告されています！
今回はそんな強力な手法であるWGCNAを使って実際にマイクロアレイのデータ解析をやってみました！
解析の概要
チュートリアルをそのままやるだけでは面白味に欠けると思ったのでGEOから全身性エリテマトーデス(SLE)のマイクロアレイデータを取ってきました！
使用したマイクロアレイのデータはGSE138458です！このデータセットを選んだ理由は2つです！
- サンプル数が十分にある（Healthy control 48サンプル、SLE 100サンプル）
- SLEの病態スコアであるSELENA-SLEDAIの情報を使用できる（←めちゃくちゃ重要）
WGCNAにはそれなりに多くのサンプル数（各条件で20サンプル以上）が必要になるため、サンプル数が十分にあることがデータセットの必要条件になります！その点、今回使用したデータセットは全部で148サンプルと十分です！
また、WGCNAはメタデータ（体重や病気の進行具合のスコアなど）を用いてそれらと遺伝子クラスターとの相関を解析することができるのが最大の強みだと私は思っています！そこで全身性エリテマトーデスの病気の進行度合いを数値化したSELENA-SLEDAIもくっついてくるこのデータセットはWGCNAの威力を最大限活かせると考えてこのデータセットを選択しました！
解析の流れはZhu et al. (2019)を参考にしました！
この記事ではデータの前処理を行いますがもし正規化済みのデータをお持ちでしたらこの記事を飛ばしてPart2 ネットワーク構築からご覧になってください！
前回のPart1ではrawデータのダウンロードから正規化を行い、サンプルと遺伝子フィルタリングまでを行いました。 [sitecard subtitle=関連記事 url= https://shiokoji11235.com[…]
発現データの前処理
まずGEOから正規化前のデータを取得します
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE138458
上のリンクの下の方にある GSE138458_non-normalized.txt.gz をダウンロードします
ダウンロードできたら適当なディレクトリの中で解凍してください。
私の場合はDataディレクトリの中で解凍しました。
以降、/Data/GSE138458_non-normalized.txt として使用します
さて、このデータを前処理していきましょう！
必要なライブラリをインポートします。
library(gplots)
library(tidyverse, warn.conflicts=F, quietly=T)
library(WGCNA)
library(GEOquery)
library(RColorBrewer)
library(dendextend)
library(ggdendro)
read_tsv("../Data/GSE138458_non-normalized.txt")%>%
rename_at(vars(starts_with("200")), function(x){paste0("SAMPLE",x)}) %>% # sampleの発現量が記載された列の列名に'SAMPLE'を追加
write_tsv("../Data/GSE138458_non-normalized_renamed.txt") #ファイル名を変更して保存
data <- read.ilmn("../Data/GSE138458_non-normalized_renamed.txt", sep = "\t", probeid ="ID_REF", expr = "SAMPLE") # ファイル読み込み
read.ilmn関数でデータを読み込みました！ この関数は発現量を記載している列名の目印となる文字列を必要とします。なので元のデータの列名を加工して’SAMPLE’という文字列を発現量の列名に追加しています。
メタデータの取得と加工
SLEスコアなどのメタデータもダウンロードして少し加工してTraitsとして保存しましょう！
gse <- getGEO("GSE138458", destdir = getwd())
phenoData <- pData(gse[[1]])
sledai_group_map <- c(0, 1, 2)
names(sledai_group_map) <- c("no val", "Low", "High")
case_control_map <- c(0, 1)
names(case_control_map) <- c("Control", "SLE Case")
Traits <- phenoData %>%
rownames_to_column("ID_REF") %>%
dplyr::select("ID_REF", "description", "case/control:ch1", "sledai group:ch1") %>%
dplyr::rename(case_control="case/control:ch1",
sledai_group="sledai group:ch1") %>%
mutate(sledai_group_levels=sledai_group_map[sledai_group],
case_control_levels=case_control_map[case_control])
Traits %>%
write_csv("../Data/Traits.csv")
以下のようなデータフレームを作成してTraits.csvとして保存しました！
読み込んだ発現データの確認
これからサンプルと遺伝子のフィルタリング（外れ値の除去）を行っていきます。そこでまずはこれから扱うことになる発現データの中身を少しだけ見てみましょう。
まずは各サンプルごとに全遺伝子の発現量の平均値を見てみましょう
meanExpressionByArray=apply(data$E, 2, mean, na.rm=T)
sizeGrWindow(9, 5)
barplot(meanExpressionByArray,
xlab = "Sample", ylab = "Mean expression",
main ="Mean expression across samples",
names.arg = c(1:dim(data$E)[[2]]), cex.names = 0.7)
167番目くらいのサンプルに平均発現量が以上に低いものがありますね。こういう外れ値だと考えられるサンプルをこのあと除去していきます！
一応各サンプルの発現量のボックスプロットも確認してみます
pdf(file = "../Output/boxplot_rawdata.pdf")
for (i in 1:(length(colnames(data$E)) %/% 50 + 1)){
print(i)
start = 50 * (i - 1) + 1
end = 50 * i
if (end > length(colnames(data$E))){
end = length(colnames(data$E))
}
boxplot(log2(data$E[,start:end]))
}
dev.off()
以下のような箱ひげ図を５０サンプルずつ描画してPDFに保存しています
こちらでもいくつかのサンプルで他とは異なるボックスプロットが見られるものがありますね。これはあとで除去します。
発現データの正規化
発現データの正規化についてはlimma User Guide 17.3 Comparing Mammary Progenitor Cell Populations with Illumina BeadChipsを参考にしています。
datExpr0 <- neqc(data)
めっちゃ簡単ですね。正規化はこれだけです。
limmaのneqc関数はNormexp法でバックグラウンド補正をしてquantile正規化をしてくれる関数です
正規化後のデータをboxplotで見てみましょう！
boxplot(datExpr0$E[,1:5])
ここまでの工程を再現するには時間がかかるのでこの段階で正規化後のデータを保存しておきます。
save(datExpr0, file = "../Output/001_NormalizeRawData/normalized_data.Rdata")
データをロードして再開するには以下のコマンドを実行します
load("../Output/001_NormalizeRawData/normalized_data.Rdata")
サンプルと遺伝子のフィルタリング
ここから発現していない遺伝子や、外れ値だと考えられるサンプルを除去していきます。まず遺伝子数とサンプル数が全部でいくつかを確認しておきます。
dim(datExpr0$E)
# [1] 47323 336
n_sample = dim(datExpr0$E)[2]
n_sample
# [1] 336
non-expressed geneのフィルタリング
ここでサンプル数全体の1/4以上で発現している（Detection p-valueが0.05以下）遺伝子という条件でフィルタリングを行います。
以下のコードはlimma User Guide 17.3.5 Normalization and filteringを参考にしています。limma User Guideでは12サンプル中３サンプル以上発現している遺伝子でフィルタリングを行っていたので今回も全体の25%以上という条件にしてみました。
expressed <- rowSums(datExpr0$other$Detection < 0.05) >= (n_sample/4)
datExpr<- datExpr0[expressed,]
dim(datExpr$E)
# [1] 17079 336
47323 → 17079と30000遺伝子ほどフィルタリングできました。
サンプルのフィルタリング
非発現遺伝子のフィルタリングが済んだので次はサンプルの外れ値を検出して除去していきます。これにはWGCNAチュートリアルでも採用されている階層的クラスタリングによる方法を用います。
Tutorial for the WGCNA package for R: I. Network analysis of liver expression data in female mice 1. Data input and cleaningを参照しています。
sampleTree = hclust(dist(t(datExpr$E)), method = "average")
save(sampleTree, file = "../Output/001_NormalizeRawData/sampleTree.Rdata")
作成したデンドログラムを描画していきます。以下のコードは階層的クラスタリングでどのサンプルがどのクラスターに分類されているのかを確認するためにSLEDAIごとの色付けなどを行っています。
dendextend, ggdendro パッケージについては以下を参考にしました。
https://cran.r-project.org/web/packages/dendextend/vignettes/dendextend.html
まずデンドログラム描画用のデータを作成します。
library(dendextend)
library(ggdendro)
dend <- as.dendrogram(sampleTree)
dend_data <- dendro_data(dend, type = "rectangle")
names(dend_data)
# [1] "segments" "labels" "leaf_labels" "class"
dend_dataの中身は枝の座標と葉のラベルなどです。
次にSLEDAIごとにカラーバーを表示するためのデータを作成していきます。ここは特に重要というわけではないので読み飛ばして頂いても大丈夫です。
sledai_group_tmp <- tibble(array_ID=labels(dend)) %>%
left_join(Traits, by=c("array_ID"="description") )
sledai_group = sledai_group_tmp$sledai_group
names(sledai_group) = sledai_group_tmp$array_ID
sledai_color_palette = c("#E34A33", "#FDBB84", "grey")
sledai_color = sledai_color_palette[unclass(as.factor(sledai_group))]
そしてcutHeight=140としてクラスターを作成します。この基準はあとで描画するデンドログラムを見て決めたものです。
cutHeight = 140
clust = cutreeStatic(sampleTree, cutHeight = cutHeight, minSize = 1)
names(clust) = sampleTree$labels
clust = clust[labels(dend)]
table(clust)
# clust
# 1 2 3 4 5 6 7 8 9 10
# 325 2 2 1 1 1 1 1 1 1
ほとんどのサンプル（325サンプル）がクラスター1に属していることがわかります。
各クラスターで色分けしてデンドログラムを描画するためのデータを作成します。
cluster_color_palette = brewer.pal(length(table(clust)), "Set3")
cluster_color = cluster_color_palette[clust]
names(cluster_color) = names(clust)
さて、やっとデンドログラムを描画できます。
dend <- dend %>%
set("labels_col", cluster_color) %>% # change color
set("labels_cex", 0) %>%
set("by_labels_branches_col", names(clust[clust == 1]),TF_values = cluster_color_palette[1])
plot(dend, ylim = c(-250, 320))
abline(h = cutHeight, col = "red")
colored_bars(dend = dend,
colors = data.frame(sledai=sledai_color, cluster=cluster_color),
rowLabels = c("sledai", "cluster"),
sort_by_labels_order = FALSE,
y_shift=-300, y_scale = 100)
結果のデンドログラムは上のようになります。左の方に外れ値となっているサンプルが集まっているのが確認できます。
デンドログラムの下のカラーバーは各サンプルのSLEDAIに基づいて色づけしています。controlがグレー、Low SLEDAIがオレンジ、High SLEDAIが赤です。この３グループで特にどのグループに外れ値が偏っているということもなさそうですのでこのまま解析を進めても問題なさそうです。
325サンプルのクラスター1が今後解析対象とするサンプルで、他のクラスターに属するものを外れ値として除去します。
keepSamples = (clust[colnames(datExpr)]==1)
notKeptSamples = colnames(datExpr)[!keepSamples]
datExpr = datExpr[,keepSamples]
nGenes = ncol(datExpr$E)
nSamples = nrow(datExpr$E)
GSE138458_series_matrix.txt中の記載によると以下の６サンプルを外れ値として解析対象から除外したそうです
“6 samples were determined to be outliers and were dropped prior to normalization: 200319680117_E 200308380073_J 200308380073_K 200308380073_L 200308380104_E 200667730029_E”
今回の基準で外れ値とされたサンプルを確認してみましょう。
notKeptSamples
# [1] "200319680117_E" "200319680117_G" "200363680045_G" "200319680125_C" "200319680125_K" "200308380073_J"
# [7] "200308380073_K" "200308380073_L" "200308380104_E" "200308380104_I" "200667730029_E"
論文の筆者らによって外れ値とされたものはすべて含まれていることが確認できます。
このことからある程度信頼できる基準で外れ値を判定できたと考えてもいいと思われます。
336→325と11サンプルを外れ値として除去しました。
遺伝子のフィルタリング
サンプルの外れ値を除去できたので、次は遺伝子のフィルタリングを行います。
今回参考にしたZhu et al. (2019)ではMean Absolute Deviation(MAD)で上位3000遺伝子を解析対象として採用しています。今回はその基準を参考に同様のフィルタリングを行っていきます。
MADは日本語だと平均絶対偏差と訳されます。つまり発現変動の大きい遺伝子が重要だと仮定してそれを解析対象としようという考えです。
WGCNAのFAQ 2.Should I filter probesets or genes?でも言及されているようにMADによる遺伝子のフィルタリングはWGCNAのための前処理として頻繁に使用されます。
gene_mad = apply(datExpr$E, MARGIN = 1, FUN = mad)
hist(gene_mad)
MADのヒストグラムを確認してみるとほとんどの遺伝子が0.5前後であることが確認できます。
ここからMADの値が上位3000位以内の遺伝子のみを抽出してみます。
gene_mad_rank = rank(-gene_mad, ties.method = "first")
keep = gene_mad_rank < 3000
hist(gene_mad[keep])
MADで1弱の値を持つ遺伝子を中心に取得することができました。もう少し基準を厳しく取得するのもありかもしれませんがまずはこの比較的緩めの基準で解析を進めることにします。WGCNAによるネットワーク構築がうまく行かなかった際にはまずここの基準を再考するべきでしょう。
datExpr <- datExpr[keep,]
dim(datExpr)
最後にgoodSamplesGenes関数で欠損値があまりに多い遺伝子やサンプルが無いかを確認してやります。
gsg = goodSamplesGenes(datExpr$E,verbose = 3)
gsg$allOK
# [1] TRUE
gsg$allOKがTRUEなので問題なさそうです！
これで前処理はすべて終了です！
あとは前処理したデータを保存しておきます。
解析対象となったサンプル以外のものをTraitsから除去しておきます。
datTraits0 <- Traits %>%
column_to_rownames("description")
datTraits <- datTraits0[rownames(datExpr),]
次の章で使用するWGCNAの関数群は行にサンプル、列に遺伝子のマトリクスを引数にするためdatExprを転置して保存します。
datExpr <- t(datExpr$E)
datTraits0 <- Traits %>%
column_to_rownames("description")
datTraits <- datTraits0[rownames(datExpr),]
save(datExpr, datTraits, file = "../Output/001_NormalizeRawData/SLE_GSE138458-01-dataInput.RData")
ここまで読んでいただいてありがとうございました。と言ってもWGCNA自体はまだ全然やってないんですが(笑)
次の記事ではいよいよWGCNAを使って遺伝子ネットワークを構築していきます！
前回のPart1ではrawデータのダウンロードから正規化を行い、サンプルと遺伝子フィルタリングまでを行いました。 [sitecard subtitle=関連記事 url= https://shiokoji11235.com[…]
参考
- Zhu et al. (2019)
- WGCNA tutorial
- WGCNA FAQ
- limma User Guide
- dendextend: https://cran.r-project.org/web/packages/dendextend/vignettes/dendextend.html
- ggdendro: http://www.sthda.com/english/wiki/beautiful-dendrogram-visualizations-in-r-5-must-known-methods-unsupervised-machine-learning