Comparing gene abundances between metagenomes

Question

My workflow until now:

Find fragments of a marker gene in unassembled metagenomes > download and assemble metagenomes > recover the gene neighborhood / gene set of interest

Right now I have a rough estimate of how abundant these genes are by the 'depth' noted on the assembled contigs (I use MEGAHIT for assembly). I was wondering if there's a more thorough/proper way to do this. I would like to compare the abundance of specific genes between a) samples in the same study, and b) different studies. I imagine that the size of the individual metagenomes should be considered in both cases, but point b) might add additional difficulties such as different sequencing techniques. I would appreciate your insights.

Answer 1

I would avoid using assemblies to answer this question, as there's no guarantee that you will be able to assemble your genes of interest; you can however estimate their abundance even if they are relatively rare.

How I understand your question as being one of estimating the abundance of either some specific genes (e.g. butyrate metabolism genes) or all genes in a microbial community across multiple samples for comparative purposes. In other words, not 16S or marker gene analysis for the purposes of estimating organismal abundance, which is a rather different problem (though in that case I would still not use an assembly).

A more standard workflow is:

1. align metagenomic reads against some existing database of genes annotated appropriately.
2. estimate the number of reads aligning against some gene or ortholog in some way (using for example KEGG Orthology or similar).
3. use counts from (2) as input to some statistical procedure, possibly summarizing across functional categories.

Some examples of how this has been done are here, here, here. I am sure that there are more recent/relevant references but I haven't been following the field closely in the last few years.