Skip to content

High-MAPQ alignments are sometimes marked as supplementary to low-MAPQ primary alignments #1350

Description

@sjl

We are aligning ONT data with minimap2. Some of our reads have chimeric alignments, with one of the alignments having a high MAPQ (e.g. 60) and the other having a very low MAPQ (e.g. 1). We've noticed that sometimes minimap2 will choose the alignment with the low MAPQ as the primary alignment, which is surprising to us. We expected that the primary alignment would be the one with the higher score (if it's a tie then it'd be arbitrary, of course).

The SAM spec says the choice is arbitrary. I've tried to figure out if minimap2 intends to pick the alignment with the higher score (and thus this is a bug in the minimap2 code), or if it's intentionally also just arbitrary (in which case we just have to live with it). There are a couple of issues in the issue tracker where the responses imply that it should be marking the better alignment as primary, but it's not 100% clear:

Can you clarify whether the classification of primary vs supplementary should be determined by the MAPQ score or if the choice is intentionally arbitrary?

I've attached example.fastq.gz with an example read that we see this happening for. I'm using the latest git head of Minimap2:

$ git show --oneline -s
e5066c7 (HEAD -> master, origin/master, origin/HEAD) added sim2bed to convert simulated read names to BED

$ make -j6
cc -c -g -Wall -O2 -Wc++-compat  -DHAVE_KALLOC  main.c -o main.o
cc -c -g -Wall -O2 -Wc++-compat  -DHAVE_KALLOC  kthread.c -o kthread.o
...rest of the cc calls elided
ar -csru libminimap2.a kthread.o kalloc.o misc.o bseq.o sketch.o sdust.o options.o index.o lchain.o align.o hit.o seed.o jump.o map.o format.o pe.o esterr.o splitidx.o ksw2_ll_sse.o ksw2_extz2_sse41.o ksw2_extd2_sse41.o ksw2_exts2_sse41.o ksw2_extz2_sse2.o ksw2_extd2_sse2.o ksw2_exts2_sse2.o ksw2_dispatch.o
ar: `u' modifier ignored since `D' is the default (see `U')
cc -g -Wall -O2 -Wc++-compat  main.o -o minimap2 -L. -lminimap2 -lm -lz -lpthread

Align the FASTQ to a pre-built index for hg38_no_alt:

./minimap2 -a -y -x map-ont -t 4 --secondary=no --eqx -Y ~/scratch/hg38-minimap-index/index.mmi example.fastq | samtools sort -o result.bam
[M::main::4.781*1.00] loaded/built the index for 195 target sequence(s)
[M::mm_mapopt_update::5.709*1.00] mid_occ = 694
[M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 195
[M::mm_idx_stat::6.251*1.00] distinct minimizers: 100167746 (38.80% are singletons); average occurrences: 5.519; average spacing: 5.607; total length: 3099922541
[M::worker_pipeline::6.258*1.00] mapped 1 sequences
[M::main] Version: 2.30-r1299-dirty
[M::main] CMD: ./minimap2 -a -y -x map-ont -t 4 --secondary=no --eqx -Y /home/sjl/scratch/hg38-minimap-index/index.mmi example.fastq
[M::main] Real time: 6.608 sec; CPU: 6.606 sec; Peak RSS: 7.714 GB

Looking at the results, the alignment with MAPQ 1 is marked as primary, and the alignment with MAPQ 60 is the supplementary alignment:

$ samtools view result.bam | awk 'BEGIN { print "id, flag, supp, mapq" } { id=$1; flag=$2; mapq=$5; supp=and(flag, 0x800) > 0; printf("%s, %x, %x, %d\n", id, flag, supp, mapq); }' | column -s, -t
id                                     flag   supp   mapq
7cfbe678-7810-4d84-bacf-00e0636b793d   800    1      60
7cfbe678-7810-4d84-bacf-00e0636b793d   10     0      1

Thanks!

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions