Entering edit mode
zen
•
0
@zen-22107
Last seen 4.0 years ago
Hello,
I am doing network analysis using WGCNA. I have successfully created the modules but got an error after relating it to the phenotype. I am using following code:
library(WGCNA)
GE.adjusted=read.csv2("Expression_File.csv", sep = ";", dec = ",", header = TRUE)
rownames(GE.adjusted)=GE.adjusted[,1]
GE.adjusted=GE.adjusted[,-1]
GE.adjusted=t(GE.adjusted)
Phenotype <- read.csv("Phenotype_File.csv", sep=";")
a=Phenotype[,1]
Phenotype=Phenotype[,-1]
Phenotype=as.matrix(Phenotype,ncol=5)
rownames(Phenotype)=a
colnames(Phenotype)=c("A", "B", "C", "D")
gsg = goodSamplesGenes(GE.adjusted, verbose = 3)
gsg$allOK
if (!gsg$allOK)
{
if (sum(!gsg$goodGenes)>0)
printFlush(paste("Removing genes:", paste(names(GE.adjusted)[!gsg$goodGenes], collapse = ", ")));
if (sum(!gsg$goodSamples)>0)
printFlush(paste("Removing samples:", paste(rownames(GE.adjusted)[!gsg$goodSamples], collapse = ", ")));
datExpr0 = GE.adjusted[gsg$goodSamples, gsg$goodGenes]
}
save(datExpr0, file="MissingData.RData")
sampleTree = hclust(dist(datExpr0), method = "average");
pdf("sampleClustering.pdf", width = 30, height = 30);
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
cex.axis = 1.5, cex.main = 2)
dev.off()
pdf("sampleClustering2.pdf", width = 30, height = 30);
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
cex.axis = 1.5, cex.main = 2)
abline(h = 100000, col = "red");
dev.off()
clust = cutreeStatic(sampleTree, cutHeight = 100000, minSize = 10)
table(clust)
keepSamples = (clust==1)
datExpr = datExpr0[keepSamples, ]
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
SamplesG = rownames(datExpr);
SamplesP = rownames(Phenotype);
traitRows = match(SamplesG, SamplesP);
datTraits = as.matrix(Phenotype[traitRows, ]);
collectGarbage();
save(datExpr, datTraits, file = "dataInput.RData")
#########Modules
options(stringsAsFactors = FALSE)
enableWGCNAThreads()
lnames = load(file = "dataInput.RData");
net_unsigned = blockwiseModules(datExpr, power = 6,
TOMType = "signed", minModuleSize = 30, maxBlockSize = 35000,
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = TRUE, pamRespectsDendro = FALSE,
saveTOMs = TRUE,
saveTOMFileBase = "TOM_signed",
verbose = 5)
pdf("Dendogram_Modules_signed.pdf", width = 30, height = 30);
mergedColors = labels2colors(net_unsigned$colors)
plotDendroAndColors(net_unsigned$dendrograms[[1]], mergedColors[net_unsigned$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
dev.off()
moduleLabels = net_unsigned$colors
moduleColors = labels2colors(net_unsigned$colors)
MEs = net_unsigned$MEs;
geneTree = net_unsigned$dendrograms[[1]];
save(MEs, moduleLabels, moduleColors, geneTree,
file = "networkConstruction-auto.RData")
#########Phenotype x Module Correlation
options(stringsAsFactors = FALSE);
load(file = "dataInput.RData");
load(file = "networkConstruction-auto.RData");
nGenes = ncol(datExpr);
nSamples = nrow(datExpr);
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
pdf("HeatMap.pdf", width = 7, height = 20);
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = names(datTraits),
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"))
At this step, I got error:
Error in labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits), :
Length of 'xLabels' must equal the number of columns in 'Matrix.'
Can you please help me to solve this? Thank you.
Thank you, I have tried your suggestion. However, I still get the same error.
In that case you probably dropped the column names somewhere in your code. Check
colnames(datTraits)
; if the result is empty or not what you expect, check how you created it.You are right. With
colnames(datTraits)
, I getNULL
. I useddatTraits = as.matrix(Phenotype[traitRows, ]);
. I am not able to find where the problem is.Try
datTraits = as.matrix(Phenotype[traitRows, , drop = FALSE])
I tried it and got this error:
Your
moduleTraitCor
may not be a matrix. Make sure it is a matrix (rundim(moduleTraitCor)
) and if it isn't, either useWGCNA::cor
instead ofcor
or manually assign dimensions tomoduleTraitCor
asdim(moduleTraitCor) = c(ncol(MEs), ncol(datTraits))
.Thank you! I tried using both of your suggestions but still got the same error. And I get this warning:
dim(moduleTraitCor)
gave33 1
where I have 33 modules for 1 phenotype.My phenotype data have some missing values as
"NA"
.Thank you, I have tried your suggestion. However, I still get the same error.