Hello, I just have been struggling with a task.

I used mmseqs2 to perform sequence alignment form metagenomic contigs to reference gene database. The output of mmseqs2 was a blast+ table (format 6).

besides that, I mapped the metagenomic reads against the contigs to elaborate a bam file.

From these data, what is the most efficient way of performing a region abundance calculation (nucleotide regions defined on the blast+ output table) form a bam file?

I have been trying converting the blast+ table to gff, then to gtf and calculate the coverage per each region using htseq but there has been problems with the tables and htseq() outputs 0 coverage per each region.

This is the format of my gtf file (showing 3 first columns)

This is the format of my of my file (showing 3 first columns)

If you recommend me another way of doing this I'll glad to hear your recommendations

If it helps I built a Granges object in R with the Granges() function