データの読み込みとクラスタリング
library(WGCNA)
# ディレクトリ名の確認
dir()
# データの読み込み
load('DEG_TB_LTBI_step13.Rdata')
# 行がサンプル名、列が遺伝子名
################################ サンプルクラスタリング ####################
datExpr = t(dataset_TB_LTBI_DEG)
# 初期クラスタリング
sampleTree = hclust(dist(datExpr), method = "average")
# サンプルツリーのプロット: 12x9インチのグラフィックウィンドウを開く
sizeGrWindow(12,9)
par(cex=0.6)
par(mar=c(0,4,2,0))
plot(sampleTree, main = "外れ値検出のためのサンプルクラスタリング", sub="", xlab="", cex.lab = 1.5,
cex.axis = 1.5, cex.main = 2)
# 結果画像をPDFでエクスポート: ファイル名=1_sampleClustering.pdf
#pdf(file='sampleCluestering.pdf',width = 12,height = 9)
# カットラインのプロット
abline(h = 87, col = "red")##カット高さは不確定なので赤線なし
dev.off()
初期クラスタリングの観察
カットラインの位置を定義して分割
この例では、右側に孤立したサンプルが見られ、削除するのに適しているため、87のカットラインを設定します。
左側も2つのブロックに分割されるため、処理が必要です。保持します。
### ライン下のクラスタを決定
clust = cutreeStatic(sampleTree, cutHeight = 87, minSize = 10)
table(clust)
#clust
#0 1 2
#5 57 40
### clust 1が保持したいサンプルを含む
keepSamples = (clust==1|clust==2)
datExpr0 = datExpr[keepSamples, ]
dim(datExpr0) #[1] 97 2813
# データの保存
save(datExpr0,file='datExpr0_cluster_filter.Rdata')
性状データの読み込み
サンプル名をマッチさせて、性状データと発現データが一致することを保証します。
#################### 性状データの読み込み ###########################
# 自分の性状データをロード
load('design_TB_LTBI.Rdata')
traitData=design
# 臨床性状データの読み込み
#traitData = read.table("trait_D.txt",row.names=1,header=T,comment.char = "",check.names=F)########性状データファイル名は実際に合わせて変更してください。もし作業ディレクトリが性状データのパスでない場合は、正しいパスを追加してください。
dim(traitData)
#names(traitData)
# 臨床性状を保持する発現データに類似したデータフレームを作成します。
fpkmSamples = rownames(datExpr0)
traitSamples =rownames(traitData)
# サンプル名をマッチさせて、性状データと発現データが一致することを保証します。
traitRows = match(fpkmSamples, traitSamples)
datTraits = traitData[traitRows,]
rownames(datTraits)
collectGarbage()
形状情報を追加して再クラスタリング
# サンプルの再クラスタリング
sampleTree2 = hclust(dist(datExpr0), method = "average")
# 性状を色表現に変換: 白は低値、赤は高値、灰色は欠損値
traitColors = numbers2colors(datTraits, signed = FALSE)
# サンプルデンドログラムとその下の色をプロットします。
#sizeGrWindow(20,20)
##pdf(file="2_Sample dendrogram and trait heatmap.pdf",width=12,height=12)
plotDendroAndColors(sampleTree2, traitColors,
groupLabels = names(datTraits),
main = "サンプルデンドログラムと性状ヒートマップ")
dev.off()
下の赤色部分は、概ね2つのクラスターに分かれており、良好な結果です。
ネットワーク構築
############################# ネットワーク構築 ################################
# WGCNA内でのマルチスレッディングを許可します。現在のところ、この呼び出しは必要です。
# ここでのエラーは無視しても構いませんが、表示される場合はWGCNAを更新することをお勧めします。
# 注意: RStudioやその他のサードパーティR環境を使用する場合は、この行をスキップしてください。
# 上記の注意を参照してください。
enableWGCNAThreads()
# ソフトスレッショルディングパワーのセットを選択
powers = c(1:15)
# ネットワークトポロジー分析関数を呼び出します
sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)
# 結果をプロットします:
sizeGrWindow(15, 9)
#pdf(file="3_Scale independence.pdf",width=9,height=5)
par(mfrow = c(1,2))
cex1 = 0.9
# ソフトスレッショルディングパワーに対するスケールフリートポロジーフィット指数
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="ソフトスレッショルディングパワー",ylab="スケールフリートポロジーモデルフィット、符号付きR^2",type="n",
main = paste("スケール独立性"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1,col="red");
# このラインはR^2カットオフhに対応します
abline(h=0.90,col="red")
# ソフトスレッショルディングパワーに対する平均接続度
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="ソフトスレッショルディングパワー",ylab="平均接続度", type="n",
main = paste("平均接続度"))
tex