WGCNA构建共表达网络问题，相关性矩阵出不来
WGCNA在进行TOMsimilarityFromExpr操作时，得到的TOM数据框出了对角线是1，其余全都是NaN，可能是哪里出了问题？chatGPT给出的几种方案都试过，没效果，下面贴出代码和结果，过程中没有报错
# 加载必须的包并做参数设置
library(MASS)
library(class)
library(cluster)
library(impute)
library(Hmisc)
library(WGCNA)
options(stringsAsFactors = F)
#load file
mycounts <- read.csv("2_4_expression count.csv", header = T)
###save thequestion!!!
is_duplicate <- duplicated(mycounts$gene_id)
mycounts <- mycounts[!is_duplicate, ]
#把格式改好
rownames(mycounts) <- mycounts[ , 1]
colnames(mycounts) <- mycounts[ , 1]
mycounts <- mycounts[ , -c(1)]
#Look at the first six lines and view the depth of data
head(mycounts)
dim(mycounts)
d
#去掉0行
mycounts_1 <- mycounts[rowSums(mycounts) != 0,]
dim(mycounts_1)
######## 从表达量矩阵中提取这些基因的表达量,将读取的csv文件变为向量，再提取，然后看数据深度对不对
common_genes <- read.csv("2_4_common_genes.csv", stringsAsFactors = T)
common_genes <- common_genes[ , 1]
selected_expr <- mycounts_1[common_genes, ]
dim(selected_expr)
# 数据初步处理
# 提取出表达数据
datExpr <- as.data.frame(t(selected_expr))
# 样本聚类检查离群值（就是树上特别不接近的）
gsg <- goodSamplesGenes(datExpr, verbose=3)
# 是true的话说明所有genes都通过了筛选
gsg$allOK
if (!gsg$allOK){
if (sum(!gsg$goodGenes)>0)
printFlush(paste("Removing genes:", paste(names(datExpr)[!gsg$goodGenes], collapse=", ")))
if (sum(!gsg$goodSamples)>0)
printFlush(paste("Removing samples:", paste(rownames(datExpr)[!gsg$goodSamples], collapse=", ")))
datExpr0 <- datExpr[gsg$goodSamples, gsg$goodGenes]
}
# 画图来看是否有特别离群的
sampleTree <- hclust(dist(datExpr0), method="average")
sizeGrWindow(12, 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)
abline(h=8000, col="red")
###求软阈值
# 允许使用最大线程
allowWGCNAThreads()
# 软阈值的预设范围
powers <- c(c(1:10), seq(from=12, to=20, by=2))
# 自动计算推荐的软阈值
sft <- pickSoftThreshold(datExpr, powerVector=powers, verbose=1, networkType="unsigned")
# 推荐值。如果是NA，就需要画图来自己挑选
sft$powerEstimate
# 作图
sizeGrWindow(9, 5)
par(mfrow=c(1, 2))
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")
abline(h=0.80, 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")
# 如果是NA，此时就需要根据图来自己指定数值了
sft$powerEstimate <- 16.0
###最后构建拓扑重叠矩阵(TOM矩阵)，目的是引入中间节点。因为构建邻接矩阵采用pearson相关，只描述线性关系，可能会漏掉一些非线性相关的基因
# 构建网络
# deepSplit是生成模块的梯度，从[0:4]中选择，越大模块越多
# minModuleSize最小模块的基因数目
# mergeCutHeight是合并相似的模块的合并系数，是通过主成分分析出来的
# mumericLabels 已数字命名模块
# nThreads 线程数，当设置为0时使用最大线程
powerEstimate <- as.numeric(sft$powerEstimate)
datExpr <- lapply(datExpr, as.numeric)
datExpr <- do.call(cbind, datExpr)
datExpr <- as.numeric(datExpr)
powerEstimate <- as.numeric(powerEstimate)
net <- blockwiseModules(datExpr, corType="pearson",
power=powerEstimate,
TOMType="unsigned", saveTOMs=TRUE, saveTOMFileBase="femaleMouseTOM",
deepSplit=2, minModuleSize=30,
reassignThreshold=0, mergeCutHeight=0.25,
numericLabels=TRUE, pamRespectsDendro=FALSE, nThreads=0,
verbose=3)
# 查看每个模块的基因数，其中0模块下为没有计算进入模块的基因数
table(net$colors)
sizeGrWindow(12, 9)
# 把模块编号转成颜色
mergedColors <- labels2colors(net$colors)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
"Module colors",
dendroLabels=F, hang=0.03,
addGuide=T, guideHang=0.05)
# 计算模块特征向量MEs
moduleLabels <- net$colors
moduleColors <- labels2colors(net$colors)
MEs <- net$MEs
geneTree <- net$dendrograms[[1]]
# 保存数据
save(MEs, moduleLabels, moduleColors, geneTree,
file="networkConstruction-auto.RData")
#########做主成分分析，用PC1代表模块的指标（用模块特征向量ME表示）。计算ME与性状的相关系数，再将相关系数用热图表示出来,这一步似乎没有做到，原因是没有足够的性状指标，年龄、性别等
nGenes <- ncol(datExpr)
nSamples <- nrow(datExpr)
# 重新计算MEs
MEs0 <- moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs <- orderMEs(MEs0)
# 主成分向量和性状的关联，pearson校正
datTraits <- t(selected_expr)
moduleTraitCor <- cor(MEs, datTraits, use="p")
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples)
# 画出模块和性状的关联热图
sizeGrWindow(10, 6)
# 把校正值和p值写在一起
textMatrix <- paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep="")
dim(textMatrix) <- dim(moduleTraitCor)
par(mar=c(6, 8.5, 3, 3))
##### 画图
#MEs 是 "Module Eigengenes" 的缩写。在基因共表达网络分析（WGCNA）中，一个"模块"是一个高度共表达的基因集，每个模块可以用该模块内的基因表达数据的第一主成分来代表，这就是所谓的"模块特征基因"（Module Eigengene）
#模块和基因的对应关系被保存在 net$colors 中。每个模块被赋予一个颜色作为标签，net$colors 是一个向量，长度等于基因的数量，每个基因被赋予其所属模块的颜色
labeledHeatmap(Matrix=moduleTraitCor,
xLabels=rep("", ncol(moduleTraitCor)),
yLabels=names(MEs),
ySymbols=names(MEs),
colorLabels=F,
colors=greenWhiteRed(50),
textMatrix=emptyTextMatrix,
setStdMargins=F,
cex.text=0.5,
zlim=c(-1,1),
main=paste("Module-trait relationships"))
#得到相关性最高的模块后，计算模块内单独的基因与模块的相关性（MM），也可以计算单独的基因与性状的相关性（GS）
## 计算模组
modNames <- substring(names(MEs), 3)
# 计算MM
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="")
# 计算GS
geneTraitSignificance <- as.data.frame(cor(datExpr, datTraits, use="p"))
GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
names(geneTraitSignificance) <- paste("GS.", names(datTraits), sep="")
names(GSPvalue) <- paste("p.GS.", names(datTraits), sep="")
# 输出
geneInfo <- cbind(geneModuleMembership, MMPvalue, geneTraitSignificance, GSPvalue)
write.table(geneInfo, file="geneInfo.txt", sep="\t", quote=F)
# power是之前的软阈值
# 然后再试一次
TOM <- TOMsimilarityFromExpr(datExpr, power = sft$powerEstimate)
dissTOM <- 1 - TOM
【以下回答由 GPT 生成】
问题: 您好！我可以帮助您优化您的问题。在您的问题中，您似乎遇到了WGCNA构建共表达网络时相关性矩阵无法生成的问题。为了更好地理解您的问题，请提供以下信息：
- 您使用的具体方法是什么？
- 请附上您的代码和结果，以便更好地理解问题。
- 是否有任何错误消息或警告？如果有，请提供具体的信息。
请提供这些信息，以便我能更好地帮助您优化问题。谢谢！
如果你已经解决了该问题, 非常希望你能够分享一下解决方案, 写成博客, 将相关链接放在评论区, 以帮助更多的人 ^-^