WGCNA module trait relationship issue
1
0
Entering edit mode
zen • 0
@zen-22107
Last seen 3.5 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.

wgcna • 3.2k views
ADD COMMENT
0
Entering edit mode
@peter-langfelder-4469
Last seen 1 day ago
United States

In the last function call, replace xLabels = names(datTraits) by xLabels = colnames(datTraits). In the tutorial, datTraits is a data.frame, whereas in your code it is a matrix which does not normally have names, only colnames.

ADD COMMENT
0
Entering edit mode

Thank you, I have tried your suggestion. However, I still get the same error.

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

You are right. With colnames(datTraits), I get NULL. I used datTraits = as.matrix(Phenotype[traitRows, ]);. I am not able to find where the problem is.

ADD REPLY
0
Entering edit mode

Try datTraits = as.matrix(Phenotype[traitRows, , drop = FALSE])

ADD REPLY
0
Entering edit mode

I tried it and got this error:

Error in colorMatrix[, c] : incorrect number of dimensions
In addition: Warning message:
In greenWhiteRed(50) :
  WGCNA::greenWhiteRed: this palette is not suitable for people
with green-red color blindness (the most common kind of color blindness).
Consider using the function blueWhiteRed instead.
ADD REPLY
0
Entering edit mode

Your moduleTraitCor may not be a matrix. Make sure it is a matrix (run dim(moduleTraitCor)) and if it isn't, either use WGCNA::cor instead of cor or manually assign dimensions to moduleTraitCor as dim(moduleTraitCor) = c(ncol(MEs), ncol(datTraits)).

ADD REPLY
0
Entering edit mode

Thank you! I tried using both of your suggestions but still got the same error. And I get this warning:

moduleTraitCor = WGCNA::cor(MEs, datTraits, use = "p");
Warning message:
In storage.mode(y) <- "double" : NAs introduced by coercion

dim(moduleTraitCor) gave 33 1 where I have 33 modules for 1 phenotype.

My phenotype data have some missing values as "NA".

ADD REPLY
0
Entering edit mode

Thank you, I have tried your suggestion. However, I still get the same error.

ADD REPLY

Login before adding your answer.

Traffic: 447 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6