Aligning short terminal exons

[MichaelHiller]

Dear Heng and team,

we are aligning polished IsoSeq reads against genomes and I noticed a few cases where the last exon of 13 bp was softclipped and not aligned. However, the exon perfectly aligns, if one creates a ~370 bp intron with GTG ... AAG consensus splice sites.

Blat correctly places the terminal 13 bp exon (second block), while minimap2 stops the alignment at the end of the first block in this image:

We are calling minimap2 even with
minimap2 --eqx -a -c -t $nThreats -ax splice:hq -uf --secondary=no -C5 -o ${P_outMinimap2}/ALL.CuP.aln.sam --junc-bed ${ref_annot} -cs long ${genome_fa} ${P_out_isoC}/ALL.CuP.fasta.gz

Is there a way to increase sensitivity to correctly align such terminal exons? I would be happy to spend a more to get the right alignment.

Thx a lot
Michael

[lh3]

Do you have the sequences? It might be possible to tune parameters to get the alignment.

[MichaelHiller]

Sure, I attach the read and the genomic context around that gene. Running this reduced example gives the same alignment.
read.fa.gz
genome.fa.gz

[lh3]

Thanks for the example. It is not possible to tune parameter to get the right alignment. I will keep this issue open and think more in future.

[MichaelHiller]

Alright, thanks a lot for looking into this.

Conceptually, HiSat2 (and other mappers?) make a list of downstream exon candidate positions that are inferred from mappings of other reads (in our case they are already given via --junc-bed ${ref_annot}) and then check if shorter terminal exons align to these candidates.
This could be a postprocessing step for reads that have unaligning terminal parts.

[nextgenusfs]

I wrote gapmm2 to address similar issue, it tries to re-align the terminal exons using edlib. I'm sure there are other ways to do this as well. https://github.com/nextgenusfs/gapmm2

[MichaelHiller]

Great. We will give it a try. Thanks for posting it