寫在前面
之前我們完成了WGCNA輸入數(shù)據(jù)的清洗虱肄，網(wǎng)絡構(gòu)建和模塊識別致板。??
但這里我們的示例數(shù)據(jù)內(nèi)所含有的基因其實是很少的，而在實際情況中咏窿，一個簡單的測序可能就要包含上萬個基因斟或，這對大家的電腦無疑是不小的壓力。??
在WGCNA的包內(nèi)其實也提供了解決方案集嵌，基本思想是分級聚類萝挤。??
1?? 首先御毅，我們使用快速但相對粗糙的聚類方法，用于將基因預聚類成大小接近的模塊怜珍，且不超過你所設定的基因最大值端蛆。??
2?? 然后我們分別在每個模塊中執(zhí)行完整的網(wǎng)絡分析。??
3?? 最后酥泛，合并特征基因高度相關的模塊今豆。??
用到的包
rm(list = ls())
library(WGCNA)
library(tidyverse)
示例數(shù)據(jù)
load("FemaleLiver-01-dataInput.RData")
軟閾值
4.1 topology analysis
首先我們還是要和之前一樣進行soft thresholding power β的計算。??
powers <- c(c(1:10), seq(from = 12, to=20, by=2))
sft <- pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
4.2 可視化
顯然揭璃，我們的結(jié)果和之前是一樣的晚凿，6亭罪。??
sizeGrWindow(9, 5)
par(mfrow = c(1,2))
cex1 = 0.9
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,col="red")
abline(h=0.90,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,col="red")
構(gòu)建網(wǎng)絡與模塊識別
5.1 網(wǎng)絡構(gòu)建
這里我們就要設置每個block的最大size是多少了瘦馍，分次計算，再合并应役，以此減少內(nèi)存的負擔情组。??
bwnet <- blockwiseModules(
datExpr, maxBlockSize = 2000,
power = 6, TOMType = "unsigned", minModuleSize = 30,
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = T,
saveTOMs = T,
saveTOMFileBase = "femaleMouseTOM-blockwise",
verbose = 3)
5.2 查看模塊數(shù)
table(bwnet$colors)
對比結(jié)果
這里我們再對比一下結(jié)果，看看2種方法得出結(jié)果的區(qū)別箩祥。??
6.1 載入之前的結(jié)果
這里我們把之前不分次計算的結(jié)果載入進來院崇，后面會用到label和colors。??
load(file = "FemaleLiver-02-networkConstruction-auto.RData")
6.2 匹配顏色
把colors和label匹配起來袍祖，嘿嘿底瓣。??
bwLabels <- matchLabels(bwnet$colors, moduleLabels);
bwModuleColors <- labels2colors(bwLabels)
6.3 分次結(jié)果可視化
sizeGrWindow(6,6)
# Block 1
plotDendroAndColors(bwnet$dendrograms[[1]], bwModuleColors[bwnet$blockGenes[[1]]],
"Module colors", main = "Gene dendrogram and module colors in block 1",
dendroLabels = F, hang = 0.03,
addGuide = T, guideHang = 0.05)
# Block 2
plotDendroAndColors(bwnet$dendrograms[[2]], bwModuleColors[bwnet$blockGenes[[2]]],
"Module colors", main = "Gene dendrogram and module colors in block 2",
dendroLabels = F, hang = 0.03,
addGuide = T, guideHang = 0.05)
對比兩種方法的結(jié)果差異
7.1 對比一下
我們以可視化的形式對比一下，分割出來的模塊差異不大蕉陋。??
sizeGrWindow(12,9)
plotDendroAndColors(geneTree,
cbind(moduleColors, bwModuleColors),
c("Single block", "2 blocks"),
main = "Single block gene dendrogram and module colors",
dendroLabels = F, hang = 0.03,
addGuide = T, guideHang = 0.05)
7.2 對比eigengenes
這里我們提取一下2種方法得到的module eigengenes捐凭。??
singleBlockMEs <- moduleEigengenes(datExpr, moduleColors)$eigengenes
blockwiseMEs <- moduleEigengenes(datExpr, bwModuleColors)$eigengenes
match之后看一下結(jié)果，嘿嘿凳鬓。??
高度一致茁肠，所以這種blockwise的方法，請放心食用吧缩举，各位垦梆。??
single2blockwise <- match(names(singleBlockMEs), names(blockwiseMEs))
signif(diag(cor(blockwiseMEs[, single2blockwise], singleBlockMEs)), 3)
如何引用
??
Langfelder, P., Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559 (2008). https://doi.org/10.1186/1471-2105-9-559
<center>最后祝大家早日不卷!~</center>
點個在看吧各位~ ?.???? ??? ?
<center> <b>?? 往期精彩 <b> </center>
?? <font size=1>?? ComplexHeatmap | 顏狗寫的高顏值熱圖代碼！</font>
?? <font size=1>?? ComplexHeatmap | 你的熱圖注釋還擠在一起看不清嗎＝龊ⅰ托猩？</font>
?? <font size=1>?? Google | 谷歌翻譯崩了我們怎么辦！辽慕？（附完美解決方案）</font>
?? <font size=1>?? scRNA-seq | 吐血整理的單細胞入門教程</font>
?? <font size=1>?? NetworkD3 | 讓我們一起畫個動態(tài)的烧拘蹋基圖吧~</font>
?? <font size=1>?? RColorBrewer | 再多的配色也能輕松搞定！~</font>
?? <font size=1>?? rms | 批量完成你的線性回歸</font>
?? <font size=1>?? CMplot | 完美復刻Nature上的曼哈頓圖</font>
?? <font size=1>?? Network | 高顏值動態(tài)網(wǎng)絡可視化工具</font>
?? <font size=1>?? boxjitter | 完美復刻Nature上的高顏值統(tǒng)計圖</font>
?? <font size=1>?? linkET | 完美解決ggcor安裝失敗方案（附教程）</font>
?? <font size=1>......</font>
本文由mdnice多平臺發(fā)布
