SureMap

SureMap is a versatile, error-tolerant, high-sensitivity read mapper. It aligns DNA sequencing reads against a reference genome and, unlike most mappers, is built to handle "difficult" reads — the ones that need many edit operations before they match anywhere, either because the sequencer was noisy or because the sampled genome genuinely diverges from the reference in that region.

Three algorithms share one index:

Algorithm	Intended input	Idea
short	reads up to ~200 bp	Search the whole read within a fixed edit budget and report every location within that budget.
hard	reads other mappers leave unmapped	Anchor with short seeds first, then re-search the read around the anchor with a much larger error tolerance.
long	PacBio / Nanopore reads	Anchor with many short seeds, then stitch a fragment-wise approximate global alignment across the whole read.

SureMap is described in:

MohammadJavad Rezaei Seraji and Seyed Abolfazl Motahari. SureMap: Versatile, Error Tolerant, and High Sensitive Read Mapper. bioRxiv 173740 (2017). https://doi.org/10.1101/173740

---

Building

Aligner — needs only a C++11 compiler and pthreads:

make aligner

Indexer — additionally needs OpenMP, GNU libstdc++ parallel mode (<parallel/algorithm>, so GCC rather than Clang), and libdivsufsort built as static 64-bit libraries. Set the following in libdivsufsort's top-level CMakeLists.txt before building it:

option(BUILD_SHARED_LIBS "Set to OFF to build static libraries" OFF)
option(BUILD_DIVSUFSORT64 "Build libdivsufsort64" ON)

Then:

make indexer

make alone builds both. The suffix array is constructed with pSAscan:

Juha Kärkkäinen, Dominik Kempa, Simon J. Puglisi. Parallel External Memory Suffix Sorting. CPM 2015, pp. 329–342.

---

Usage

1. Index the reference

./SureMap-Indexer.o ref.fa /path/to/index/

This writes four files next to each other, all sharing the prefix /path/to/index/ref.fa:

File	Contents
.fm	FM-index of the reference
.rev.fm	FM-index of the reversed reference (used to extend a match rightwards)
.rinfo	the reference itself, packed at 2 bits per base
.cinfo	chromosome names and their offsets in the concatenated reference

Indexing the human genome takes a while and needs roughly 3 GB of working memory plus scratch disk space next to the current directory.

2. Align reads

./SureMap-Aligner.o <algorithm> [options] <ref.fa> <reads.fastq>

<ref.fa> is the prefix the index files were written under, i.e. exactly the path you passed to the indexer plus the reference filename. Options must come before the two positional arguments.

# short reads, 10 threads, expected error rate 10%, fast mode
./SureMap-Aligner.o short -t 10 -e 0.1 -m fast -o output.sam /path/to/index/hg19.fa illumina.fastq

# hard-to-map reads: up to 4 edits including indels, report at most 10 locations
./SureMap-Aligner.o hard -v 4 -k 10 -g -m sensitive -o output.sam /path/to/index/hg19.fa bwa-unmapped.fastq

# long reads, 500 bp fragments in the approximate global alignment
./SureMap-Aligner.o long -t 30 -L 500 -o output.sam /path/to/index/hg19.fa pacbio.fastq

Options

Option	Meaning	Default
-t <int>	worker threads (capped at 40)	1
-o <file>	output SAM file	report.sam
-k <int>	maximum alignments reported per read	1
-a	report every location found	off
-u	report only reads that map to exactly one location	off
-b	report the best alignment only	off
-g	allow insertions and deletions, not just mismatches	off
-e <float>	expected error rate; the edit budget becomes e × read length	0.05
-v <int>	absolute edit budget per read; overrides -e	—
-m <mode>	fast, normal, sensitive, very-sensitive	normal
-L <int>	fragment length of the long-read global alignment (long only)	500

-e, -v, -g, -k, -a, -u and -b apply to short and hard; long derives its own budget from its fixed 30% assumed error rate.

Output is SAM. Reads that could not be placed are emitted with flag 4. MAPQ is always reported as 255 — SureMap reports a list of candidate locations rather than a posterior over them.

---

How it works

Divide and merge. For a read R and edit budget d, split R into halves R₁R₂. Any alignment of R within d edits must align one half within ⌊d/2⌋ edits, because the edits distribute across the two halves. So SureMap maps each half within ⌊d/2⌋, then extends each surviving half across its sibling. Applied recursively the budget halves at every level, and at the leaves (d = 0) the problem is plain exact matching against the FM-index.

Extending left and right would report the same alignment twice, so left extensions are additionally required to cost strictly more than ⌊d/2⌋ — everything cheaper has already been found by the mirrored right extension.

Extension. Extending a matched half means walking the FM-index one character at a time while filling a banded edit-distance table, one new row per character. A row is dropped as soon as none of its cells can still reach the budget. Crucially, one BWT interval stands for every genomic position that shares the matched prefix, so a single DP row covers all of them at once, and the traversal is a tree over the extending characters rather than one pass per genomic position.

Row memoisation. The same row reaches the same point of the search over and over across that tree. Rows are packed to ⌈log₂(d+2)⌉ bits per cell so they fit in six 64-bit words, and each transition is memoised in a hash table keyed on the packed row plus the identity of the DP step (which node of the recursion tree, which row of its table, which character is being appended). Each worker owns its table and the table is scoped to a single read, so it needs no locking and a read is aligned identically no matter which thread picks it up.

Long reads. Short seeds are drawn across the read and searched with the machinery above. Once three seeds land on the genome with roughly the same spacing they have on the read, the read is anchored. The read is then cut into -L-sized fragments that are aligned one after another, each starting where the previous one ended; a second pass re-aligns a window straddling every fragment boundary to repair the drift, which is what lets SureMap emit a real CIGAR for the whole read instead of clipping it.

---

Benchmarks

The numbers below are from the paper, not re-measured in this repository.

Short reads — sensitivity and speed

10 M simulated reads of length 50 from hg19, mismatch errors only, 20 threads.

Read error rate	Aligner	Time (s)	Reads aligned
1 %	Bowtie -v 3	500	99.83 %
	SOAP2 -v 3	500	98.65 %
	SureMap -v 3	700	99.83 %
5 %	Bowtie -v 3	480	76.9 %
	SOAP2 -v 3	210	55.65 %
	SureMap -v 3	600	76.9 %
	SureMap -v 5	7 500	96.3 %
10 %	SureMap -v 5	10 100	63 %
	SureMap -v 7	16 500	88 %

Up to 10 locations are reported per read. Peak memory: Bowtie 2.4 GB, SOAP2 5.4 GB, SureMap 8.7 GB.

The trade-off is explicit: at a budget other mappers also support, SureMap matches them and costs ~40 % more time. The reason to reach for it is the rows they cannot produce at all — 96 % aligned at a 5 % error rate, 88 % at 10 %.

Short reads — uniquely aligned fraction

Read error rate	Bowtie -v 3	SureMap -v 3	SureMap -v 5	SureMap -v 7
1 %	81 %	81 %	—	—
5 %	65 %	65 %	77 %	—
10 %	—	—	52 %	70 %

Downstream analyses routinely discard multi-mapping reads, so the uniquely-aligned fraction sets the coverage they actually get to work with.

PacBio

25 249 reads, 180 Mbp total, mapped to hg19 on a 32-core machine (30 threads). v is the variation rate between a read and the region it aligned to; each column is the percentage of all reads aligned at or below that rate. Coverage is the total length of reads aligned with v ≤ 0.3.

Aligner	Option	Time	v ≤ 0.3	v ≤ 0.25	v ≤ 0.2	v ≤ 0.15	Coverage (Mbp)
SureMap	ℓ=200	95 s	85.5	79.8	68.2	44.9	130.0
SureMap	ℓ=300	107 s	88.0	82.2	70.5	46.3	138.8
SureMap	ℓ=650	162 s	90.8	84.9	72.3	47.2	150.7
SureMap	ℓ=1100	251 s	91.4	85.2	72.5	45.1	154.1
SureMap	ℓ=1600	365 s	92.0	85.6	72.5	47.2	156.5
SureMap	ℓ=3000	760 s	92.4	85.4	72.5	48.0	159.1
BWA-MEM	-x pacbio	180 s	90.9	84.0	70.7	46.1	159.8
BLASR	default	755 s	94.6	90.1	77.8	51.6	161.6
LAST	default	751 s	94.4	90.4	72.1	38.3	148.5

Peak memory: SureMap 8.9 GB, BWA-MEM 5.4 GB, BLASR 14.5 GB, LAST 12.5 GB.

-L buys accuracy with time almost linearly. At ℓ=650 SureMap is faster than BWA-MEM and aligns a slightly larger fraction of reads at every quality threshold; at ℓ=200 it is ~2× faster than BWA-MEM while still covering 81 % of what BWA-MEM covers.

Nanopore

30 143 reads, 176 Mbp total (ERR1676721), same setup.

Aligner	Option	Time	v ≤ 0.3	v ≤ 0.25	v ≤ 0.2	v ≤ 0.15	Coverage (Mbp)
SureMap	ℓ=200	187 s	90.0	87.1	80.0	61.2	145.4
SureMap	ℓ=300	204 s	91.9	89.5	82.8	63.3	155.8
SureMap	ℓ=650	230 s	93.1	91.1	84.8	65.0	164.1
SureMap	ℓ=1100	307 s	93.3	91.3	85.2	65.3	165.0
SureMap	ℓ=1600	411 s	93.4	91.3	85.2	65.3	165.4
SureMap	ℓ=3000	700 s	93.6	91.4	85.3	65.4	166.1
BWA-MEM	-x ont2d	441 s	94.4	91.7	84.3	64.1	168.0
BLASR	default	615 s	94.2	92.9	88.1	67.2	172.0

At ℓ=650 SureMap is ~1.9× faster than BWA-MEM and ~2.7× faster than BLASR, and beats both at the strictest thresholds (v ≤ 0.2 and v ≤ 0.15).

Variant calling on real data

Sample SRR099983 from the 1000 Genomes Project. BWA-MEM and Bowtie2 (very-sensitive) together map 95 % of the reads. Running SureMap with -u over the remaining ~12 M reads placed 2.5 % of them (~300 000 reads), and those reads cluster: VCFtools at depth ≥ 8 and quality > 80 calls 500 new variants that neither mapper had the coverage to support.

---

Repository layout

src/
  aligner_main.cpp        SureMap-Aligner entry point
  indexer_main.cpp        SureMap-Indexer entry point
  index/
    BitArray.h            bit vector with constant-time rank
    CompressedString.h    2-bits-per-base DNA string
    bwt.h                 FM-index: load and search
    bwt_builder.h         FM-index construction (pSAscan)
  align/
    common.h              shared aliases and limits
    types.h               BWT intervals, alignment records, DP steps
    globals.h/.cpp        index, options, read batch, per-thread state
    dp_row.h              banded DP row packed into 64-bit words
    row_cache.h           per-thread memo table for DP row transitions
    fm_search.h/.cpp      FM-index primitives and the pruning bounds
    pairwise_align.h/.cpp Hirschberg alignment, used to build CIGARs
    search.h/.cpp         divide-and-merge search: the core of the mapper
    cigar.h/.cpp          CIGAR construction
    sam_output.h/.cpp     SAM records and output ordering
    mappers.h             the three algorithms ...
    mapper_short.cpp        ... short
    mapper_hard.cpp         ... hard
    mapper_long.cpp         ... long
    options.h/.cpp        command line
    runner.h/.cpp         batching, worker threads, index loading
psascan_src/              vendored pSAscan, unmodified
tests/                    end-to-end smoke test on a synthetic genome

Per-thread scratch space is allocated for -t workers only, so a single-threaded run does not pay for the maximum of 40.

tests/run_tests.sh covers every algorithm end to end, including a multi-threaded run, and checks that reads land where they were simulated from.

Testing

tests/ builds a small synthetic genome, indexes it without needing pSAscan or libdivsufsort, simulates reads whose true position is encoded in the read name, and checks that the aligner places them there:

make test

Known limitations

• long shuffles its seed order from a clock-seeded RNG, so its output varies between runs. short and hard are deterministic and produce identical output at any -t.
• MAPQ is not computed; every aligned record reports 255.
• Paired-end reads are not supported; mates must be aligned separately.

Availability

SureMap is released under the MIT/X11 license. The latest source is at https://github.com/MohammadJRS/SureMap/.