Part1写在前面

WGCNA的教程，我们在前期的推文中已经推出好久了。今天再结合前期的教程进行优化。只是在现有的教程基础上，进行修改。其他的并无改变。

Part2前期WGCNA教程

• WGCNA分析 | 全流程分析代码 | 代码一

• WGCNA分析 | 全流程分析代码 | 代码二

• WGCNA分析 | 全流程代码分享 | 代码三

Part3本次教程的优化点

注意：本次教程在教程二的基础上修改。

1. 教程代码更规范
2. 教程添加了过滤数值流程
3. 教程添加merge模块后图形绘制（注：在教程二的基础上）
4. 教程添加更多的注释信息
5. 教程后期添加视频讲解
6. 若已经购买教程二的同学，可以直接向小杜索要本教程代码（注意：请在教程二中留言。若购买此教程，留言可获得本教程全部代码）

Part4如何加入社群

小杜的生信笔记，目前有两个社群（QQ群和微信群）。1. QQ群：刚刚建立，免费加入社群。主要是社群同学自己探讨问题，自己不会过多参与和涉及，会不定期在此群发送相关学习资料。加入方式：2. 微信群：付费社群。添加小杜好友，加友请知：加友须知！！

3. 欢迎投稿

小杜一直在分享自己平时学习笔记，因此分享的内容大多数是与自己相关的，局限性比较大。我一直在倡导大家一起分享自己学习笔记或教程。分享内容不限于生信教程，可以是文章or文献、遇到的问题及解决方案、学习感悟等等。

Part5WGCNA教程 | 代码四

本期教程输出所有的文件信息。

1设置文件位置及导入包

##'@加权基因共表达分析(WGCNA)
##'@2023.09.02
##'@整理者：小杜的生信笔记
##'@前期教程网址：https://mp.weixin.qq.com/s/Ln9TP74nzWhtvt7obaMp1A
##'
##'@本教程主要主要是为了优化前期代码，在前期的基础上进行修改。
##'@若您的数据量较大，我们建议WGCNA在服务器上跑。
##'

##==============================================================================

setwd("E:\\小杜的生信筆記\\2023\\20230217_WGCNA\\WGCNA_04")
rm(list = ls())
#install.packages("WGCNA")
#BiocManager::install('WGCNA')
library(WGCNA)
options(stringsAsFactors = FALSE)
#enableWGCNAThreads(60) ## 打开多线程
#Read in the female liver data set

2导入数据

我们这里给出了两种不同文件的导入方式，txt和csv。

#'@导入数据
#'@导入txt格式数据
# WGCNA.fpkm = read.table("ExpData_WGCNA.txt",header=T,
#                         comment.char = "",
#                         check.names=F)
##'@导入csv文件格式
WGCNA.fpkm = read.csv("ExpData_WGCNA.csv", header = T, check.names = F)
WGCNA.fpkm[1:10,1:10]

3检查数据缺失值

gsg = goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK
if (!gsg$allOK)
{
  if (sum(!gsg$goodGenes)>0)
    printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")))
  if (sum(!gsg$goodSamples)>0)
    printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")))
  # Remove the offending genes and samples from the data:
  datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}

4过滤数值 [optional]

@过滤数值（optional）,此步根据你自己的数据进行过滤数值，过滤的数值大小根据自己需求修改

meanFPKM=0.5  ###--过滤标准，可以修改
n=nrow(datExpr0)
datExpr0[n+1,]=apply(datExpr0[c(1:nrow(datExpr0)),],2,mean)
datExpr0=datExpr0[1:n,datExpr0[n+1,] > meanFPKM]
# for meanFpkm in row n+1 and it must be above what you set--select meanFpkm>opt$meanFpkm(by rp)
filtered_fpkm=t(datExpr0)
filtered_fpkm=data.frame(rownames(filtered_fpkm),filtered_fpkm)
names(filtered_fpkm)[1]="sample"
head(filtered_fpkm)
###'@输出过滤的文件
write.table(filtered_fpkm, file="WGCNA.filter.txt",
            row.names=F, col.names=T,quote=FALSE,sep="\t")

5样本聚类

###'@样本聚类
sampleTree = hclust(dist(datExpr0), method = "average")
pdf("1_sample clutering.pdf", width = 6, height = 4)
par(cex = 0.7);
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)
dev.off()

样本sample05与整体数据差异较大，过滤此数据。

6过滤样本

pdf("2_sample clutering_delete_outliers.pdf", width = 6, height = 4)
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, 
     cex.axis = 1.5, cex.main = 2) +
  abline(h = 1500, col = "red")    ###'@"h = 1500"参数为你需要过滤的参数高度
dev.off()

##'@过滤离散样本
##'@"cutHeight"为过滤参数，与上述图保持一致
clust = cutreeStatic(sampleTree, cutHeight = 1500, minSize = 10)
keepSamples = (clust==1)
datExpr = datExpr0[keepSamples, ]
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)

两个样本直接数据比较

7载入性状数据

##'@导入csv格式
traitData = read.csv("TraitData.csv",row.names=1)
# ##'@导入txt格式
# traitData = read.table("TraitData.txt",row.names=1,header=T,comment.char = "",check.names=F)
head(traitData)
allTraits = traitData
dim(allTraits)
names(allTraits)
# 形成一个类似于表达数据的数据框架
fpkmSamples = rownames(datExpr)
traitSamples =rownames(allTraits)
traitRows = match(fpkmSamples, traitSamples)
datTraits = allTraits[traitRows,]
rownames(datTraits)
collectGarbage()

8Re-cluster samples

# Re-cluster samples
sampleTree2 = hclust(dist(datExpr), method = "average")
# 
traitColors = numbers2colors(datTraits, signed = FALSE)
# Plot the sample dendrogram and the colors underneath.

9Sample dendrogram and trait heatmap

#sizeGrWindow(12,12)
pdf(file="3_Sample_dendrogram_and_trait_heatmap.pdf",width=8 ,height= 6)
plotDendroAndColors(sampleTree2, traitColors,
                    groupLabels = names(datTraits),
                    main = "Sample dendrogram and trait heatmap",cex.colorLabels = 1.5, cex.dendroLabels = 1, cex.rowText = 2)
dev.off()

数据处理结束，继续后续网络构建分析！

#'@打开多线程分析
enableWGCNAThreads(5)

10获得最佳阈值（softpower）

#'@获得soft阈值
#powers = c(1:30)
powers = c(c(1:10), seq(from = 12, to=30, by=2))
#'@调用网络拓扑分析功能
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

11绘制softpower图形

#'@绘制soft power plot 
pdf(file="4_软阈值选择.pdf",width=12, height = 8)
par(mfrow = c(1,2))
cex1 = 0.85
# Scale-free topology fit index as a function of the soft-thresholding power
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");
# this line corresponds to using an R^2 cut-off of h
abline(h=0.8,col="red")
# Mean connectivity as a function of the soft-thresholding power
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")
dev.off()

12选择最佳的softpower

softpower是WGCNA分析影响最大的参数，因此，选择最佳的softpower参数直接影响你后续的结果。但是如何选择呢？？

1. 系统默认最佳阈值

softPower =sft$powerEstimate

但是，系统默认也不一定可以满足我们的需求，以及返回结果可能是"NA"，这是正常的表现，或是返回结果用于后面的分析，获得模型结果不符合自己的预期，因此，softpower值可以自己进行更改。

2. 自己手动更改

我们根据上图，自己手动选择相对较好的softpower即可，但是，也是需要不断的尝试。

13构建网络模块

1. 方法一：使用blockwiseModules()构建
2. 方法二：使用WGCNA教程一的方法进行构建

##'@直接使用blockwiseModules()函数构建网络模块
## 模块可视化
net = blockwiseModules(datExpr, power = 20 ,#手动改softpower
                       #signed, unsigned
                       TOMType = "signed",
                       ##'@模块数量，可更改，一般为：10，20，30
                       minModuleSize = 30,
                       reassignThreshold = 0 
                       #'@Module merging options
                       mergeCutHeight = 0.4,
                       corType = "pearson",   ## 相关性系数
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,maxBlockSize = 20000,
                       saveTOMFileBase = "MyTOM",
                       verbose = 3,    ##'@冗长程度的整数。零表示无声，数值越大，输出就越冗长
                       nThreads = 4     ##'@计算时使用线程数
                       )
table(net$colors)

blockwiseModules()函数更多参数，可以自己进行查看。

??blockwiseModules()

注意：

1. 使用blockwiseModules（）函数可以直接合并merge模块的参数

mergeCutHeight: 及时合并模块的参数，默认为0.25

Dynamic Tree

mergedColors = labels2colors(net$colors)
table(mergedColors)
pdf(file="5_Dynamic Tree Cut.pdf",width= 8,height=6)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
dev.off()

生信知识库网址，前期订阅无需重复订阅：