Problem for understanding the number of modules in construction network based on WGCNA algorithm

Question

[modarzi]

Hi, for using Co-expression network, I use WGCNA algorithm. based on WGCNA tutorial, I use blockwiseModules() and detect my module.So, I would like to know how many module do I have? in other words, by which query or function I can find my answer?

I appreciate if anybody share his/her with me.

Best Regards,

Mohammad

Answer 1

[Jake Warner]

Hi,

blockwiseModules stores the module assignments in colors.

So to count the modules you can see how many levels are in the colors vector:

my_network <-  = blockwiseModules(datExpr, ...)
module_assignments <- my_network$colors

#count the modules
nlevels(factor(module_assignments))

[modarzi]

Dear Dr. Warner

Thanks for your comment. I have another problem. As you know for preparing output file for gene ontology analysis in WGCNA, I need GeneAnnotation.csv file for my Gene Expression profile. My Gene Expression data set is RNA-seq data from The Genome Cancer Atlas (TCGA). For that purpose I downloaded "gencode.v22.genes.csv" but the attributes name of "gencode.v22.genes" is not similar to GeneAnnotaion.csv file in Tutorials for the WGCNA package.

I don't know how should I handle Annotation file for my dataset?

And my second problem is that below code is based on MicroArray and my data set is RNA-seq.

# Read in the probe annotation
GeneAnnotation=read.csv(file="GeneAnnotation.csv")
# Match probes in the data set to those of the annotation file
probes = names(datExprFemale)
probes2annot = match(probes,GeneAnnotation$substanceBXH)
# data frame with gene significances (cor with the traits)
datGS.Traits=data.frame(cor(datExprFemale,datTraits,use="p"))
names(datGS.Traits)=paste("cor",names(datGS.Traits),sep=".")
datOutput=data.frame(ProbeID=names(datExprFemale),
                     GeneAnnotation[probes2annot,],moduleColorsFemale,datKME,datGS.Traits)
# save the results in a comma delimited file
write.table(datOutput,"FemaleMouseResults.csv",row.names=F,sep=",")

How can I customize that code for RNA-seq?

I appreciate if you share your comment with me.

[Jake Warner]

Hi Modzari,

This is more of an R problem and less of a WGCNA or RNAseq/array problem. What the code you pasted does, is import an annotation file which is a series of columns with gene name, probe name, genomic coordinates etc, and associate those annotations with the results of the correlation test performed as cor(datExprFemale,datTraits,use="p"). To adapt this you would 1. Run the correlation test on your own trait data and 2. associate the results with your own annotation file. To associate the results you just need to use a key value (gene ID, probe ID, transcript ID) that exists in both your datExpr dataframe and your annotation file. Then you use any number of functions to put them together including merge() from base, or join(), from dplyer to combine the tables. The tutorial uses match() followed by a square bracket subset but I find that clunky. I prefer to use the join family of functions from dplyr.

I'm shooting in the dark without running the whole tutorial myself but a simplified version of the above would look like this:

#read anntoations
GeneAnnotation=read.csv(file="GeneAnnotation.csv")

#get the gene names. Looks like datExprFemale is wide??
gene = names(datExprFemale)

# run the correlation
datGS.Traits=data.frame(cor(datExprFemale,datTraits,use="p"))

#flip the dataframe if it is indeed wide
datGS.Traits = t(datGS.Traits)

#add the genes
datGS.Traits$gene = gene

#make sure you have a 'gene' in your annotation file for the join to work
datGS.Traits.annotated = dplyr::full_join(datExprFemale, datGS.Traits, by='gene')

[modarzi]

Dear Dr Warner

Thanks for comment and code. As you know, I use RNA-seq data set as an input in WGCNA . Now, I want to import my network in Cytoscape for visualization. based on WGCNA tutorial, for that purpose I have to run below code that is based on MicroArray data set:

# select modules modules = c("blue","brown") 
# Select module probes 
inModule=is.finite(match(moduleColorsFemale,modules)) 
modProbes=probes[inModule] 
match1=match(modProbes,GeneAnnotation$substanceBXH) 
modGenes=GeneAnnotation$gene_symbol[match1] 
# Select the corresponding Topological Overlap modTOM = TOM[inModule, inModule] 
dimnames(modTOM) = list(modProbes, modProbes) 
# Export the network into edge and node list files for Cytoscape 
cyt = exportNetworkToCytoscape(modTOM, 
edgeFile=paste("CytoEdge",paste(modules,collapse="-"),".txt",sep=""), nodeFile=paste("CytoNode",paste(modules,collapse="-"),".txt",sep=""), 
weighted = TRUE, threshold = 0.02,nodeNames=modProbes, 
altNodeNames = modGenes, nodeAttr = moduleColorsFemale[inModule])

How can I customize that code for RNA-seq?

I appericite if you share your comment with me.

Best Regards,

Mohammad