Dear community members,

I am currently working with the microbiome analysis of a plant using the RNASeq data. I have 3 replication with 2 treatment data. I have normalized count data of each microbe in the above data sets. I was wondering if I could take the count data and use edgeR to calculate the foldchange of the organism? Say, I may get staphylococcus aureus is 3 log2FC upregulated in treatment compared to control? Instead of gene, I would be saying particular organism is abundant or decreased abundance with respect to Fold change.

Please do clarify and suggest how it has to be carried out.

Hi Megha,

This is doable and I have done it in Alzheimer's disease context, wherein for 628 canonical genesets, pathways (KEGG, Biocarta, Wikipathways etc.). I did it in the following way:

1. Identify gene <--> gene-set/pathway association. For instance, MAPK signalling --> MAPK1, ERK1, MAPK12
2. Summarise average gene expression on to the gene-set/pathway level, thereby deriving pathway level expression
3. Apply edgeR/limma based differential analysis. You would need to test it with different p-value thresholds as pathway expression is quite robust compared to gene expression.

Ultimately you end up with differentially expressed pathways and you can have logFC between your 2 organisms. However, from point (3) these fold changes would be quite small, so bear that in mind.