Is there an alignment program that allows you to sort results by coverage (amount of reference sequence covered)?

Question

[CC]

I am wanting to align reads to a database of reference sequences. This database contains thousands of reference sequences. I would like to only see the reference sequences that the reads cover a certain percentage, eg. only see the reference sequence if reads cover >50% of it.

Currently, I am doing my alignment with BWA MEM, and the output is really hard to analyse, because I'll have hundreds of reference sequences with reads aligning to them, but 80% of these only have reads aligning to one little spot of the reference sequence, which doesn't really make it a match.

Is there an alignment program that will allow me align reads to multiple sequences, and then sort results by the amount of the reference sequence that is covered? Or is there a way that I can this information out of a SAM or BAM file?

I'm sorry if this question has been answered before but I couldn't find a post that answered this for me.

Thank you for your help.

[lieven.sterck]

that is information you could retrieve from a blast search but if you are working with really big sets it might not be the most efficient approach. MegaBlast might already be better suited.

Answer 1

[Pierre Lindenbaum]

Using my custom sam comparator samcustomsortjdk : http://lindenb.github.io/jvarkit/SamCustomSortJdk.html and using function Cigar.getReferenceLength . https://github.com/samtools/htsjdk/blob/master/src/main/java/htsjdk/samtools/Cigar.java#L76

 java -jar dist/samcustomsortjdk.jar  --body -e 'private int score(final SAMRecord R) { if(R.getReadUnmappedFlag() || R.getCigar()==null) return 0; return R.getCigar().getReferenceLength();} @Override public int compare(SAMRecord a,SAMRecord b) { return Integer.compare(score(a),score(b));}' input.bam