Error in TOMplot (WGCNA)
1
0
Entering edit mode
@abhishek-singh-4725
Last seen 10 months ago
France

Hi,

I am getting this error while running TOMplot.

TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
Error in .heatmap(as.matrix(dissim), Rowv = as.dendrogram(dendro, hang = 0.1),  :
  row dendrogram ordering gave index of wrong length

geneTree and plotTOM were made exactly the way as described in the manual https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/index.html

 Here is the complete code:


library(WGCNA)
femData = read.csv("lg2FPKM.csv");
dim(femData);
names(femData);
datExpr0 = as.data.frame(t(femData[, -1]))
names(datExpr0) = femData$substanceBXH;
rownames(datExpr0) = names(femData)[-1];
gsg = goodSamplesGenes(datExpr0, verbose = 3);
gsg$allOK

sampleTree = hclust(dist(datExpr0), method = "average");
# Plot the sample tree: Open a graphic output window of size 12 by 9 inches
# The user should change the dimensions if the window is too large or too small.
sizeGrWindow(12,9)
#pdf(file = "Plots/sampleClustering.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)

Plot a line to show the cut
abline(h = 15, col = "red");
# Determine cluster under the line
clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
table(clust)
# clust 1 contains the samples we want to keep.
keepSamples = (clust==0)
datExpr = datExpr0[keepSamples, ]
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)

traitData = read.csv("ClinicalTraits.csv");
dim(traitData)
names(traitData)

#allTraits = traitData[, -c(31, 16)];
allTraits = traitData[, c(2, 3,4) ];
dim(allTraits)
names(allTraits)

femaleSamples = rownames(datExpr);
traitRows = match(femaleSamples, allTraits$Pig);
datTraits = allTraits[traitRows, -1];
rownames(datTraits) = allTraits[traitRows, 1];
collectGarbage()

# Re-cluster samples
sampleTree2 = hclust(dist(datExpr), 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.
plotDendroAndColors(sampleTree2, traitColors, groupLabels = names(datTraits), main = "Sample dendrogram and trait heatmap")


lnames = load(file = "FemaleLiver-01-dataInput.RData");
lnames
powers = c(c(1:10), seq(from = 12, to=20, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
# Plot the results:
sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;

# Scale-free topology fit index as a function of the soft-thresholding power
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.90,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")

net = blockwiseModules(datExpr, power = 6,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "femaleMouseTOM",
                       verbose = 3)

# open a graphics window
sizeGrWindow(12, 9)
# Convert labels to colors for plotting
mergedColors = labels2colors(net$colors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs;
geneTree = net$dendrograms[[1]];
save(MEs, moduleLabels, moduleColors, geneTree,
     file = "FemaleLiver-02-networkConstruction-auto.RData")

lnames = load(file = "FemaleLiver-01-dataInput.RData");
#The variable lnames contains the names of loaded variables.
lnames
# Load network data saved in the second part.
lnames = load(file = "FemaleLiver-02-networkConstruction-auto.RData");
lnames
enableWGCNAThreads()
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)

# Calculate topological overlap anew: this could be done more efficiently by saving the TOM
# calculated during module detection, but let us do it again here.
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);
# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
plotTOM = dissTOM^7;
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA;
# Call the plot function
sizeGrWindow(9,9)
TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")

 

Could please help me in fixing this error.

I tried looking through Search engines but it was not helpful.

 

Thank you

WGCNA TOMplot • 5.7k views
ADD COMMENT
0
Entering edit mode

How do you build the dendrogram? How do you create the plotTOM and the geneTree?

ADD REPLY
0
Entering edit mode

Hi,

Please  have a look at the updated question.

ADD REPLY
0
Entering edit mode

Have you checked the dimension of the data? I mean you seem to mix files from the tutorial and yours.

ADD REPLY
0
Entering edit mode

Yes, the function is certainly using my files and I did not mixed it up with the tutorial files. However, the names are same.


 

ADD REPLY
1
Entering edit mode
@abhishek-singh-4725
Last seen 10 months ago
France

Hi,

 

The reduction in number of genes from 12000 to 3500 did the trick for me. I still don't know the reason but will wait for the authors to reply on this.

 

ADD COMMENT
0
Entering edit mode

Abhishek, 

Did you ever get a response from the authors? Thanks!

 

 

 

 

ADD REPLY
0
Entering edit mode

Hello,

I also met this error when I doing this process. 
After searching for the solution of this problem, I found someone post a possible solution of this error. 
It told that the major problem of this error may come from the length difference of this three variables( plotTOM, geneTree and moduleColors.) (geneTree is less than the other two)

Back to the source of geneTree:

geneTree = net$dendrograms[[1]]

net = blockwiseModules(datExpr, power = 6, TOMType = "unsigned", minModuleSize = 30, reassignThreshold = 0, mergeCutHeight = 0.25, numericLabels = TRUE, pamRespectsDendro = FALSE, saveTOMs = TRUE, saveTOMFileBase = "femaleMouseTOM", verbose = 3)

The main reason of it comes from the default set of maxBlockSize in blockwisModules function is 5000. 
When our data size is over than the default size, this function may automatically divide your data into several parts.

(maxBlockSize: integer giving maximum block size for module detection. Ignored if blocks above is non-NULL. Otherwise, if the number of genes in datExpr exceeds maxBlockSize, genes will be pre-clustered into blocks whose size should not exceed maxBlockSize.https://www.rdocumentation.org/packages/WGCNA/versions/1.61/topics/blockwiseModules)

The possible solution that it suggested is to change the maxBlockSize value over the value of our own data.

Hope this will help.

ref site: http://www.biotrainee.com/thread-205-1-1.html

ADD REPLY
0
Entering edit mode

In addition, I just found that there is already some caution related to this problem in the WGCNA tutorials.

Here is the content in the tutorials:

"A word of caution 
for the readers who would like to adapt this code for their own data. The function blockwiseModules has many parameters, and in this example most of them are left at their default value. We have attempted to provide reasonable default values, but they may not be appropriate for the particular data set the reader wishes to analyze. We encourage the user to read the help file provided within the package in the R envi- ronment and experiment with tweaking the network construction and module detection parameters. The potential reward is, of course, better (biologically more relevant) results of the analysis.

A second word of caution concerning block size. 
In particular, the parameter maxBlockSize tells the function how large the largest block can be that the reader’s computer can handle. The default value is 5000 which is appropriate for most modern desktops. Note that if this code were to be used to analyze a data set with more than 5000 probes, the function blockwiseModules will split the data set into several blocks. This will break some of the plotting code below, that is executing the code will lead to errors. Readers wishing to analyze larger data sets need to do one of the following: • If the reader has access to a large workstation with more than 4 GB of memory, the parameter maxBlockSize can be increased. A 16GB workstation should handle up to 20000 probes; a 32GB workstation should handle perhaps 30000. A 4GB standard desktop or a laptop may handle up to 8000-10000 probes, depending on operating system and ihow much memory is in use by other running programs. • If a computer with large-enough memory is not available, the reader should follow Section 2.c, Dealing with large datasets, and adapt the code presented there for their needs. In general it is preferable to analyze a data set in one block if possible, although in Section 2.c we present a comparison of block-wise and single-block analysis that indicates that the results are very similar."

Ref: Tutorial for the WGCNA package for R: I. Network analysis of liver expression data in female mice 2.a Automatic network construction and module detection Peter Langfelder and Steve Horvath November 25, 2014

ADD REPLY

Login before adding your answer.

Traffic: 1009 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