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
