Tree cut issue in WGCNA

library(WGCNA);
# The following setting is important, do not omit.
library(flashClust)
enableWGCNAThreads()
options(stringsAsFactors = FALSE);
#Read in the female liver data set
femData = read.csv("CTRL_WGCNA.txt",sep = "\t");

head(femData)
# Read in the male liver data set
maleData = read.csv("BPD_WGCNA.txt",sep = "\t");
# Take a quick look at what is in the data sets (caution, longish output):
dim(femData)
names(femData)
dim(maleData)
names(maleData)

#=====================================================================================
#
#  Code chunk 2
#
#=====================================================================================

# We work with two sets:
nSets = 2;
# For easier labeling of plots, create a vector holding descriptive names of the two sets.
setLabels = c("Female liver", "Male liver")
shortLabels = c("Female", "Male")
# Form multi-set expression data: columns starting from 9 contain actual expression data.
multiExpr = vector(mode = "list", length = nSets)

multiExpr[[1]] = list(data = as.data.frame(t(femData[-c(1)])));
names(multiExpr[[1]]$$data) = femData$$Gene;
rownames(multiExpr[[1]]$$data) = names(femData)[-c(1)]; multiExpr[[2]] = list(data = as.data.frame(t(maleData[-c(1)]))); names(multiExpr[[2]]$$data) = maleData$$Gene; rownames(multiExpr[[2]]$$data) = names(maleData)[-c(1)];
# Check that the data has the correct format for many functions operating on multiple sets:
exprSize = checkSets(multiExpr)

#=====================================================================================
#
#  Code chunk 3
#
#=====================================================================================

# Check that all genes and samples have sufficiently low numbers of missing values.
gsg = goodSamplesGenesMS(multiExpr, verbose = 3);
gsg$allOK #===================================================================================== # # Code chunk 4 # #===================================================================================== if (!gsg$$allOK) { # Print information about the removed genes: if (sum(!gsggoodGenes) > 0) printFlush(paste("Removing genes:", paste(names(multiExpr[[1]]data)[!gsggoodGenes], collapse = ", "))) for (set in 1:exprSizenSets) { if (sum(!gsggoodSamples[[set]])) printFlush(paste("In set", setLabels[set], "removing samples", paste(rownames(multiExpr[[set]]data)[!gsggoodSamples[[set]]], collapse = ", "))) # Remove the offending genes and samples multiExpr[[set]]data = multiExpr[[set]]data[gsggoodSamples[[set]], gsg$$goodGenes]; } # Update exprSize exprSize = checkSets(multiExpr) } #===================================================================================== # # Code chunk 5 # #===================================================================================== sampleTrees = list() for (set in 1:nSets) { sampleTrees[[set]] = hclust(dist(multiExpr[[set]]$data), method = "average")
}

#=====================================================================================
#
#  Code chunk 6
#
#=====================================================================================

pdf(file = "SampleClustering.pdf", width = 12, height = 12);
par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
xlab="", sub="", cex = 0.7);
dev.off();

#=====================================================================================
#
#  Code chunk 7
#
#=====================================================================================

# Choose the "base" cut height for the female data set
baseHeight = 80
# Adjust the cut height for the male data set for the number of samples
cutHeights = c(80, 100*exprSize$$nSamples[2]/exprSize$$nSamples[1]);
# Re-plot the dendrograms including the cut lines
pdf(file = "SampleClustering.pdf", width = 12, height = 12);
par(mfrow=c(2,1))
par(mar = c(0, 4, 2, 0))
for (set in 1:nSets)
{
plot(sampleTrees[[set]], main = paste("Sample clustering on all genes in", setLabels[set]),
xlab="", sub="", cex = 0.7);
abline(h=cutHeights[set], col = "red");
}
dev.off();

#=====================================================================================
#
#  Code chunk 8
#
#=====================================================================================

for (set in 1:nSets)
{
# Find clusters cut by the line
labels = cutreeStatic(sampleTrees[[set]], cutHeight = cutHeights[set])
# Keep the largest one (labeled by the number 1)
keep = (labels==1)
multiExpr[[set]]$$data = multiExpr[[set]]$$data[keep, ]
}
collectGarbage();
# Check the size of the leftover data
exprSize = checkSets(multiExpr)
exprSize

$nSets [1] 2$nGenes
[1] 8213

$nSamples [1] 0 0$structureOK
[1] TRUE


I think the issue is with the tree cut height because when i check the label all of them are zeroes i guess that is issue so far which is why im getting exprSize nsamples as 0 0instead of my actual sample size.

I have tried it with different cut height but so far no success may its a pretty simple issue but im not able to figure it out.

Data which im using here is the link

Any suggestion or help would be highly appreciated

1 Answer

It is likely that the cut heights set in the code (80 and 100*ratio of sample numbers) are too low. Look at the sample trees plotted into a pdf file (SampleClustering.pdf) and note the merging heights. You need to set the cut heights appropriately so that most samples go into a single cluster.

Peter

• "(SampleClustering.pdf) and note the merging heights. You need to set the cut heights appropriately so that most samples go into a single cluster." so my y axix upper limit is 120 .Should i set it to 120 as its the highest height .I will attach the fig which was generate for the clustering – krushnach Chandra Oct 17 at 19:16
• The issue is resolved now – krushnach Chandra Oct 18 at 8:07