This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Infer order of supplementary alignments [pysam]

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?

pysam alignment
• 2,607 views • link
updated by IP IP • written by WouterDeCoster WouterDeCoster

1 answer

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

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,

• 0 views • link
updated by WouterDeCoster WouterDeCoster • written by IP IP

Thanks a lot, that sounds very useful. I haven't looked at PAF alignments yet, only in bam format. What I ended up doing was to write a Python script and parse the CIGAR strings to get the order of the supplementary alignments. A preliminary script can be found here.

I'll definitely take a look at your approach the next time I need to do a similar task.

• 0 views • link
written by WouterDeCoster WouterDeCoster

Looks cool! I just had a quick look at the code and it seems that it is a better approach for high throughput processing of the nanopore complex reads. Thanks a lot for sharing!

• 0 views • link
written by IP IP

Comparing the results with your method would definitely be the next step, as I'm not 100% confident I did not mess up somewhere!

• 0 views • link
written by WouterDeCoster WouterDeCoster

Log in to answer this question.

More posts like this