minimap2 instead of bowtie2?

[nick-youngblut]

Given that minimap2 is 2-3X faster than bowtie2 for mapping Illumina reads, I'd prefer to use minimap2. Have you done any testing to compare minimap2 versus bowtie2 in regards to the accuracy of CoPTR?

[tyjo]

I looked into minimap2 using the example data. It appears to produce fewer alignments than bowtie2. Here are the number of mapped reads from each aligner:

bowtie2

• ERR969281.bam: 6852
• ERR969282.bam: 6377
• ERR969283.bam: 12628
• ERR969285.bam: 10416
• ERR969286.bam: 10722
• ERR969428.bam: 25409
• ERR969429.bam: 57803
• ERR969430.bam: 78986

minimap2 (using the -x sr flag for short reads)

• ERR969281.bam: 4141
• ERR969282.bam: 3827
• ERR969283.bam: 7756
• ERR969285.bam: 6498
• ERR969286.bam: 6649
• ERR969428.bam: 24041
• ERR969429.bam: 54927
• ERR969430.bam: 73973

Note that CoPTR does filter reads after alignment. Even with the filtering step, bowtie2 produces more alignments.

Given the discrepancy, I would be hesitant to use minimap2 with CoPTR.

[nick-youngblut]

Good to know! Thanks for the info

[cdiener]

Hey, I think the discrepancies you see are due to minimap2 not returning secondary alignments by default. You get much better agreement with --secondary=yes (about 10% more mappings). For example:

minimap2 -ax sr --secondary=yes ref-db/db.fna fastq/ERR969281.fastq.gz | samtools view -c -F 4
[M::mm_idx_gen::0.334*0.99] collected minimizers
[M::mm_idx_gen::0.396*1.30] sorted minimizers
[M::main::0.396*1.30] loaded/built the index for 190 target sequence(s)
[M::mm_mapopt_update::0.396*1.30] mid_occ = 1000
[M::mm_idx_stat] kmer size: 21; skip: 11; is_hpc: 0; #seq: 190
[M::mm_idx_stat::0.415*1.28] distinct minimizers: 1070677 (99.05% are singletons); average occurrences: 1.015; average spacing: 6.006; total length: 6527538
[M::worker_pipeline::0.533*1.41] mapped 10818 sequences
[M::main] Version: 2.18-r1015
[M::main] CMD: minimap2 -ax sr --secondary=yes ref-db/db.fna fastq/ERR969281.fastq.gz
[M::main] Real time: 0.547 sec; CPU: 0.765 sec; Peak RSS: 0.057 GB
6049

And you can get more mappings when lowering k and w:

minimap2 -ax sr --secondary=yes -k 15 -w 10 ref-db/db.fna fastq/ERR969281.fastq.gz | samtools view -c -F 4
[M::mm_idx_gen::0.412*1.00] collected minimizers
[M::mm_idx_gen::0.518*1.40] sorted minimizers
[M::main::0.519*1.40] loaded/built the index for 190 target sequence(s)
[M::mm_mapopt_update::0.519*1.40] mid_occ = 1000
[M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 190
[M::mm_idx_stat::0.542*1.38] distinct minimizers: 1183441 (97.71% are singletons); average occurrences: 1.030; average spacing: 5.354; total length: 6527538
[M::worker_pipeline::0.676*1.55] mapped 10818 sequences
[M::main] Version: 2.18-r1015
[M::main] CMD: minimap2 -ax sr --secondary=yes -k 15 -w 10 ref-db/db.fna fastq/ERR969281.fastq.gz
[M::main] Real time: 0.695 sec; CPU: 1.066 sec; Peak RSS: 0.066 GB
7149

This is a bit more consistent with results from the literature where minimap2, bwa-mem and bowtie tend to perform similarly. However, AFAICT you can use any bam with coptr so I don't see an issue to default to bowtie which is more probably more sensitive and can be pretty memory efficient with the --mm flag.