Forcing continuous alignments

[rlorigro]

Hi,

I am using minimap2 to map long ONT reads (up to 400Kbp), and I am finding that the supplementary chaining is producing large rearrangements in the read sequence, flipping strands, and jumping across regions. I have tried using the -m parameter to limit chaining, and several other parameters, but none of them seem to force the primary alignment to be linear and unidirectional.

I should emphasize that the reference I am aligning to is the same individual that the reads are derived from, so I expect there to be no SVs. Here is one example:

Read       Length Start  Stop      Region RefLength  Start      Stop
2626312    448022 112    922    +  chr2   242696747  130488450  130489279
2626312    448022 7047   7861   +  chr2   242696747  130489870  130490701
2626312    448022 41024  447988 +  chr2   242696747  130533922  130950797
2626312    448022 3105   3985   +  chr2   242696747  130885322  130886221
2626312    448022 12041  40959  -  chr2   242696747  131672947  131702229
2626312    448022 5227   5854   -  chr2   242696747  131710435  131711099
2626312    448022 2039   2944   -  chr2   242696747  131707727  131708642

Is there a parameter (or combination of parameters) that would help with this? I think what is needed is to allow the primary alignment to extend further.

Thanks

[rlorigro]

I should also add that this is a repetitive region, and that I have been using winnowmap to accommodate for that.

[lh3]

What is the winnowmap2 output?

[rlorigro]

The example I posted was using winnowmap already. In that case, I aligned HG002 reads to CHM13 v1.0, but in future runs I aligned HG002 reads to the HiFiasm HG002 assembly (~190x CCS coverage from here) and observed a similar problem, in other reads. More detailed information:

Command:

meryl count k=15 output merylDB /home/ryan/data/hifi/human/hg002/HG002.mat.fa
meryl print greater-than distinct=0.9998 merylDB > repetitive_k15.txt

winnowmap \
-W /home/ryan/data/hifi/human/hg002/repetitive_k15.txt \
-x map-ont \
-t 7 \
--secondary=no \
/home/ryan/data/hifi/human/hg002/HG002.mat.fa \
/home/ryan/data/nanopore/human/test/overlap_guppy_360_tangle_subset_1/HG002-Guppy-3.6.0-UL-run3-subset1.fasta \
> HG002-Guppy-3.6.0-UL-run3-subset1_VS_HG002_mat_primary.paf

stderr:

[M::mm_idx_gen::0.035*0.04] reading downweighted kmers
[M::mm_idx_gen::0.065*0.24] collected downweighted kmers, no. of kmers read=66775
[M::mm_idx_gen::0.065*0.24] saved the kmers in a bloom filter: hash functions=2 and size=960072 
[M::mm_idx_gen::98.885*1.08] collected minimizers
[M::mm_idx_gen::100.291*1.16] sorted minimizers
[M::main::100.291*1.16] loaded/built the index for 546 target sequence(s)
[M::main::100.291*1.16] running winnowmap in SV-aware mode
[M::main::100.291*1.16] stage1-specific parameters minP:1000, incP:4.00, maxP:16000, sample:1000, min-qlen:10000, min-qcov:0.5, min-mapq:5, mid-occ:5000
[M::main::100.291*1.16] stage2-specific parameters s2_bw:2000, s2_zdropinv:25
[M::mm_idx_stat] kmer size: 15; skip: 50; is_hpc: 0; #seq: 546
[M::mm_idx_stat::100.479*1.16] distinct minimizers: 23589668 (41.44% are singletons); average occurrences: 5.194; average spacing: 25.013
[M::worker_pipeline::111.965*1.78] mapped 2273 sequences
[M::main] Version: 2.0, pthreads=7, omp_threads=3
[M::main] CMD: winnowmap -W /home/ryan/data/hifi/human/hg002/repetitive_k15.txt -x map-ont -t 7 --secondary=no /home/ryan/data/hifi/human/hg002/HG002.mat.fa /home/ryan/data/nanopore/human/test/overlap_guppy_360_tangle_subset_1/HG002-Guppy-3.6.0-UL-run3-subset1.fasta
[M::main] Real time: 112.009 sec; CPU: 199.570 sec; Peak RSS: 3.884 GB

And some weird alignments I observed:

Jumping between contigs (not bridging contig termini):

736840  41280   2057    39942   +   h2tg000102l 36105583    464829  502652  7336    38234   60  tp:A:P  cm:i:541    s1:i:7250   s2:i:802    rl:i:0
736840  41280   75  968 +   h2tg000005l 111658231   629986  630886  327 910 60  tp:A:P  cm:i:27 s1:i:324    s2:i:0  rl:i:0

Reversing strand on the same contig:

207028  74088   33102   73900   +   h2tg000102l 36105583    3378    44916   7102    41787   60  tp:A:P  cm:i:510    s1:i:6921   s2:i:0  rl:i:0
207028  74088   70  36902   -   h2tg000102l 36105583    43512   80920   5749    37648   60  tp:A:P  cm:i:415    s1:i:5587   s2:i:194    rl:i:0

Jumping intra + inter contig

2027405 47862   34  37932   -   h2tg000005l 111658231   106444  145003  6510    38787   60  tp:A:P  cm:i:481    s1:i:6327   s2:i:968    rl:i:0
2027405 47862   39262   39765   -   h2tg000005l 111658231   712065  712577  90  519 56  tp:A:P  cm:i:6  s1:i:85 s2:i:0  rl:i:0
2027405 47862   41040   47808   -   h2tg000102l 36105583    146693  153675  614 7023    60  tp:A:P  cm:i:43 s1:i:561    s2:i:0  rl:i:0

[rlorigro]

There could still be SVs because the reads are not binned, and this is the maternal haplotype from HiFiasm, but in general I would just like to forbid chaining with large gaps and favor extending the primary alignment.

[lh3]

If you can give me just these reads, I can have a look. Nonetheless, I suspect those are the best you can get.

I should emphasize that the reference I am aligning to is the same individual that the reads are derived from, so I expect there to be no SVs.

Assembly can be wrong. Reads can be chimeric both locally and globally.

[rlorigro]

OK, here are the 3 reads from above: https://rlorigro-public-files.s3-us-west-1.amazonaws.com/reads/guppy_360/misaligned_reads_guppy_3.6.0_HG002.fasta

There are many more examples like these so I would be surprised if they are all chimeras. These reads are also dumped from Shasta, which is usually done after some filtering for chimerism. I think the issue is that there are repeats in this region, and the primary/supplementary segments are being split amongst the repeats.

Also I noticed that winnowmap is running with "sv-aware mode", but I couldn't find more info about this. Does this affect chaining and supplementaries? Maybe a question for @cjain7

[lh3]

I guess --sv-aware implements the algorithm in the winnowmap2 preprint.

[rlorigro]

Ok i will take a closer look at the paper. I just tried rerunning with --sv-off and it appears to have reduced these rearrangements.

[rlorigro]

With sv-off I do still get some strange ones though, unfortunately:

2497450 60236   146 29307   -   h2tg000005l 111658231   848222  878809  3055    30728   46  tp:A:P  cm:i:219    s1:i:2764   s2:i:2331   rl:i:636
2497450 60236   12923   60193   -   h2tg000102l 36105583    312140  361719  8282    49837   60  tp:A:P  cm:i:616    s1:i:7919   s2:i:543    rl:i:636

After looking some more at the paper, I am wondering if disabling sv-aware mode will defeat the purpose of using winnowmap over minimap.

[ekg]

@rlorigro wfmash -m -N -p 90 will generate contiguous approximate mappings of nanopore reads. You can remove -m to get a base-level alignment obtained by global pairwise alignment of each mapping result. wfmash lacks mapping quality, but it could provide you a baseline here. We're just beginning tests on non-assembly inputs, so please let me know if you try it and it works for you.

[rlorigro]

Thanks @ekg, this is interesting. As it turns out, I've been digging into these issues more and have found that @lh3 was right about many of them being chimeric. A majority of the chimeras were palindromes, where they reverse direction about halfway through the sequence, and a smaller portion of them were conventional chimeras, where a chunk of the read appears to have been a concatenation of totally unrelated sequence. For the palindromes, it seems that nanopore has issues with pulling in the complementary strand directly after the first strand finishes translocating, and the signal segmenter doesn't catch the transition, so both strands of the signal are sent to the basecaller as a single read.

There are occasionally still some unexplained artifacts, where supplementaries are almost entirely contained inside other alignments, but fortunately the mapQ score for these is very low or 0 in the cases I have observed so far.

So overall I would say the mappings were not issues with winnowmap (if SV mode is turned off), though I think it should still be possible to force more strict chaining.

[rlorigro]

To be more specific about the remaining non-chimeric artifacts I mentioned, here are 2 examples:

14617   59497   29  59436   +   h2tg000052l 87420292    61492027    61545226    18078   59881   60  tp:A:P  cm:i:1359   s1:i:17961  s2:i:255    rl:i:1688
14617   59497   38972   44839   +   h2tg000043l 49469036    33482998    33488884    1103    5906    2   tp:A:P  cm:i:83 s1:i:1097   s2:i:1085   rl:i:1688
14617   59497   36265   41733   -   h2tg000036l 139977256   27376259    27381830    75  5589    14  tp:A:P  cm:i:5  s1:i:50 s2:i:0  rl:i:1688

1181401 72122   105 72059   -   h2tg000052l 87420292    61507216    61573955    10694   73316   60  tp:A:P  cm:i:775    s1:i:10400  s2:i:40 rl:i:1116
1181401 72122   43191   48674   -   h2tg000052l 87420292    43644377    43649923    488 5563    0   tp:A:P  cm:i:33 s1:i:472    s2:i:472    rl:i:1116
1181401 72122   40607   46177   -   h2tg000077l 81620956    67987710    67993283    60  5590    13  tp:A:P  cm:i:4  s1:i:51 s2:i:0  rl:i:1116

They aren't overlapping by more than 50% in the reference (as specified by the -M parameter), but in the query coordinate space, some of these mappings are contained entirely inside other mappings. The final supplementary alignment of read 1181401 has a surprisingly high mapQ of 13 despite this.