Hello all , thanks for biobakery’s perfect workflow, I can easily calculate gene families abundance from metagenomic seq. But now I met a question. HUMANN provides a genefamilies.tsv to show the abundance of gene families and the contributions of species. But when i compare the gene number of MAG annotated by PROKKA and gene families with non-zero value from HUMANN in the sample sample, I surprisingly found that gene number is far less than gene families with non-zero value (eg Escherichia_coli MAG contains about 4k genes, but gene families has about 20K ).How could I explain this difference? Or is there any way to valuate the gene absence or presence for each species from shotgun metagenome data? Any valuable answer is welcomed.
This is one of the challenges that can arise when profiling genes from short reads (vs. assembling and annotating a genome). If three genes are present in the reference database, AACCGGTT + AACC + GGTT, it can be difficult to differentiate the presence of the former from the presence of the latter two (and vice versa), and so all tend to be reported, which inflates the number of detected genes. This can be particularly problematic for large pangenomes, which are more likely to contain variants of gene sequence families.
I find that when you group such genes into broader functions (e.g. EC terms), this behavior tends to smooth out. E.g. the three genes above will all group into the same broader function OR only the first (real?) gene will group, and the latter two (fragments?) will be discarded.
1 Like
Thanks for your advice! I’ll try it later