Dear @misialq,
Functional gene classification of contigs using qiime moshpit eggnog-diamond-search works well, running emapper.py under the hood.
However, the output produced by the moshpit implementation only includes the , emapper.seed_orthologs tsv file.
While this is sufficient input for a subsequent annotation using qiime moshpit eggnog-annotate, this does not allow calculating read abundance for the functional genes.
In my projects, which include paired MGX and MTX datasets, it is most relevant to be able quantify the (relative) expression of functional genes, rather than just detecting them.
What I imagine as a straightforward way is to let qiime2 output the .gff file that is created by emapper.py, additional to the emapper.seed_orthologs output. Supposedly the .gff contains the information of the location of predicted genes on the contigs.
Converting the .gff file to a .bed file using bedtools *gff2bed* and then intersecting that with the .bam (sorted?) alignment file produced by qiime moshpit map-reads-to-contigs using bedtools *coverage* should produce the read count per gene. Could extracting gene length from emapper.seed_orthologs file even enable normalization to RPKM?
Digging through my slurm stdout file, I observed that emapper.py defaults to producing the .gff file, but that it has been writing to /temp and that, in contrast to the emapper.seed_orthologs file, the .gff is discarded as the q2 action completes.
I'm not a software engineer, but I can imagine it is a fairly simple adjustment to let the .gff file be written to qza output.
Thanks in advance for your thoughts.
Best,
Pieter
