0
$\begingroup$
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(!gsg$goodGenes) > 0)
    printFlush(paste("Removing genes:", paste(names(multiExpr[[1]]$data)[!gsg$goodGenes], 
                                              collapse = ", ")))
  for (set in 1:exprSize$nSets)
  {
    if (sum(!gsg$goodSamples[[set]]))
      printFlush(paste("In set", setLabels[set], "removing samples",
                       paste(rownames(multiExpr[[set]]$data)[!gsg$goodSamples[[set]]], collapse = ", ")))
    # Remove the offending genes and samples
    multiExpr[[set]]$data = multiExpr[[set]]$data[gsg$goodSamples[[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

$\endgroup$
2
$\begingroup$

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

$\endgroup$
  • $\begingroup$ "(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 $\endgroup$ – krushnach Chandra Oct 17 at 19:16
  • $\begingroup$ The issue is resolved now $\endgroup$ – krushnach Chandra Oct 18 at 8:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.