WGCNA的使用

WGCNA的使用

WGCNA算是一种很经典的分析方法。可以将临床特征与基因联系起来，但是怎么联系的呢。这就涉及了相关性分析，也就是聚类。这里很重要，要将聚类理解成相关性分析，那么生信学起来就会容易很多。

后面有空我会写一篇关于聚类的文章，就会发现其实生信就是在聚类成的一类一类中去摸索我们想要的。

WGCNA将临床特征与基因联系起来的过程中，引入了一个新的东西，叫做模块。这个模块也就是将基因聚类

那这样我们就有几种研究方法呢。这就是简单的排列组合而已了。

可以去研究基因与模块的关系，那就是把基因聚类成模块。这部分做了吗？做了的，也就是使用软阈值构建模块可以去研究模块与临床特征的关系。也是做了的。也就是将模块与临床特征的关系以热图的形式展示出来，而且哈，去研究这两者的关系运用的是cor函数，这个函数特别重要，研究的是皮尔森相关性。可以去研究基因与临床特征的关系，做了吗，也是做了的。运用的也是cor函数。这些研究完还有别的吗，还有！就是将三者放在一起研究。怎么研究呢？ 可以将基因与模块的关系矩阵与模块与临床特征的关系矩阵放在一起可以将基因与模块的关系矩阵与基因与临床特征的关系矩阵放在一起，这个研究的多一点，因为涉及到基因。可以将模块与临床特征的关系矩阵与基因与临床特征的关系矩阵放在一起。

回头看一下，是不是已经有6种研究方法了，其实还不止，稍微将一个变量变一下，又有几种方法出来。应该反复去阅读思考上面的内容。

了解完WGCNA能做的事情之后那学起来就会容易很多。

下面一步一步讲

1.处理数据。下面是详细的代码

rm(list = ls())
library(WGCNA)# 读取基因表达矩阵数据
fpkm = read.table("data/fpkm.txt", header = T, row.names = 1, check.names = F)
head(fpkm)### 选取基因方法 ##### 第一种，通过标准差选择
## 计算每个基因的标准差
fpkm_sd = apply(fpkm,1,sd)#1是对每一行，2是对每一列
## 使用标准差对基因进行降序排序
fpkm_sd_sorted = order(fpkm_sd, decreasing = T)
## 选择前5000个标准差较大的基因
fpkm_num = fpkm_sd_sorted[1:5000]
## 从表达矩阵提取基因
fpkm_filter = fpkm[fpkm_num,]
## 对表达矩阵进行转置
WGCNA_matrix = t(fpkm_filter)#变成行名是样本，列名是基因
## 保存过滤后的数据
save(WGCNA_matrix, file = "data/Step01-fpkm_sd_filter.Rda")## 第二种，使用绝对中位差选择，推荐使用绝对中位差
WGCNA_matrix = t(fpkm[order(apply(fpkm,1,mad), decreasing = T)[1:5000],])#mad代表绝对中位差
save(WGCNA_matrix, file = "data/Step01-fpkm_mad_filter.Rda")## 第三种，使用全基因
WGCNA_matrix = t(fpkm)
save(WGCNA_matrix, file = "data/Step01-fpkm_allgene.Rda")### 去除缺失值较多的基因/样本 ###
rm(list = ls())
# 加载表达矩阵
load(file = "data/Step01-fpkm_mad_filter.Rda")
# 加载WGCNA包
library(WGCNA)
# 判断是否缺失较多的样本和基因
datExpr0 = WGCNA_matrix
gsg = goodSamplesGenes(datExpr0, verbose = 3)# 是否需要过滤，TRUE表示不需要，FALSE表示需要
gsg$allOK
# 当gsg$allOK为TRUE，以下代码不会运行，为FALSE时，运行以下代码过滤样本
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 = ", ")));# 提取保留的基因和样本datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}### 通过样本聚类识别离群样本，去除离群样本 ###
sampleTree = hclust(dist(datExpr0), method = "average");#使用hclust函数进行均值聚类
# 绘制样本聚类图确定离群样本
sizeGrWindow(30,9)
pdf(file = "figures/Step01-sampleClustering.pdf", width = 30, height = 9);
par(cex = 0.6);
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)
# 根据上图判断，需要截取的高度参数h
abline(h = 120, col = "red")#在120的地方画条线
dev.off()# 去除离群得聚类样本，cutHeight参数要与上述得h参数值一致
clust = cutreeStatic(sampleTree, cutHeight = 120, minSize = 10)
table(clust)
# clust
# 0   1  0就是要去除的，1就是要保存的
# 15 162# clust 1聚类中包含着我们想要得样本，将其提取出来
keepSamples = (clust==1)
datExpr = datExpr0[keepSamples, ]# 记录基因和样本数，方便后续可视化
nGenes = ncol(datExpr)#基因数
nSamples = nrow(datExpr)#样本数
save(datExpr, nGenes, nSamples,file = "data/Step01-WGCNA_input.Rda")

如果有一点R语言基础会发现上面代码就做了三件事，选择基因，过滤基因，基因聚类。为什么要做这些事情，下面来一件一件说。

选择基因中有三种方法。其实这个部分是很灵活的，三种方法都可以，如果计算机资源可以，那完全可以采用所有基因进行后续分析。

过滤基因。这相当于数据的清洗。

基因聚类。这一部分特别重要。有人说可以不做，但自我感觉重要，因为我们后面做的都是与聚类相关的，如果不做这一步直接做后面相关的聚类可信度会少很多。

筛选软阈值。为什么要筛选软阈值？那就需要了解一下软阈值的意义。这是WGCNA的底层思维。筛选到软阈值之后就能构建模块了，也就是构建共表达网络。有没有发现，这是一步一步递进的。

主要运用的函数是，这个函数有必要好好研究一下，源代码我会放在最后，有时间我再来写篇文章。

这个函数运行的过程我大致讲一下，要好好理解才行。

第一步

在得到一个m × n 表达矩阵（m个样本，n个基因）后，第一步是计算两基因在多样本的表达水平相关性（例如 ，值一般在-1 ~ 1之间），从而得到表达相似度矩阵 matix，简记为S。（在这里要思考一下为什么要计算表达水平相关性，举一反三，计算别的相关性可以吗？）

根据是否考虑相关性的正负性，有两种处理方式。

不要把sign和弄混了。前者表示相关的正负性，后者表示方向性。使用WGCNA建立的都是网络。

第二步

使用 将基因表达相似矩阵转为共表达邻接矩阵。根据 的特点分为hard与soft两种。

第三步

鉴定模块及相关分析

# 清空所有变量
rm(list = ls())
# 加载包
library(WGCNA)
# 允许多线程运行
enableWGCNAThreads()
# 加载表达矩阵
load("data/Step01-WGCNA_input.Rda")### 选择软阈值 ###
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# 进行网络拓扑分析
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)#β=power,就是软阈值# 可视化结果
sizeGrWindow(9, 5)
pdf(file = "figures/Step02-SoftThreshold.pdf", width = 9, height = 5);par(mfrow = c(1,2))#一个画板上，画两个图，一行两列
cex1 = 0.9;
# 无尺度网络阈值得选择
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],xlab="Soft Threshold (power)",#x轴ylab="Scale Free Topology Model Fit,signed R^2",type="n",#y轴main = paste("Scale independence"));#标题
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],labels=powers,cex=cex1,col="red");# 用红线标出R^2的参考值
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")dev.off()# 无尺度网络检验，验证构建的网络是否是无尺度网络
softpower=sft$powerEstimateADJ = abs(cor(datExpr,use="p"))^softpower#相关性取绝对值再幂次
k = as.vector(apply(ADJ,2,sum,na.rm=T))#对ADJ的每一列取和，也就是频次pdf(file = "figures/Step02-scaleFree.pdf",width = 14)
par(mfrow = c(1,2))hist(k)#直方图
scaleFreePlot(k,main="Check scale free topology")dev.off()### 一步构建网络 ###
net = blockwiseModules(datExpr, #处理好的表达矩阵power = sft$powerEstimate,#选择的软阈值TOMType = "unsigned", #拓扑矩阵类型，none表示邻接矩阵聚类，unsigned最常用，构建无方向minModuleSize = 30,#网络模块包含的最少基因数reassignThreshold = 0, #模块间基因重分类的阈值mergeCutHeight = 0.25,#合并相异度低于0.25的模块numericLabels = TRUE, #true，返回模块的数字标签 false返回模块的颜色标签pamRespectsDendro = FALSE,#调用动态剪切树算法识别网络模块后，进行第二次的模块比较，合并相关性高的模块saveTOMs = TRUE,#保存拓扑矩阵saveTOMFileBase = "data/Step02-fpkmTOM", verbose = 3)#0，不反回任何信息，＞0返回计算过程
# 保存网络构建结果
save(net, file = "data/Step02-One_step_net.Rda")# 加载网络构建结果
load(file = "data/Step02-One_step_net.Rda")
# 打开绘图窗口
sizeGrWindow(12, 9)
pdf(file = "figures/Step02-moduleCluster.pdf", width = 12, height = 9);
# 将标签转化为颜色
mergedColors = labels2colors(net$colors)
# 绘制聚类和网络模块对应图
plotDendroAndColors(dendro = net$dendrograms[[1]], #hclust函数生成的聚类结果colors = mergedColors[net$blockGenes[[1]]],#基因对应的模块颜色groupLabels = "Module colors",#分组标签dendroLabels = FALSE, #false,不显示聚类图的每个分支名称hang = 0.03,#调整聚类图分支所占的高度addGuide = TRUE, #为聚类图添加辅助线guideHang = 0.05,#辅助线所在高度main = "Gene dendrogram and module colors")
dev.off()# 加载TOM矩阵
load("data/Step02-fpkmTOM-block.1.RData")
# 网络特征向量
MEs = moduleEigengenes(datExpr, mergedColors)$eigengenes
# 对特征向量排序
MEs = orderMEs(MEs)
# 可视化模块间的相关性
sizeGrWindow(5,7.5);
pdf(file = "figures/Step02-moduleCor.pdf", width = 5, height = 7.5);par(cex = 0.9)
plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle = 90)
dev.off()## TOMplotdissTOM = 1-TOMsimilarityFromExpr(datExpr, power = sft$powerEstimate); #1-相关性=相异性nSelect = 400
# 随机选取400个基因进行可视化，设置seed值，保证结果的可重复性
set.seed(10);#设置随机种子数
select = sample(nGenes, size = nSelect);#从5000个基因选择400个
selectTOM = dissTOM[select, select];#选择这400*400的矩阵
# 对选取的基因进行重新聚类
selectTree = hclust(as.dist(selectTOM), method = "average")#用hclust重新聚类
selectColors = mergedColors[select];#提取相应的颜色模块
# 打开一个绘图窗口
sizeGrWindow(9,9)
pdf(file = "figures/Step02-TOMplot.pdf", width = 9, height = 9);
# 美化图形的设置
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
# 绘制TOM图
TOMplot(plotDiss, #拓扑矩阵，该矩阵记录了每个节点之间的相关性selectTree, #基因的聚类结果selectColors, #基因对应的模块颜色main = "Network heatmap plot, selected genes")
dev.off()

多步法构建模块

# 清空环境变量
rm(list = ls())
#加载R包
library(WGCNA)
# 允许多线程运行
enableWGCNAThreads()
# 加载表达矩阵
load("data/Step01-WGCNA_input.Rda")## 选择软阈值
sizeGrWindow(9,5);
par(mfrow = c(1,2));
powers = c(c(1:10), seq(from = 12, to=20, by=2))
RpowerTable=pickSoftThreshold(datExpr, powerVector=powers)[[2]]
cex1=0.7
pdf(file="figures/Step03-softThresholding.pdf",width = 14)
par(mfrow=c(1,2))
plot(RpowerTable[,1], -sign(RpowerTable[,3])*RpowerTable[,2],xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",main = paste("Scale independence"))
text(RpowerTable[,1], -sign(RpowerTable[,3])*RpowerTable[,2], labels=powers,cex=cex1,col="red")abline(h=0.9,col="red")
plot(RpowerTable[,1], RpowerTable[,5],xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",main = paste("Mean connectivity"))
text(RpowerTable[,1], RpowerTable[,5], labels=powers, cex=cex1,col="red")
dev.off()## 无尺度网络验证
softpower=7ADJ = abs(cor(datExpr,use="p"))^softpower
k = as.vector(apply(ADJ,2,sum,na.rm=T))pdf(file = "figures/Step03-scaleFree.pdf",width = 14)
par(mfrow = c(1,2))
hist(k)
scaleFreePlot(k,main="Check scale free topology")
dev.off()## 计算邻接矩阵
adjacency = adjacency(datExpr,power=softpower)
## 计算TOM拓扑矩阵
TOM = TOMsimilarity(adjacency)
## 计算相异度
dissTOM = 1- TOM#模块初步聚类分析
library(flashClust)
geneTree = flashClust(as.dist(dissTOM),method="average")
#绘制层次聚类树
pdf(file = "figures/Step03-GeneClusterTOM-based.pdf")
plot(geneTree,xlab="",sub="",main="Gene clustering on TOM-based",labels=FALSE,hang=0.04)
dev.off()#构建初步基因模块
#设定基因模块中至少30个基因
minModuleSize=30
# 动态剪切树识别网络模块
dynamicMods = cutreeDynamic(dendro = geneTree,#hclust函数的聚类结果distM = dissTOM,#deepSplit = 2,pamRespectsDendro = FALSE,minClusterSize = minModuleSize)#设定基因模块中至少30个基因
# 将标签转换为颜色
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)#看聚类到哪些模块，哪些颜色pdf(file="figures/Step03-DynamicTreeCut.pdf",width=9,height=5)
plotDendroAndColors(dendro = geneTree, colors = dynamicColors, groupLabels = "Dynamic Tree Cut",dendroLabels = FALSE, hang = 0.03,addGuide = TRUE,main = "Gene dendrogram and module colors")
dev.off()#前面用动态剪切树聚类了一些模块，现在要对这步结果进一步合并，合并相似度大于0.75的模块，降低网络的复杂度
#计算基因模块特征向量
MEList = moduleEigengenes(datExpr,colors=dynamicColors)#计算特征向量
MEs = MEList$eigengenes;#提取特征向量
MEDiss1 = 1-cor(MEs);#计算相异度
METree1 = flashClust(as.dist(MEDiss1),method="average")#对相异度进行flashClust聚类
#设置特征向量相关系数大于0.75
MEDissThres = 0.25;#相异度在0.25以下，也就是相似度大于0.75，对这些模块合并
#合并模块
merge = mergeCloseModules(datExpr, #合并相似度大于0.75的模块dynamicColors,cutHeight = MEDissThres, verbose=3)
mergedColors = merge$colorstable(dynamicColors)#动态剪切树的模块颜色
table(mergedColors)#合并后的模块颜色，可以看到从18个模块变成了14个模块mergedMEs = merge$newMEs#合并后的14个模块#重新命名合并后的模块
moduleColors = mergedColors;
colorOrder = c("grey",standardColors(50));
moduleLabels = match(moduleColors,colorOrder)-1;
MEs = mergedMEs;
MEDiss2 = 1-cor(MEs);#计算相异度
METree2 = flashClust(as.dist(MEDiss2),method="average");#对合并后的模块进行聚类#绘制聚类结果图
pdf(file="figures/Step03-MECombined.pdf",width=12,height=5)
par(mfrow=c(1,2))
plot(METree1,xlab="",sub="",main="Clustering of ME before combined")# METree1是动态剪切树形成的模块
abline(h=MEDissThres,col="red")#相异度为0.25
plot(METree2,xlab="",sub="",main="Clustering of ME after combined")# METree2是合并后的模块
dev.off()pdf(file="figures/Step03-MergedDynamics.pdf",width=10,height=4)
plotDendroAndColors(dendro = geneTree,#剪切树colors = cbind(dynamicColors,mergedColors),#将两种方法形成的模块颜色合并在一起groupLabels = c("Dynamic Tree Cut","Merged Dynamics"),dendroLabels = FALSE, hang = 0.03,addGuide=TRUE,guideHang=0.05,main="Gene Dendrogram and module colors")
dev.off()# 模块中基因数
write.table(table(moduleColors),"data/Step03-MEgeneCount.txt",quote = F,row.names = F)# 保存构建的网络信息
moduleColors=mergedColors
colorOrder=c("grey", standardColors(50))
moduleLabels=match(moduleColors, colorOrder)-1
MEs=mergedMEs
save(MEs, moduleLabels, moduleColors, geneTree, file="data/Step03-Step_by_step_buildnetwork.rda")## 绘制样本间的相关性
load("data/Step01-WGCNA_input.Rda")
load(file="data/Step03-Step_by_step_buildnetwork.rda")MEs = orderMEs(MEs)sizeGrWindow(5,7.5);
pdf(file = "figures/Step03-moduleCor.pdf", width = 5, height = 7.5);par(cex = 0.9)
plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle= 90)
dev.off()## TOMplotdissTOM = 1-TOMsimilarityFromExpr(datExpr, power = softpower); nSelect = 400
# 随机选取400个基因进行可视化，s设置seed值，保证结果的可重复性
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
# 对选取的基因进行重新聚类
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = mergedColors[select];
# 打开一个绘图窗口
sizeGrWindow(9,9)
pdf(file = "figures/Step03-TOMplot.pdf", width = 9, height = 9);
# 美化图形的设置
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
dev.off()

第四步

# 清空环境变量
rm(list = ls())
# 加载包
library(WGCNA)
# 加载表达矩阵
load("data/Step01-WGCNA_input.Rda")# 读入临床信息
clinical = read.table("data/clinical.txt",stringsAsFactors = TRUE, header = T,row.names = 1,na.strings = "",sep = "\t")
# 查看临床信息
head(clinical)
# 对表达矩阵进行预处理
datTraits = as.data.frame(do.call(cbind,lapply(clinical, as.numeric)))
rownames(datTraits) = rownames(clinical)# 对样本进行聚类
sampleTree2 = hclust(dist(datExpr), method = "average")# 将临床信息转换为颜色，白色表示低，红色表示高，灰色表示缺失
traitColors = numbers2colors(datTraits, signed = FALSE)pdf(file = "figures/Step04-Sample_dendrogram_and_trait_heatmap.pdf", width = 24);
# 样本聚类图与样本性状热图
plotDendroAndColors(sampleTree2, traitColors,groupLabels = names(datTraits), main = "Sample dendrogram and trait heatmap")
dev.off()#### 网络的分析
###### 基因模块与临床信息的关系
# 加载构建的网络
load(file = "data/Step03-Step_by_step_buildnetwork.rda")
# 对模块特征矩阵进行排序
MEs=orderMEs(MEs)
#计算模型特征矩阵和样本信息矩阵的相关度。
moduleTraitCor=cor(MEs, datTraits, use="p")
write.table(file="data/Step04-modPhysiological.cor.xls",moduleTraitCor,sep="\t",quote=F)
moduleTraitPvalue=corPvalueStudent(moduleTraitCor, nSamples)
write.table(file="data/Step04-modPhysiological.p.xls",moduleTraitPvalue,sep="\t",quote=F)#使用labeledHeatmap()将上述相关矩阵和p值可视化。
pdf(file="figures/Step04-Module_trait_relationships.pdf",width=9,height=7)
textMatrix=paste(signif(moduleTraitCor,2),"\n(",signif(moduleTraitPvalue,1),")",sep="")
dim(textMatrix)=dim(moduleTraitCor)
# 基因模块与临床信息相关性图
labeledHeatmap(Matrix=moduleTraitCor,#模块和表型的相关性矩阵，这个参数最重要，其他可以不变xLabels=colnames(datTraits),yLabels=names(MEs),ySymbols=names(MEs),colorLabels=FALSE,colors=blueWhiteRed(50),textMatrix=textMatrix,setStdMargins=FALSE,cex.text=0.5,cex.lab=0.5,zlim=c(-1,1),main=paste("Module-trait relationships"))
dev.off()## 单一模块与某一表型相关性
M_stage = as.data.frame(datTraits$M_stage)
# 分析自己感兴趣的临床信息，此处以M_stage为示例
names(M_stage) = "M_stage"
# 模块对应的颜色
modNames = substring(names(MEs), 3)# 计算基因模块特征
geneModuleMembership = as.data.frame(cor(datExpr, 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="");# 计算M分期基因特征显著性
geneTraitSignificance = as.data.frame(cor(datExpr, M_stage, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));# 对结果进行命名
names(geneTraitSignificance) = paste("GS.", names(M_stage), sep="")
names(GSPvalue) = paste("p.GS.", names(M_stage), sep="")# 设置需要分析的模块名称，此处为brown模块
module = "brown"
# 提取brown模块数据
column = match(module, modNames);
moduleGenes = moduleColors==module;# 可视化brown模块与M分期的相关性分析结果
sizeGrWindow(7, 7);
pdf(file="figures/Step04-Module_membership_vs_gene_significance.pdf")
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 M Stage",main = paste("Module membership vs. gene significance\n"),cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
dev.off()## 单一特征与所有模块关联分析
GS = as.numeric(cor(datTraits$M_stage,datExpr, use="p"))
GeneSignificance = abs(GS);
pdf(file="figures/Step04-Gene_significance_for_M_stage_across_module.pdf",width=9,height=5)
plotModuleSignificance(GeneSignificance, moduleColors, ylim=c(0,0.2), main="Gene significance for M stage across module");
dev.off()## 模块中的hub基因
#### 为每一个模块寻找hub基因
HubGenes <- chooseTopHubInEachModule(datExpr,#WGCNA分析输入的表达矩阵moduleColors)#模块颜色信息
# 保存hub基因结果
write.table (HubGenes,file = "data/Step04-HubGenes_of_each_module.xls",quote=F,sep='\t',col.names = F)#### 与某种特征相关的hub基因
NS = networkScreening(datTraits$M_stage,#M分期MEs,#datExpr)#WGCNA分析输入的表达矩阵
# 将结果写入到文件
write.table(NS,file="data/Step04-Genes_for_M_stage.xls",quote=F,sep='\t')## 模块GO/KEGG分析
# 加载R包
library(anRichment)
library(clusterProfiler)
##### GO分析
# 构建GO背景基因集
GOcollection = buildGOcollection(organism = "human")
geneNames = colnames(datExpr)
# 将基因SYMBOL转换为ENTREZID基因名
geneID = bitr(geneNames,fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db", drop = FALSE)
# 将基因名对应结果写入文件中
write.table(geneID, file = "data/Step04-geneID_map.xls", sep = "\t", quote = TRUE, row.names = FALSE)
# 进行GO富集分析
GOenr = enrichmentAnalysis(classLabels = moduleColors,#基因所在的模块信息identifiers = geneID$ENTREZID,refCollection = GOcollection,useBackground = "given",threshold = 1e-4,thresholdType = "Bonferroni",getOverlapEntrez = TRUE,getOverlapSymbols = TRUE,ignoreLabels = "grey");
# 提取结果，并写入结果到文件
tab = GOenr$enrichmentTable
names(tab)
write.table(tab, file = "data/Step04-GOEnrichmentTable.xls", sep = "\t", quote = TRUE, row.names = FALSE)# 提取主要结果，并写入文件
keepCols = c(1, 3, 4, 6, 7, 8, 13)
screenTab = tab[, keepCols]
# 小数位为2位
numCols = c(4, 5, 6)
screenTab[, numCols] = signif(apply(screenTab[, numCols], 2, as.numeric), 2)# 给结果命名
colnames(screenTab) = c("module", "GOID", "term name", "p-val", "Bonf", "FDR", "size")
rownames(screenTab) = NULL# 查看结果
head(screenTab)
# 写入文件中
write.table(screenTab, file = "data/Step04-GOEnrichmentTableScreen.xls", sep = "\t", quote = TRUE, row.names = FALSE)##### KEGG富集分析
# AnRichment没有直接提供KEGG数据的背景集
# 这里使用MSigDBCollection构建C2通路数据集
KEGGCollection = MSigDBCollection("data/msigdb_v7.1.xml", MSDBVersion = "7.1",organism = "human",excludeCategories = c("h","C1","C3","C4","C5","C6","C7")) 
# KEGG分析
KEGGenr = enrichmentAnalysis(classLabels = moduleColors,identifiers = geneID$ENTREZID,refCollection = KEGGCollection,useBackground = "given",threshold = 1e-4,thresholdType = "Bonferroni",getOverlapEntrez = TRUE,getOverlapSymbols = TRUE,ignoreLabels = "grey")
# 提取KEGG结果，并写入文件
tab = KEGGenr$enrichmentTable
names(tab)
write.table(tab, file = "data/Step04-KEGGEnrichmentTable.xls", sep = "\t", quote = TRUE, row.names = FALSE)# 提取主要结果并写入文件
keepCols = c(1, 3, 4, 6, 7, 8, 13)
screenTab = tab[, keepCols]
# 取两位有效数字
numCols = c(4, 5, 6)
screenTab[, numCols] = signif(apply(screenTab[, numCols], 2, as.numeric), 2)# 对结果表格进行重命名
colnames(screenTab) = c("module", "ID", "term name", "p-val", "Bonf", "FDR", "size")
rownames(screenTab) = NULL
# 查看结果
head(screenTab)
# 写入文件中
write.table(screenTab, file = "data/Step04-KEGGEnrichmentTableScreen.xls", sep = "\t", quote = TRUE, row.names = FALSE)### 输出cytoscape可视化
# 重新计算TOM，power值设置为前面选择好的
TOM = TOMsimilarityFromExpr(datExpr, power = 7)
# 输出全部网络模块
cyt = exportNetworkToCytoscape(TOM,edgeFile = "data/Step04-CytoscapeInput-edges-all.txt",#基因间的共表达关系nodeFile = "data/Step04-CytoscapeInput-nodes-all.txt",#weighted = TRUE,threshold = 0.1,nodeNames = geneID$SYMBOL,altNodeNames = geneID$ENTREZID,nodeAttr = moduleColors)
# 输出感兴趣网络模块
modules = c("brown", "red")
# 选择上面模块中包含的基因
inModule = is.finite(match(moduleColors, modules))
modGenes = geneID[inModule,]
# 选择指定模块的TOM矩阵
modTOM = TOM[inModule, inModule]# 输出为Cytoscape软件可识别格式
cyt = exportNetworkToCytoscape(modTOM,edgeFile = paste("data/Step04-CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),nodeFile = paste("data/Step04-CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),weighted = TRUE,threshold = 0.2,nodeNames = modGenes$SYMBOL,altNodeNames = modGenes$ENTREZID,nodeAttr = moduleColors[inModule])