I have the following problem with
WGCNA - http://labs.genetics.ucla.edu/horvath/htdocs/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/
Working on Section 1.6, Export of networks to external software (Cytoscape)
I'm currently trying to perform WGCNA on a set of genes and I'm having trouble getting the top x hub genes for each module. I am trying to export a network to Cytoscape and used the same method for getting the top x hub genes as outlined for exporting to VisANT.
# Select modules (only interested in one for now)
modules = c("greenyellow")
# Select module probes
probes = names(datExpr)
inModule = is.finite(match(bwModuleColors, modules))
modProbes = probes[inModule]
modGenes = annot$gene_symbols[match(modProbes, annot$geneID)]
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule]
dimnames(modTOM) = list(modProbes, modProbes)
# Restrict the network to the top 30 genes
nTop = 30
IMConn = softConnectivity(datExpr[, modProbes]
top = (order(-IMConn) <= nTop)
# Export the network into a fomat that Cytoscape can read
cyt = exportNetworkToCytoscape(modTOM[top, top],
edgeFile = paste("CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep = ""),
nodeFile = paste("CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep = ""),
weight = TRUE,
threshold = 0.02,
nodeNames = modProbes,
altNodeNames = modGenes,
nodeAttr = bwModuleColors[inModule])
I've written a short loop to count the number of connections to each gene and it works as expected, but the xth gene consistently has zero connections. Let's say that x is 30. If I increase the cutoff to 31 hub genes, the 30th gene now shows connections to the others in the network, but the 31st gene shows nothing. In addition, this change increases AND decreases some of the number of connections to other genes in the network. This really bothers me, because connections should only be added, since the network is getting bigger by one gene, and the changes should be accounted for by the 30th gene, but this is not the case for the output.
# Split the cytoscape file into two parts: edge and node
node <- cyt$nodeData
edge <- cyt$edgeData
# The limit covers all of the connections in the edge file by determining the length of the column ‘fromNode’
limit <- length(edge$fromNode)
# Create an empty list to store the counts for each gene
counts = list()
# Loop for the genes going from 1 to the number of genes specified for the network, ‘nTop’
for (i in 1:nTop) {
# Reset the count for each new gene and specify the names of the gene of interest and the matching genes
name = node$nodeName[[i]]
count = 0
# Nested loop that searches for matches to the gene in question in both the ‘fromNode’ and ‘toNode’columns, and adds one to the count for each match.
for (j in 1:limit) {
matchName1 = edge$fromNode[[j]]
matchName2 = edge$toNode[[j]]
if (name == matchName1 || name == matchName2)
{count = count + 1}
}
# Create a string for the attribute in the correct format
attribute <- paste(name, "=", count)
# Adds the count to the list
counts <- c(counts, attribute)
}
# End of loop
The loop seems to be working as expected, so I'm thinking that the problem is with the network construction. I'm currently referring back to what I know about linear algebra, matrices and topology to try to see if the problem is the way they're being sorted or something like that, but it might just be the way that the exportNetworkToCytoscape() function works.
