library(WGCNA)
setwd("D:\\julei")
enableWGCNAThreads()
## 读取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)
dim(WGCNA.fpkm)
names(WGCNA.fpkm)
datExpr0 = as.data.frame(t(WGCNA.fpkm[,-1]))
names(datExpr0) = WGCNA.fpkm$sample;##########如果第一行不是ID命名，就写成fpkm[,1]
rownames(datExpr0) = names(WGCNA.fpkm[,-1])
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]
}
##filter
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="mRNA.filter.txt",
row.names=F, col.names=T,quote=FALSE,sep="\t")
sampleTree = hclust(dist(datExpr0), method = "average")
pdf(file = "1.sampleClustering.pdf", width = 15, height = 8)
par(cex = 0.6)
par(mar = c(0,6,6,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 2,
cex.axis = 1.5, cex.main = 2)
### Plot a line to show the cut
#abline(h = 180, col = "red")##剪切高度不确定，故无红线
dev.off()
traitData = read.table("TraitData.txt",row.names=1,header=T,comment.char = "",check.names=F)
allTraits = traitData
dim(allTraits)
names(allTraits)
## 形成一个类似于表达数据的数据框架
fpkmSamples = rownames(datExpr0)
traitSamples =rownames(allTraits)
traitRows = match(fpkmSamples, traitSamples)
datTraits = allTraits[traitRows,]
rownames(datTraits)
collectGarbage()
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)
pdf(file="2.Sample_dendrogram_and_trait_heatmap.pdf",width=20,height=12)
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()
# 设置soft-thresholding powers的数量
powers = c(1:30)
sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)
cex1=0.9
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()
softPower =sft$powerEstimate
adjacency = adjacency(datExpr0, power = softPower)
TOM = TOMsimilarity(adjacency);
dissTOM = 1-TOM
geneTree = hclust(as.dist(dissTOM), method = "average");
pdf(file="4_Gene clustering on TOM-based dissimilarity.pdf",width=24,height=18)
plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04)
dev.off()
minModuleSize = 30
# Module identification using dynamic tree cut:
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
deepSplit = 2, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize);
table(dynamicMods)
# Convert numeric lables into colors
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
# Plot the dendrogram and colors underneath
#sizeGrWindow(8,6)
pdf(file="5_Dynamic Tree Cut.pdf",width=8,height=6)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors")
dev.off()
MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
MEs = MEList$eigengenes
# Calculate dissimilarity of module eigengenes
MEDiss = 1-cor(MEs);
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average")
# Plot the result
#sizeGrWindow(7, 6)
pdf(file="6_Clustering of module eigengenes.pdf",width=7,height=6)
plot(METree, main = "Clustering of module eigengenes",
xlab = "", sub = "")
######剪切高度可修改
MEDissThres = 0.4
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")
dev.off()
mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColors = merge$colors
# Eigengenes of the new merged modules:
mergedMEs = merge$newMEs
table(mergedColors)
在最后运行出现 Error in merge$colors : object of type 'closure' is not subsettable，有什么办法能解决这个问题吗，查了很多代码都要用到，没法不用，应该也不存在包没有引用和大小写之类的错误
原文来自https://blog.csdn.net/kanghua_du/article/details/129232234
直接的错误信息是说merge$colors那里，这个是封闭类型的数据，不能进行子集操作。说白了就是调用方式不对，请检查下你的变量名或函数名有没有写错，比如你的merge$colors这个变量哪里来的，你程序中都没有，这些都要检查下。
这个错误通常表示您尝试使用一个函数或闭包对象的返回值，但是该返回值并不是您期望的类型。
您可以尝试检查一下代码，查看是否有任何意外的赋值或类型转换，或者尝试重启R并重新运行您的代码。
如果问题仍然存在，您可以尝试使用调试器来诊断问题。在发生错误的行之前，添加以下代码：
options(error = recover)
这将允许您进入调试器并检查运行时环境中的变量和对象。您可以使用traceback()函数来获取有关错误的更多详细信息。
该回答引用ChatGPT：这个错误通常是由于在merge$colors中试图访问闭包对象而导致的。可能的原因是在之前的代码中已经定义了与merge名称相同的对象，从而覆盖了merge函数。
您可以尝试使用以下方法来解决此错误：
1.重命名所有与merge名称相同的变量和对象。
2.在调用merge函数时，将其指定为base包的函数，例如base::merge()。
3.使用detach(package:WGCNA)取消加载WGCNA软件包，然后再次加载并运行代码。
希望这些方法能够解决您的问题。
出现"Error in merge$colors : object of type 'closure' is not subsettable"错误，通常是因为在使用WGCNA包时，调用了merge函数而参数传递错误引起的。这个错误提示中的"closure"指的是函数，也就是说，出现了函数不应该出现的地方，导致了错误。
解决这个问题，可以尝试以下几个方法：
检查参数传递是否正确：在调用merge函数时，需要传递正确的参数，包括待合并的数据对象和颜色向量等。需要检查一下传递的参数是否正确。
检查版本是否兼容：WGCNA包的不同版本之间可能存在差异，导致函数调用出错。可以尝试卸载当前版本的WGCNA包，重新安装一个兼容的版本。
检查其他代码是否出错：错误提示中提到了merge函数，但是实际上这个错误可能是由其他代码引起的。需要细心地检查代码中是否存在其他类型的对象被错误地当做了向量或矩阵来处理，或者是否存在其他语法错误等。