I have an alignment of an assembled contig with the reference genome. The contigs are very long (PacBio + HiC assembly), and one "query" will generate multiple alignments with the "reference", as a consequence of structural variation, which breaks up the alignments.

I want to figure out if those alignments are in *the right order* relative to the contig they originate from. Next step would be to check if any of these alignments have deletions - leading to gaps, or change direction (due to inversions). I can't find any pysam attributes which seem do to the trick, right now, and I'm not sure if it's encoded in some obscure tag of the bam?

Suggestions?
Hi,

We (in our lab) had a similar problem, were we have reads with many supplementary alignments. Apparently, more people has the same problem as [you and me][1]

In order to get the supplementary tag as a function of the query (read), what I did was the following:

 1. Align the reads with minimap2 (I have nanopore data), and output in paf format
 2. Sort the data by column 1 and then by column 2 (sort -k 1,1 -k 2,2). In this case, you are sorting the alignments first by query name (you get all the alignments for a read one after each other) and then when you sort by column two you get the supplementary alignments one after each other **as a function of the query**


Let me know if this helps,


 


  [1]: https://github.com/lh3/minimap2/issues/386