Hello,

I am using WGCNA to identify co-expressed gene networks from an RNAseq data set in which 24 mice were assigned to one of 4 experimental treatment groups. I have cleaned my data and there are 13861 genes that were used for the network construction using a signed network type, power = 18 (recommended from the network fit indices), deepSplit = 4, minModuleSize = 30, mergeCutHeight = 0.25. I have 19 modules ranging in size from 52 - 5486 genes. Now I would like to determine if there are any modules are related to experimental condition. Does anyone have a recommendation as to how I should compare modules across groups?

Thus far, I have run one-way ANOVAs on each module eigengene and corrected using FDR < 0.05 with the following loop. Note: datTraits[, 1] contains my grouping variable

Filename = net4_30_18 #This is the name of my network 

Unlisted = data.frame()

i = 1

while (i <= NCOL(Filename$MEs)){           
    Unlisted = rbind(Unlisted, unlist(summary(aov(Filename$MEs[ , i] ~ datTraits[ , 1]))))
  i = i+1  
}

names(Unlisted) = c("DF1", "DF2", "SS1", "SS2", "MS1", "MS2", "F1", "F2", "PR1", "PR2")
row.names(Unlisted) = names(Filename$MEs)

Cutoff = 0.05/NCOL(Filename$MEs)

summary(Unlisted$PR1 <= 0.1)
summary(Unlisted$PR1 <= 0.05)
summary(Unlisted$PR1 <= 0.005)
summary(Unlisted$PR1 <= Cutoff)

summary(p.adjust(Unlisted$PR1, method = 'fdr') <= 0.05)
summary(p.adjust(Unlisted$PR1, method = 'fdr') <= 0.1)

Unlisted[p.adjust(Unlisted$PR1, method = 'fdr') <= 0.05, ]

Does anyone have any recommendations on how they would approach this?

Thanks,

Mike

Well, this is really a question on how to compare 4 groups that has nothing to do with WGCNA. I think you did the WGCNA part well, now you want to test changes of a few (19, to be precise) variables between 4 treatment groups. There are multiple ways to test differences among 4 groups and you may want to consult a statistician (rather than the BioC mailing list).

Thank you, Peter. Confirmation that my preliminary steps are on the right track is reassuring.