

生信小课堂(4) R中轻松进行WGCNA分析

R语言数据分析指南

2023-10-12

欢迎关注R语言数据分析指南

❝
本节分享如何基于表达矩阵的结果来做WGCNA分析以及模块结果的导出，过程仅供参考后期可以根据导出的结果对图表进行精细的加工，希望对各位观众老爷能有帮助，数据和代码已经被打包并上传到小编的2023VIP会员交流群。已经加群的朋友可以自行下载。如果你需要，可以参考文末的方式进行购买。

实际运行中要根据自己的实际数据来调整部分参数，运行所需内存至少需要30G+
❞

关注下方公众号下回更新不迷路

加载R包

library(WGCNA)
library(tidyverse)
library(cowplot)

导入数据

gene_exp <- read.table(file ='heatmap.txt',row.names=1)
datExpr0 <- t(gene_exp) # 转置

数据过滤

# 缺失数据及无波动数据过滤
gsg <- goodSamplesGenes(datExpr0,minFraction = 1/2) #基因的缺失数据比例阈值
datExpr <- datExpr0[gsg$goodSamples, gsg$goodGenes]

# 通过聚类，查看是否有明显异常样本, 如果有需要剔除掉
plot(hclust(dist(datExpr)),cex.lab = 1.5,cex.axis = 1.5,cex.main = 2)

enableWGCNAThreads(nThreads =40) # 设置线程数
# 通过对 power 的多次迭代，确定最佳 power
sft <- pickSoftThreshold(datExpr, powerVector = 1:20, # 尝试 1 到 20
                         networkType="unsigned")

fig_power1 <- ggplot(data = sft$fitIndices,aes(x = Power,y = SFT.R.sq)) +
  geom_point(color = 'red') +
  geom_text_repel(aes(label = Power)) +
  geom_hline(aes(yintercept = 0.85), color = 'red') +
  labs(title = 'Scale independence',
       x = 'Soft Threshold (power)',
       y = 'Scale Free Topology Model Fit,signed R^2') +
  theme_few() +
  theme(plot.title = element_text(hjust = 0.5))

fig_power2 <- ggplot(data = sft$fitIndices,
                     aes(x = Power,
                         y = mean.k.)) +
  geom_point(color = 'red') +
  geom_text_repel(aes(label = Power)) +
  labs(title = 'Mean connectivity',
       x = 'Soft Threshold (power)',
       y = 'Mean Connectivity') +
  theme_few()+
  theme(plot.title = element_text(hjust = 0.5))

plot_grid(fig_power1, fig_power2)

sft$powerEstimate  #  最佳B值

net <- blockwiseModules(datExpr,  # 输入数据
  corType = "pearson", # 相关系数算法，pearson|bicor
  power = 5, # 前面得到的 soft power
  networkType = "unsigned", # unsigned | signed | signed hybrid
  # 计算 TOM 矩阵
  TOMType = "unsigned", # none | unsigned | signed
  saveTOMs = TRUE,
  saveTOMFileBase = "blockwiseTOM",
  # 聚类并划分模块
 # deepSplit = 0, # 0|1|2|3|4, 值越大得到的模块就越多越小，不设置则使用默认值
  minModuleSize =30,
  mergeCutHeight = 0.25,  #  对距离小于mergeCutHeight的模块进行合并
  numericLabels = FALSE, # 以数字命名模块
  nThreads = 0, # 0 则使用所有可用线程
  maxBlockSize = 100000 # 需要大于基因的数目
)

table(net$colors) # 查看每个模块包含基因数目
save(net,file="net.Rdata") # 结果保存

         1416          4980          4271           139           118           107 
   darkorange       darkred darkturquoise         green   greenyellow          grey 
          100           124           112          2119           151          1014 
       grey60     lightcyan    lightgreen   lightyellow       magenta  midnightblue 
          137           138           136           133           175           138 
       orange paleturquoise          pink        purple           red     royalblue 
          101            58          1249           166          1720           127 
  saddlebrown        salmon       skyblue     steelblue           tan     turquoise 
           80           143            89            69           148          5221 
       violet         white        yellow 
           36            96          3898

WGCNA的结果导出

wgcna_result <- data.frame(gene_id = names(net$colors),module = net$colors)

wgcna_result %>% write.table(file="wgcna_result.xls",row.names = F,sep="\t",quote = F)

绘制聚类图与模块

plotDendroAndColors(
    dendro = net$dendrograms[[1]], 
    colors = net$colors,
    groupLabels = "Module colors",
    dendroLabels = FALSE, 
    addGuide = TRUE)

rdatTraits <- read.table("group.xls",sep="\t",header = T,row.names = 1)
design=model.matrix(~0+ datTraits$group)
colnames(design)=levels(factor(datTraits$group))
rownames(design) <- rownames(datTraits)

design %<>% as.data.frame()
moduleTraitCor <- cor(net$MEs, design,use = "p",method = 'pearson')
  
  # 计算 Pvalue
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor,nrow(datExpr))

sizeGrWindow(10,6)
# 连接相关性和 pvalue
textMatrix <- paste(signif(moduleTraitCor, 2), "\n(",
                      signif(moduleTraitPvalue, 1), ")", sep = "");
  dim(textMatrix) <- dim(moduleTraitCor)
  par(mar = c(6,10,3,3))
  labeledHeatmap(Matrix = moduleTraitCor,
  xLabels = names(design),yLabels = names(net$MEs),
  ySymbols = names(net$MEs),colorLabels = FALSE,
  colors = blueWhiteRed(50),textMatrix = textMatrix,
  setStdMargins = FALSE,cex.text = 0.5,zlim = c(-1,1),
  main = paste("Module-trait relationships"))

moduleTraitCor %>% write.table(file="moduleTraitCor.xls",sep="\t",quote = F)

模块的导出

my_modules <- c('yellow')
# 提取该模块的表达矩阵
m_wgcna_result <- filter(wgcna_result, module %in% my_modules)
m_datExpr <- datExpr[, m_wgcna_result$gene_id]
# 计算该模块的 TOM 矩阵  
m_TOM <- TOMsimilarityFromExpr(m_datExpr,
    power = 20,networkType = "unsigned",
    TOMType = "unsigned")
dimnames(m_TOM) <- list(colnames(m_datExpr), colnames(m_datExpr))

# 导出 Cytoscape 输入文件
cyt <- exportNetworkToCytoscape(m_TOM,
    edgeFile = "CytoscapeInput-yellow-network.txt",
    weighted = TRUE,threshold = 0.2)

MEs <- net$MEs # 修改名称
colnames(MEs) <- str_remove(colnames(MEs),"ME")
save(MEs,design,file="df.Rdata") # 结果保存

❝
本节内容介绍到此结束，过程仅供参考；有需要学习个性化数据可视化的朋友，欢迎到小编的「淘宝店铺」 「R语言数据分析指南」购买「2023年度会员文档」同步更新中「售价149元」，内容主要包括各种「高分论文的图表分析复现以及一些个性化图表的绘制」均包含数据+代码；按照往年数据小编年产出约在150+以上
❞

购买后微信发小编订单截图即邀请进新的会员交流群，小编的文档为按年售卖，只包含当年度的「除系列课程外」的文档，有需要往年文档的朋友也可下单购买，需要了解更多信息的朋友欢迎淘宝店铺交流咨询。