GEO生信数据挖掘（九）WGCNA分析

本文主要是介绍GEO生信数据挖掘（九）WGCNA分析，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

第六节，我们使用结核病基因数据，做了一个数据预处理的实操案例。例子中结核类型，包括结核，潜隐进展，对照和潜隐，四个类别。第七节延续上个数据，进行了差异分析。第八节对差异基因进行富集分析。本节进行WGCNA分析。

加载数据，进行聚类

初次聚类观察

自己定义红线位置，进行切割划分

载入性状数据

增加形状信息后，再次聚类

网络构建

选取soft-thresholding powers

基于tom的差异的基因聚类，绘制聚类树

根据聚类情况，设置颜色

计算eigengenes

模块的自动合并

模块与临床形状的关系热图（关键数据）

红色模块样本表达情况（相关性大）

产生了很多数据（各个模块的和临床性状）

后续挖掘核心基因时，需要用到Cytoscape，生成绘图所需要的数据

加载数据，进行聚类

library(WGCNA)
#读取目录名称，方便复制粘贴
dir()
#加载数据
load('DEG_TB_LTBI_step13.Rdata')#这里行为样品名，列为基因名
################################样品聚类#################### 
datExpr = t(dataset_TB_LTBI_DEG)
#初次聚类
sampleTree = hclust(dist(datExpr), method = "average")
# Plot the sample tree: Open a graphic output window of size 20 by 15 inches
# The user should change the dimensions if the window is too large or too small.
sizeGrWindow(12,9)
#pdf(file='sampleCluestering.pdf',width = 12,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)#结果图片自己导出PDF，文件名=1_sampleClustering.pdf### Plot a line to show the cut
abline(h = 87, col = "red")##剪切高度不确定，故无红线dev.off()

初次聚类观察

自己定义红线位置，进行切割划分

本例发现右侧有些样本孤立，适合被剔除，设置红线87切割。

左侧也被切成两块，需要做处理，保留。


### Determine cluster under the line
clust = cutreeStatic(sampleTree, cutHeight = 87, minSize = 10)
table(clust)
#clust
#0  1  2 
#5 57 40### clust 1 contains the samples we want to keep.
keepSamples = (clust==1|clust==2)
datExpr0 = datExpr[keepSamples, ]
dim(datExpr0) #[1]   97 2813
#保存数据
save(datExpr0,file='datExpr0_cluster_filter.Rdata')

载入性状数据

匹配样本名称，性状数据与表达数据保证一致


#################### 载入性状数据###########################
#加载自己的性状数据
load('design_TB_LTBI.Rdata')
traitData=design
#Loading clinical trait data
#traitData = read.table("trait_D.txt",row.names=1,header=T,comment.char = "",check.names=F)########trait file name can be changed######性状数据文件名，根据实际修改，如果工作路径不是实际性状数据路径，需要添加正确的数据路径
dim(traitData)
#names(traitData)
# remove columns that hold information we do not need.
#allTraits = traitData
dim(traitData)
names(traitData)# Form a data frame analogous to expression data that will hold the clinical traits.
fpkmSamples = rownames(datExpr0)
traitSamples =rownames(traitData)
#匹配样本名称，性状数据与表达数据保证一致
traitRows = match(fpkmSamples, traitSamples)
datTraits = traitData[traitRows,]
rownames(datTraits) 
collectGarbage()

增加形状信息后，再次聚类


# Re-cluster samples
sampleTree2 = hclust(dist(datExpr0), method = "average")
# Convert traits to a color representation: white means low, red means high, grey means missing entry
traitColors = numbers2colors(datTraits, signed = FALSE)
# Plot the sample dendrogram and the colors underneath.#sizeGrWindow(20,20)
##pdf(file="2_Sample dendrogram and trait heatmap.pdf",width=12,height=12)
plotDendroAndColors(sampleTree2, traitColors,groupLabels = names(datTraits),main = "Sample dendrogram and trait heatmap")dev.off()

下方红色，大致分成了两类，效果不错。