写在前面
最近遇到了非常不讲理的病人和家属，真是忍不住想要吐槽两句。😤
为了让病人做上手术，求了输血科半天才给备了血，家属也答应去献血。😒
万万没想到，术后都2天了还是没去献血，无论是问病人还是问家属，都说家里没有一个人适合献血的。🤒
用别人的血就天经地义、理所当然，让你去献血就是大夫在害你了，我也是没得话讲，农夫与蛇的现代翻版。😷
言归正传了，继续之前的WGCNA教程吧，希望坚持给大家分享知识也能得到大家的理解。😞
用到的包
rm(list = ls())
library(tidyverse)
library(WGCNA)
示例数据
load("./Consensus-dataInput.RData")
提取数据集数量
提取一下我们的数据集个数，后面会用到。🤓
nSets <- checkSets(multiExpr)$nSets
构建网络与模块识别（一步法）
5.1 创建power向量
powers <- c(seq(4,10,by=1), seq(12,20, by=2))
5.2 计算power
接着我们为每个数据集计算一下power，挑选软阈值。😋
powerTables <- vector(mode = "list", length = nSets)
for (set in 1:nSets){
powerTables[[set]] = list(data = pickSoftThreshold(multiExpr[[set]]$data,
powerVector=powers,
verbose = 2)[[2]])
collectGarbage()
}
5.3 可视化参数设置
我们可视化一下结果吧，先设置一下参数。🤓
colors = c("black", "red")
plotCols = c(2,5,6,7)
colNames = c("Scale Free Topology Model Fit", "Mean connectivity", "Median connectivity", "Max connectivity")
ylim = matrix(NA, nrow = 2, ncol = 4);
for (set in 1:nSets)
{for (col in 1:length(plotCols))
{
ylim[1, col] = min(ylim[1, col], powerTables[[set]]$data[, plotCols[col]], na.rm = T);
ylim[2, col] = max(ylim[2, col], powerTables[[set]]$data[, plotCols[col]], na.rm = T);
}
}
5.4 可视化
ok，可以正式可视化了。😗
sizeGrWindow(8, 6)
par(mfcol = c(2,2))
par(mar = c(4.2, 4.2 , 2.2, 0.5))
cex1 = 0.7
for (col in 1:length(plotCols)) for (set in 1:nSets)
{
if (set==1)
{
plot(powerTables[[set]]$data[,1], -sign(powerTables[[set]]$data[,3])*powerTables[[set]]$data[,2],
xlab="Soft Threshold (power)",ylab=colNames[col],type="n", ylim = ylim[, col],
main = colNames[col]);
addGrid();
}
if (col==1)
{
text(powerTables[[set]]$data[,1], -sign(powerTables[[set]]$data[,3])*powerTables[[set]]$data[,2],
labels=powers,cex=cex1,col=colors[set]);
} else
text(powerTables[[set]]$data[,1], powerTables[[set]]$data[,plotCols[col]],
labels=powers,cex=cex1,col=colors[set]);
if (col==1)
{
legend("bottomright", legend = setLabels, col = colors, pch = 20) ;
} else
legend("topright", legend = setLabels, col = colors, pch = 20) ;
}
5.5 识别模块
net <- blockwiseConsensusModules(multiExpr, power = 6, minModuleSize = 30,
deepSplit = 2, pamRespectsDendro = F,
mergeCutHeight = 0.25,
# maxBlockSize = 20000,
numericLabels = T, minKMEtoStay = 0,
saveTOMs = T,
verbose = 5)
5.6 查看网络
names(net)
5.7 转换label
consMEs <- net$multiMEs
moduleLabels <- net$colors
moduleColors <- labels2colors(moduleLabels)
consTree <- net$dendrograms[[1]]
可视化网络及模块结果
sizeGrWindow(8,6)
plotDendroAndColors(consTree, moduleColors, "Module colors",
dendroLabels = F, hang = 0.03,
addGuide = T, guideHang = 0.05,
main = "Consensus gene dendrogram and module colors")
Save一下
保存一下data吧，毕竟运行了那么久。🥳
save(consMEs, moduleLabels, moduleColors, consTree,
file = "./Consensus-NetworkConstruction-auto.RData")
<center>最后祝大家早日不卷!~</center>
点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰
<center> <b>📍 往期精彩 <b> </center>
📍 <font size=1>🤩 WGCNA | 值得你深入学习的生信分析方法！~</font>
📍 <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 | 让我们一起画个动态的桑基图吧~</font>
📍 <font size=1>🤩 RColorBrewer | 再多的配色也能轻松搞定！~</font>
📍 <font size=1>🧐 rms | 批量完成你的线性回归</font>
📍 <font size=1>🤩 CMplot | 完美复刻Nature上的曼哈顿图</font>
📍 <font size=1>🤠 Network | 高颜值动态网络可视化工具</font>
📍 <font size=1>🤗 boxjitter | 完美复刻Nature上的高颜值统计图</font>
📍 <font size=1>🤫 linkET | 完美解决ggcor安装失败方案（附教程）</font>
📍 <font size=1>......</font>
本文由mdnice多平台发布
