WGCNA分析_main_ppl

WGCNA分析的步骤通常包括以下几个关键步骤：

1. 数据预处理：首先，对原始的基因表达数据进行预处理，包括数据清洗、归一化和筛选差异表达基因等。

2. 构建共表达网络：基于预处理后的表达数据，计算基因之间的相关性或相关系数，构建一个加权的共表达网络。常用的相关性度量方法包括Pearson相关系数和Spearman相关系数。

3. 模块发现：通过聚类算法（如层次聚类或k-means聚类）将高度相关的基因分组为模块。常用的模块发现算法是动态树切割方法，它基于网络的拓扑结构和模块的内部连接性来确定模块的划分。

4. 模块注释和功能分析：对于每个模块，可以进行进一步的注释和功能分析，例如富集分析，以确定与该模块相关的生物学过程或功能通路。

5. 基因集对比：可以将不同条件下的样本进行比较，找出在特定条件下不同模块的变化，并通过基因集对比分析来研究关键基因和通路的变化。

6. 网络可视化：最后，可以使用可视化工具将共表达网络和模块可视化，以便更好地理解整个数据集的结构和模式。

WGCNA分析可以帮助揭示基因表达数据中的隐藏模式和生物学机制，有助于深入理解基因调控网络和相关生物学过程。它在基因组学研究和系统生物学中被广泛应用。

library(WGCNA)
library(reshape2)
library(stringr)
library(stringi)
library(dplyr)

######1.1 导入数据

setwd("./DK_drought")
dataExpr <- read.csv(file = '../3.DK-Drought-gene.csv', sep=',',row.names=1, header=T)#这里把sep=','的/改成逗号
#dim查看
dim(dataExpr)

#35775 69

dataExpr<-dataExpr[which(rowSums(dataExpr) > 0),]

#colnames函数作用 重命名“列”
colnames(dataExpr)

#mad中位绝对偏差
m.mad <- apply(dataExpr,1,mad)
dataExprVar <- dataExpr[which(m.mad > 
                                max(quantile(m.mad, probs=seq(0, 1, 0.25),)[2],0.01)),]
dataExpr <- as.data.frame(t(dataExprVar))


dataExpr<-dataExpr*1.0 #注意要有 *1.0


#转换为样品在行，基因在列的矩阵
dataExpr <- as.data.frame(t(dataExpr))

#####1.2 检测缺失值
gsg = goodSamplesGenes(dataExpr, verbose = 3)
if (!gsg$allOK)
{
  # 打印删除的基因和样本名称
  if (sum(!gsg$goodGenes)>0) 
    printFlush(paste("Removing genes:", 
                     paste(names(dataExpr)[!gsg$goodGenes], collapse = ",")));
  if (sum(!gsg$goodSamples)>0) 
    printFlush(paste("Removing samples:", 
                     paste(rownames(dataExpr)[!gsg$goodSamples], collapse = ",")));
  
  #从数据中删除有问题的基因和样本
  dataExpr = dataExpr[gsg$goodSamples, gsg$goodGenes]
}

#统计基因和样本
nGenes = ncol(dataExpr)
nSamples = nrow(dataExpr)

dim(dataExpr)#dim函数查看

#####1.3 接下来，对样本进行聚类（与随后的基因聚类相比），看看是否有明显的异常值。

###首先，查看是否有离群样本
# hclusts聚类算法, dist计算基因之间的距离
sampleTree = hclust(dist(dataExpr), method = "average")

# 出图---样本聚类。 main标题，sub:子标题，xlab:x轴标题

pdf("1_sampletree.pdf",height=20,width=30)
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="")
dev.off()

pdf("2_sampletree.pdf",height=20,width=30)
#1.4 去除聚类高度高于1500000的样本
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)+abline(h = 1500000, col = "red");
dev.off()

###############二、软阈值的选择：网络拓扑分析
#2.1步进法 有助于自定义各类参数
#构建一个加权基因网络需要选择软阈值幂β来计算邻接矩阵权重参数,即将基因间的相关系数进行乘方运算来表征其相关性，
#首先需要确定乘方的值
# 2.1.1 确定power（β）值。 c(1:10)表示1到10；seq(from = 12, to=20, by=2)表示从12开始，间隔两个数到20
powers = c(c(1:10), seq(from = 12, to=30, by=2))

# 调用网络拓扑分析函数
# verbose表示输出结果详细程度，获得各个阈值下的 R方 和 平均连接度
sft = pickSoftThreshold(dataExpr, powerVector=powers, 
                        networkType='unsigned', verbose=5)

# 设置图的显示一行两列
pdf("3_line.pdf")
par(mfrow = c(1,2))
cex1 = 0.9

## 生成阈值和网络的特征之间的关系函数
# x横轴是Soft threshold (power)，y纵轴是无标度网络的评估参数，数值越高，
#则网络越符合无标度特征 (non-scale)
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")

#设置筛选标准，R-square(即R^2)=多少，根据图来设置。
# 这里选择R-square=0.89，即h=0.89。
abline(h=0.98,col="red")


# sft$fitIndices 保存了每个power构建的相关性网络中的连接度的统计值，k就是连接度值，
#每个power值提供了max, median, max3种连接度的统计量
# 对连接度的均值进行可视化
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")

#左图：Soft Threshold (power)表示权重，纵坐标表示连接度k与p(k)的相关性。
#右图：Soft Threshold (power)表示权重，纵坐标表示平均连接度。
#一般要求k与p(k)的相关性达到0.85时的power作为β值，（由于我们这里选的相关性是0.89，
#所以在左图纵坐标相关性为0.89时，所对应的横坐标power是9。即β=9）

dev.off()


#####一步构建网络与模块检测，把所有的基因分为不同的基因模块
#确定好power值即可构建基因网络，构建基因网络和识别模块，现在是一个简单的函数调用：

#将上面的dataExpr赋值于dataExpr1
dataExpr1<-apply(dataExpr,2,as.numeric)


# datExpr表达数据，TOMType拓扑重叠矩阵计算方式，minModuleSize用于模块检测的最小模块尺寸,
# reassignThreshold 是否在模块之间重新分配基因的p值比率阈值，mergeCutHeight 树状图切割高度
# numericLabels 返回的模块应该用颜色（FALSE）还是数字（TRUE）标记,pamRespectsDendro树状图相关参数
# saveTOMs 字符串的向量，saveTOMFileBase 包含包含共识拓扑重叠文件的文件名库的字符串
#此处选择的软阈值为9；设置模块中包含的基因个数最小为30（因为我们喜欢大的模块，因此这个值应该设置的尽可能大一些），
#一个中等的敏感度（deepSplit=2）
#参数mergeCutHeight是模块融合的阈值。这个函数返回的是数值，不是模块的颜色标签。
net = blockwiseModules(dataExpr1, power = 9, maxBlockSize = nGenes,
                       TOMType ='unsigned', minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs=TRUE, saveTOMFileBase = "drought",
                       verbose = 3)

#class(dataExpr)
#str(dataExpr)

#我们现在回到网络分析。要查看标识了多少个模块以及模块大小

#使用table(net＄colors)可以显示网络总共有多少个模块（分组），每个模块（组）中包含有多少个基因。
#根据模块中基因数目的多少，降序排列，依次编号为“1~最大模块数”
table(net$colors)

#层级聚类树展示各个模块
##**0（grey）**  表示 **没有** 分到任何模块的基因# 灰色的为**未分类**到模块的基因。

moduleLabels = net$colors #net$colors显示了数据集中的每个基因被分配到什么组中，不同的组用不同的数字表示，比如1表示组1，2表示组2。。。
moduleColors = labels2colors(moduleLabels)
# Plot the dendrogram and the module colors underneath

pdf("4.pdf")
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
dev.off()

# module eigengene, 可以绘制线图，作为每个模块的基因表达趋势的展示
MEs = net$MEs #net＄MEs：网络中每个模块的特征值。


library(stringi)
MEs_col = MEs
colnames(MEs_col) = paste0("ME", labels2colors(
  as.numeric(str_replace_all(colnames(MEs),"ME",""))))
library(stringr)
library(stringi)
MEs_col = orderMEs(MEs_col)

# 根据基因表达量进行聚类所得到的各模块间的相关性图
# marDendro/marHeatmap 设置下、左、、上右的边距
pdf("5.pdf")
plotEigengeneNetworks(MEs_col, "Eigengene adjacency heatmap", 
                      marDendro = c(3,3,2,4),
                      marHeatmap = c(3,4,2,2), plotDendrograms = T, 
                      xLabelsAngle = 90)
dev.off()

# TOM = TOMsimilarityFromExpr(dataExpr, power=power, corType=corType, networkType=type)
load(net$TOMFiles[1], verbose=T)
TOM <- as.matrix(TOM)

# 如果有表型数据，也可以跟ME数据放一起，一起出图
dissTOM = 1-TOM
# Transform dissTOM with a power to make moderately strong 
# connections more visible in the heatmap
plotTOM = dissTOM^7
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA
# Call the plot function

png("6.png")
#####这一部分特别耗时!!!!!行列同时做层级聚类
TOMplot(plotTOM, net$dendrograms, moduleColors, 
        main = "Network heatmap plot, all genes")
dev.off()





#dataTraits=read.csv(file = 'C:/Users/23997/Desktop/DK WGCNA zyy/4.干旱处理下表型数据.csv',row.names = 1)

dataTraits=read.csv(file = '/public/home/fengting/work/zyyy/wgcna/DK-WGCNA-feng/Hainan-WGCNA-feng/4.Drought-Phenotype.csv',row.names = 1)

#注意，上面这一步骤，输入的表型文件的“行数目”要与基因型文件的“列数目”要对应。也就是
#说基因型文件有多少列，表型文件就要有多少行。
#注意，这里不要求  基因型文件每一列对应的名字要与 表型文件每一行的名字一一相对应。
head(dataTraits)
dataTraits=dataTraits[,]


#识别关键模块
#定义基因和样本的数量
nGenes = ncol(dataExpr); 
nSamples = nrow(dataExpr);
#用颜色标签重新计算MEs(计算给定的单个数据集中模块的特征元素（第一主成分)
MEs0 = moduleEigengenes(dataExpr, moduleColors)$eigengenes #列名模块颜色
MEs = orderMEs(MEs0)
head(MEs)
head(dataTraits)
#View(MEs)
moduleTraitCor = cor(MEs, dataTraits, use = "p"); #查看模块与特征信息的相关性
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);#计算给定相关性的t test渐近p值
sizeGrWindow(10,6)
# Will display correlations and their p-values 
textMatrix = paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")", sep = ""); 
dim(textMatrix) = dim(moduleTraitCor) 

pdf("8.pdf",height=25,width=25)
par(mar = c(6, 8.5, 3, 3)); 



# Display the correlation values within a heatmap plot 
labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(dataTraits), 
               yLabels = names(MEs), ySymbols = names(MEs), colorLabels = FALSE, 
               colors = greenWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE, 
               cex.text = 0.5, zlim = c(-1,1), main = paste("Module-trait relationships"))


V2=as.data.frame(dataTraits$V2);

EFF=as.data.frame(dataTraits$LRS)

#这里把“  dataTraits$  ”后面的 “N—1”改成表型每一列的名称。然后跑完代码，出散点图。
#然后回到这一步，再改成“N—2”依次改完，表型的名称

dataTraits$screening.efficiency

dev.off()


pdf('PH.pdf',)
# names (colors) of the modules 
modNames = substring(names(MEs), 3)

#基因和模块的相关系数
geneModuleMembership = as.data.frame(cor(dataExpr, MEs, use = "p"));
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples));
names(geneModuleMembership) = paste("MM", modNames, sep=""); names(MMPvalue) = paste("p.MM", modNames, sep="");

#基因和性状的关系
geneTraitSignificance = as.data.frame(cor(dataExpr, EFF, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.", names(EFF), sep="");
names(GSPvalue) = paste("p.GS.", names(EFF), sep="");


#模内分析：鉴定具有高GS和MM的基因
#使用GS和MM度量，我们可以鉴定出对EFF以及在感兴趣的模块中具有较高模块成员性具有重要意义的基因。 
#例如，我们看一下与EFF关联最高的黄绿色模块。 我们在黄绿色模块中绘制了基因重要性与模块成员关系的散点图。 在此模块中，GS和MM之间存在高度显着的相关性。

#module = "greenyellow"
# column = match(module, modNames); 
moduleGenes = moduleColors==module;


# Select module
module = "brown";#这里要每一次都把模块-性状图里的对应的“最红色部分模块”对应的颜色改到这里，比如对应的是MEbrown，那就改成brown
# Select module probes
probes = colnames(dataExpr) 
inModule = (moduleColors==module);
modProbes = probes[inModule];
FilterGenes= abs(geneTraitSignificance[moduleGenes, 1])> .8 & abs(geneModuleMembership[moduleGenes, column])>.8  
modProbes[FilterGenes]


module = "brown"#这里要每一次都把模块-性状图里的对应的“最红色部分模块”对应的颜色改到这里，比如对应的是MEbrown，那就改成brown
column = match(module, modNames); 
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7); par(mfrow = c(1,1)); 
#画散点图
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   abs(geneTraitSignificance[moduleGenes, 1]), 
                   xlab = paste("Module Membership in", module, "module"), 
                   ylab = "Gene significance for LRS", #把"Gene significance for" 后面的 改成要研究的表型，如PH
                   main = paste("Module membership vs. gene significance\n"), 
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
dev.off()


options(stringsAsFactors = FALSE)
allowWGCNAThreads()
nGenes = ncol(dataExpr)
nSamples = nrow(dataExpr)
# 指定datTrait中感兴趣的一个性状，这里选择GY
GY = as.data.frame(dataTraits$GY)
names(GY) = "GY"
#  各基因模块的名字（颜色）modNames = substring(names(MEs), 3)
# 计算MM的P值
geneModuleMembership = as.data.frame(cor(dataExpr, MEs, use = "p"))

modNames = substring(names(MEs), 3)

MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))

names(geneModuleMembership) = paste("MM", modNames, sep="")
names(MMPvalue) = paste("p.MM", modNames, sep="")
# 计算性状和基因表达量之间的相关性（GS）

geneTraitSignificance = as.data.frame(cor(dataExpr, GY, use = "p"))

GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))

names(geneTraitSignificance) = paste("GS.", names(GY), sep="")
names(GSPvalue) = paste("p.GS.", names(GY), sep="")
module = "lightyellow"

column = match(module, modNames); 

moduleGenes = moduleColors==module;

verboseScatterplot(abs(geneModuleMembership[moduleGenes, 
column]),abs(geneTraitSignificance[moduleGenes, 1]), xlab = 
 paste("Module Membership in", module, "module"), ylab = "Gene 
   significance for TL", main = paste("Module membership 
   vs. gene significance"), pch = 20,col="grey")
abline(h=0.2,v=0.8,col="red",lwd=1.5)

table(moduleGenes)

light_yellow_module<-as.data.frame(dimnames(data.frame(dataExpr))[[2]][moduleGenes]) 
names(light_yellow_module)="genename"


MM<-abs(geneModuleMembership[moduleGenes,column])
GS<-abs(geneTraitSignificance[moduleGenes, 1])
c<-as.data.frame(cbind(MM,GS))
rownames(c)=light_yellow_module$genename

light_yellow_hub <-abs(c$MM)>0.8&abs(c$GS)>0.2

table(light_yellow_hub)


##dimnames(data.frame(dataExpr))[[2]][light_yellow_hub]

write.csv(light_yellow_hub, "hubgene_Normal_light_yellow.csv")



light_yellow_hub<-read.csv("hubgene_Normal_light_yellow.csv")
# hub基因和module内全部基因进行匹配，匹配成功返回1，没有匹配到的返回0
# match<- light_yellow_module$genename %in% light_yellow_hub$x
# 将匹配信息添加到散点图矩阵最后一列
c$group<-match
c$group<-light_yellow_hub
write.csv(c, "hubgene_Normal_light_yellow.csv")



library(ggplot2)
pdf("MM vs. GS_light_yellow_GY.pdf",width = 7,height = 7)
ggplot(data=c, aes(x=MM, y=GS, color=group))+geom_point(size=1.5)+scale_colour_manual(values=c("grey60", "#DE6757"))+ theme_bw()+  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())+  labs(x="Module Membership in module", y="Gene significance ",title = "Module membership vs. gene significance ")+theme(axis.title.x =element_text(size=14), axis.title.y=element_text(size=14),axis.text = element_text(size = 12),axis.text.x = element_text(colour = "black"),axis.text.y = element_text(colour = "black"),plot.title = element_text(hjust = 0.5,size = 16,face = "bold"),plot.margin = unit(rep(2,4),'lines')) +theme(legend.position = 'none')+geom_hline(aes(yintercept=0.2),colour="#5B9BD5",lwd=1,linetype=5)+geom_vline(aes(xintercept=0.8),colour="#5B9BD5",lwd=1,linetype=5)
dev.off()

GO分析网站：
Functional Annotation - Comprehensive Annotation of Rice Multi-Omics (sibs.ac.cn)

参考：
WGCNA得到模块之后如何筛选模块里面的hub基因 - 简书 (jianshu.com)
WGCNA（6）：筛选 hub gene - 简书 (jianshu.com)