Release minimap2-2.21 (r1071)

Give the CIGAR constants names to clarify the code. So that ksw2.h
remains self-contained, define KSW_* versions of the CIGAR operators
it needs for use within ksw2.h. Other code should in general use the
full set of MM_CIGAR_* constants in minimap.h.

NEWS.md

@@ -1,3 +1,26 @@
+Release 2.21-r1071 (6 July 2021)
+--------------------------------
+
+This release fixed a regression in short-read mapping introduced in v2.19
+(#776). It also fixed invalid comparisons of uninitialized variables, though
+these are harmless (#752). Long-read alignment should be identical to v2.20.
+
+(2.21: 6 July 2021)
+
+
+
+Release 2.20-r1061 (27 May 2021)
+--------------------------------
+
+This release fixed a bug in the Python module and improves the command-line
+compatibiliity with v2.18. In v2.19, if `-r` is specified with an `asm*` preset,
+users would get alignments more fragmented than v2.18. This could be an issue
+for existing pipelines specifying `-r`. This release resolves this issue.
+
+(2.20: 27 May 2021, r1061)
+
+
+
 Release 2.19-r1057 (26 May 2021)
 --------------------------------
 
@@ -33,7 +56,7 @@ unimap for contig alignment. Unimap will remain an experimental project and is
 no longer recommended over minimap2. Sorry for reverting the recommendation in
 short time.
 
-(2.20: 26 May 2021, r1057)
+(2.19: 26 May 2021, r1057)
 
 
 

README.md

@@ -14,7 +14,7 @@ cd minimap2 && make
 # use presets (no test data)
 ./minimap2 -ax map-pb ref.fa pacbio.fq.gz > aln.sam       # PacBio CLR genomic reads
 ./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam         # Oxford Nanopore genomic reads
-./minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz > aln.sam # PacBio HiFi/CCS genomic reads (v2.19 or lateer)
+./minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz > aln.sam # PacBio HiFi/CCS genomic reads (v2.19 or later)
 ./minimap2 -ax asm20 ref.fa pacbio-ccs.fq.gz > aln.sam    # PacBio HiFi/CCS genomic reads (v2.18 or earlier)
 ./minimap2 -ax sr ref.fa read1.fa read2.fa > aln.sam      # short genomic paired-end reads
 ./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam       # spliced long reads (strand unknown)
@@ -74,8 +74,8 @@ Detailed evaluations are available from the [minimap2 paper][doi] or the
 Minimap2 is optimized for x86-64 CPUs. You can acquire precompiled binaries from
 the [release page][release] with:
 ```sh
-curl -L https://github.com/lh3/minimap2/releases/download/v2.19/minimap2-2.19_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.19_x64-linux/minimap2
+curl -L https://github.com/lh3/minimap2/releases/download/v2.21/minimap2-2.21_x64-linux.tar.bz2 | tar -jxvf -
+./minimap2-2.21_x64-linux/minimap2
 ```
 If you want to compile from the source, you need to have a C compiler, GNU make
 and zlib development files installed. Then type `make` in the source code
@@ -96,7 +96,7 @@ with the ARM related command lines given above.
 
 Without any options, minimap2 takes a reference database and a query sequence
 file as input and produce approximate mapping, without base-level alignment
-(i.e. no CIGAR), in the [PAF format][paf]:
+(i.e. coordinates are only approximate and no CIGAR in output), in the [PAF format][paf]:
 ```sh
 minimap2 ref.fa query.fq > approx-mapping.paf
 ```

align.c

@@ -53,16 +53,16 @@ static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq,
 	// find the score and the region where score drops most along diagonal
 	for (k = 0, score = 0; k < n_cigar; ++k) {
 		uint32_t l, op = cigar[k]&0xf, len = cigar[k]>>4;
-		if (op == 0) {
+		if (op == MM_CIGAR_MATCH) {
 			for (l = 0; l < len; ++l) {
 				score += mat[tseq[i + l] * 5 + qseq[j + l]];
 				update_max_zdrop(score, i+l, j+l, &max, &max_i, &max_j, opt->e, &max_zdrop, pos);
 			}
 			i += len, j += len;
-		} else if (op == 1 || op == 2 || op == 3) {
+		} else if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL || op == MM_CIGAR_N_SKIP) {
 			score -= opt->q + opt->e * len;
-			if (op == 1) j += len; // insertion
-			else i += len;         // deletion
+			if (op == MM_CIGAR_INS) j += len;
+			else i += len;
 			update_max_zdrop(score, i, j, &max, &max_i, &max_j, opt->e, &max_zdrop, pos);
 		}
 	}
@@ -98,12 +98,12 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 	for (k = 0; k < p->n_cigar; ++k) { // indel left alignment
 		uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4;
 		if (len == 0) to_shrink = 1;
-		if (op == 0) {
+		if (op == MM_CIGAR_MATCH) {
 			toff += len, qoff += len;
-		} else if (op == 1 || op == 2) { // insertion or deletion
+		} else if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL) {
 			if (k > 0 && k < p->n_cigar - 1 && (p->cigar[k-1]&0xf) == 0 && (p->cigar[k+1]&0xf) == 0) {
 				int l, prev_len = p->cigar[k-1] >> 4;
-				if (op == 1) {
+				if (op == MM_CIGAR_INS) {
 					for (l = 0; l < prev_len; ++l)
 						if (qseq[qoff - 1 - l] != qseq[qoff + len - 1 - l])
 							break;
@@ -116,9 +116,9 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 					p->cigar[k-1] -= l<<4, p->cigar[k+1] += l<<4, qoff -= l, toff -= l;
 				if (l == prev_len) to_shrink = 1;
 			}
-			if (op == 1) qoff += len;
+			if (op == MM_CIGAR_INS) qoff += len;
 			else toff += len;
-		} else if (op == 3) {
+		} else if (op == MM_CIGAR_N_SKIP) {
 			toff += len;
 		}
 	}
@@ -128,13 +128,13 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 			uint32_t l, s[3] = {0,0,0};
 			for (l = k; l < p->n_cigar; ++l) { // count number of adjacent I and D
 				uint32_t op = p->cigar[l]&0xf;
-				if (op == 1 || op == 2 || p->cigar[l]>>4 == 0)
+				if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL || p->cigar[l]>>4 == 0)
 					s[op] += p->cigar[l] >> 4;
 				else break;
 			}
 			if (s[1] > 0 && s[2] > 0 && l - k > 2) { // turn to a single I and a single D
-				p->cigar[k]   = s[1]<<4|1;
-				p->cigar[k+1] = s[2]<<4|2;
+				p->cigar[k]   = s[1]<<4|MM_CIGAR_INS;
+				p->cigar[k+1] = s[2]<<4|MM_CIGAR_DEL;
 				for (k += 2; k < l; ++k)
 					p->cigar[k] &= 0xf;
 				to_shrink = 1;
@@ -154,9 +154,9 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 			else p->cigar[k+1] += p->cigar[k]>>4<<4; // add length to the next CIGAR operator
 		p->n_cigar = l;
 	}
-	if ((p->cigar[0]&0xf) == 1 || (p->cigar[0]&0xf) == 2) { // get rid of leading I or D
+	if ((p->cigar[0]&0xf) == MM_CIGAR_INS || (p->cigar[0]&0xf) == MM_CIGAR_DEL) { // get rid of leading I or D
 		int32_t l = p->cigar[0] >> 4;
-		if ((p->cigar[0]&0xf) == 1) {
+		if ((p->cigar[0]&0xf) == MM_CIGAR_INS) {
 			if (r->rev) r->qe -= l;
 			else r->qs += l;
 			*qshift = l;
@@ -174,7 +174,7 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 	if (r->p == 0) return;
 	for (k = 0; k < r->p->n_cigar; ++k) {
 		uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
-		if (op == 0) {
+		if (op == MM_CIGAR_MATCH) {
 			while (len > 0) {
 				for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) {} // run of "="; TODO: N<=>N is converted to "="
 				if (l > 0) { ++n_EQX; len -= l; toff += l; qoff += l; }
@@ -183,11 +183,11 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 				if (l > 0) { ++n_EQX; len -= l; toff += l; qoff += l; }
 			}
 			++n_M;
-		} else if (op == 1) { // insertion
+		} else if (op == MM_CIGAR_INS) {
 			qoff += len;
-		} else if (op == 2) { // deletion
+		} else if (op == MM_CIGAR_DEL) {
 			toff += len;
-		} else if (op == 3) { // intron
+		} else if (op == MM_CIGAR_N_SKIP) {
 			toff += len;
 		}
 	}
@@ -195,7 +195,7 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 	if (n_EQX == n_M) {
 		for (k = 0; k < r->p->n_cigar; ++k) {
 			uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
-			if (op == 0) r->p->cigar[k] = len << 4 | 7;
+			if (op == MM_CIGAR_MATCH) r->p->cigar[k] = len << 4 | MM_CIGAR_EQ_MATCH;
 		}
 		return;
 	}
@@ -209,25 +209,25 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 	toff = qoff = m = 0;
 	for (k = 0; k < r->p->n_cigar; ++k) {
 		uint32_t op = r->p->cigar[k]&0xf, len = r->p->cigar[k]>>4;
-		if (op == 0) { // match/mismatch
+		if (op == MM_CIGAR_MATCH) {
 			while (len > 0) {
 				// match
 				for (l = 0; l < len && qseq[qoff + l] == tseq[toff + l]; ++l) {}
-				if (l > 0) p->cigar[m++] = l << 4 | 7;
+				if (l > 0) p->cigar[m++] = l << 4 | MM_CIGAR_EQ_MATCH;
 				len -= l;
 				toff += l, qoff += l;
 				// mismatch
 				for (l = 0; l < len && qseq[qoff + l] != tseq[toff + l]; ++l) {}
-				if (l > 0) p->cigar[m++] = l << 4 | 8;
+				if (l > 0) p->cigar[m++] = l << 4 | MM_CIGAR_X_MISMATCH;
 				len -= l;
 				toff += l, qoff += l;
 			}
 			continue;
-		} else if (op == 1) { // insertion
+		} else if (op == MM_CIGAR_INS) {
 			qoff += len;
-		} else if (op == 2) { // deletion
+		} else if (op == MM_CIGAR_DEL) {
 			toff += len;
-		} else if (op == 3) { // intron
+		} else if (op == MM_CIGAR_N_SKIP) {
 			toff += len;
 		}
 		p->cigar[m++] = r->p->cigar[k];
@@ -248,7 +248,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
 	r->blen = r->mlen = 0;
 	for (k = 0; k < p->n_cigar; ++k) {
 		uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4;
-		if (op == 0) { // match/mismatch
+		if (op == MM_CIGAR_MATCH) {
 			int n_ambi = 0, n_diff = 0;
 			for (l = 0; l < len; ++l) {
 				int cq = qseq[qoff + l], ct = tseq[toff + l];
@@ -260,7 +260,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
 			}
 			r->blen += len - n_ambi, r->mlen += len - (n_ambi + n_diff), p->n_ambi += n_ambi;
 			toff += len, qoff += len;
-		} else if (op == 1) { // insertion
+		} else if (op == MM_CIGAR_INS) {
 			int n_ambi = 0;
 			for (l = 0; l < len; ++l)
 				if (qseq[qoff + l] > 3) ++n_ambi;
@@ -268,7 +268,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
 			s -= q + e * len;
 			if (s < 0) s = 0;
 			qoff += len;
-		} else if (op == 2) { // deletion
+		} else if (op == MM_CIGAR_DEL) {
 			int n_ambi = 0;
 			for (l = 0; l < len; ++l)
 				if (tseq[toff + l] > 3) ++n_ambi;
@@ -276,7 +276,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
 			s -= q + e * len;
 			if (s < 0) s = 0;
 			toff += len;
-		} else if (op == 3) { // intron
+		} else if (op == MM_CIGAR_N_SKIP) {
 			toff += len;
 		}
 	}
@@ -333,7 +333,7 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 		int i;
 		fprintf(stderr, "score=%d, cigar=", ez->score);
 		for (i = 0; i < ez->n_cigar; ++i)
-			fprintf(stderr, "%d%c", ez->cigar[i]>>4, "MIDN"[ez->cigar[i]&0xf]);
+			fprintf(stderr, "%d%c", ez->cigar[i]>>4, MM_CIGAR_STR[ez->cigar[i]&0xf]);
 		fprintf(stderr, "\n");
 	}
 }
@@ -730,7 +730,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 					if (qseq[j] >= 4 || tseq[j] >= 4) ez->score += opt->e2;
 					else ez->score += qseq[j] == tseq[j]? opt->a : -opt->b;
 				}
-				ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, 0, qe - qs);
+				ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, MM_CIGAR_MATCH, qe - qs);
 			} else { // perform normal gapped alignment
 				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
 			}

cookbook.md

@@ -31,8 +31,8 @@ To acquire the data used in this cookbook and to install minimap2 and paftools,
 please follow the command lines below:
 ```sh
 # install minimap2 executables
-curl -L https://github.com/lh3/minimap2/releases/download/v2.19/minimap2-2.19_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.19_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+curl -L https://github.com/lh3/minimap2/releases/download/v2.21/minimap2-2.21_x64-linux.tar.bz2 | tar jxf -
+cp minimap2-2.21_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
 export PATH="$PATH:"`pwd`                               # put the current directory on PATH
 # download example datasets
 curl -L https://github.com/lh3/minimap2/releases/download/v2.10/cookbook-data.tgz | tar zxf -

example.c

@@ -47,7 +47,7 @@ int main(int argc, char *argv[])
 				printf("%s\t%d\t%d\t%d\t%c\t", ks->name.s, ks->seq.l, r->qs, r->qe, "+-"[r->rev]);
 				printf("%s\t%d\t%d\t%d\t%d\t%d\t%d\tcg:Z:", mi->seq[r->rid].name, mi->seq[r->rid].len, r->rs, r->re, r->mlen, r->blen, r->mapq);
 				for (i = 0; i < r->p->n_cigar; ++i) // IMPORTANT: this gives the CIGAR in the aligned regions. NO soft/hard clippings!
-					printf("%d%c", r->p->cigar[i]>>4, "MIDNSH"[r->p->cigar[i]&0xf]);
+					printf("%d%c", r->p->cigar[i]>>4, MM_CIGAR_STR[r->p->cigar[i]&0xf]);
 				putchar('\n');
 				free(r->p);
 			}

format.c

@@ -144,8 +144,8 @@ static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 	if (write_tag) mm_sprintf_lite(s, "\tcs:Z:");
 	for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
 		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
-		assert((op >= 0 && op <= 3) || op == 7 || op == 8);
-		if (op == 0 || op == 7 || op == 8) { // match
+		assert((op >= MM_CIGAR_MATCH && op <= MM_CIGAR_N_SKIP) || op == MM_CIGAR_EQ_MATCH || op == MM_CIGAR_X_MISMATCH);
+		if (op == MM_CIGAR_MATCH || op == MM_CIGAR_EQ_MATCH || op == MM_CIGAR_X_MISMATCH) {
 			int l_tmp = 0;
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
@@ -166,12 +166,12 @@ static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 				} else mm_sprintf_lite(s, ":%d", l_tmp);
 			}
 			q_off += len, t_off += len;
-		} else if (op == 1) { // insertion to ref
+		} else if (op == MM_CIGAR_INS) {
 			for (j = 0, tmp[len] = 0; j < len; ++j)
 				tmp[j] = "acgtn"[qseq[q_off + j]];
 			mm_sprintf_lite(s, "+%s", tmp);
 			q_off += len;
-		} else if (op == 2) { // deletion from ref
+		} else if (op == MM_CIGAR_DEL) {
 			for (j = 0, tmp[len] = 0; j < len; ++j)
 				tmp[j] = "acgtn"[tseq[t_off + j]];
 			mm_sprintf_lite(s, "-%s", tmp);
@@ -192,8 +192,8 @@ static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 	if (write_tag) mm_sprintf_lite(s, "\tMD:Z:");
 	for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
 		int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
-		assert((op >= 0 && op <= 3) || op == 7 || op == 8);
-		if (op == 0 || op == 7 || op == 8) { // match
+		assert((op >= MM_CIGAR_MATCH && op <= MM_CIGAR_N_SKIP) || op == MM_CIGAR_EQ_MATCH || op == MM_CIGAR_X_MISMATCH);
+		if (op == MM_CIGAR_MATCH || op == MM_CIGAR_EQ_MATCH || op == MM_CIGAR_X_MISMATCH) {
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
 					mm_sprintf_lite(s, "%d%c", l_MD, "ACGTN"[tseq[t_off + j]]);
@@ -201,15 +201,15 @@ static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
 				} else ++l_MD;
 			}
 			q_off += len, t_off += len;
-		} else if (op == 1) { // insertion to ref
+		} else if (op == MM_CIGAR_INS) {
 			q_off += len;
-		} else if (op == 2) { // deletion from ref
+		} else if (op == MM_CIGAR_DEL) {
 			for (j = 0, tmp[len] = 0; j < len; ++j)
 				tmp[j] = "ACGTN"[tseq[t_off + j]];
 			mm_sprintf_lite(s, "%d^%s", l_MD, tmp);
 			l_MD = 0;
 			t_off += len;
-		} else if (op == 3) { // reference skip
+		} else if (op == MM_CIGAR_N_SKIP) {
 			t_off += len;
 		}
 	}
@@ -271,7 +271,7 @@ double mm_event_identity(const mm_reg1_t *r)
 	if (r->p == 0) return -1.0f;
 	for (i = 0; i < r->p->n_cigar; ++i) {
 		int32_t op = r->p->cigar[i] & 0xf, len = r->p->cigar[i] >> 4;
-		if (op == 1 || op == 2)
+		if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL)
 			++n_gapo, n_gap += len;
 	}
 	return (double)r->mlen / (r->blen + r->p->n_ambi - n_gap + n_gapo);
@@ -325,7 +325,7 @@ void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const
 		uint32_t k;
 		mm_sprintf_lite(s, "\tcg:Z:");
 		for (k = 0; k < r->p->n_cigar; ++k)
-			mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDNSHP=XB"[r->p->cigar[k]&0xf]);
+			mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, MM_CIGAR_STR[r->p->cigar[k]&0xf]);
 	}
 	if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_MD)))
 		write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1);
@@ -382,7 +382,7 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 			assert(clip_len[0] < qlen && clip_len[1] < qlen);
 			if (clip_len[0]) mm_sprintf_lite(s, "%d%c", clip_len[0], clip_char);
 			for (k = 0; k < r->p->n_cigar; ++k)
-				mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDNSHP=XB"[r->p->cigar[k]&0xf]);
+				mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, MM_CIGAR_STR[r->p->cigar[k]&0xf]);
 			if (clip_len[1]) mm_sprintf_lite(s, "%d%c", clip_len[1], clip_char);
 		}
 	}

krmq.h

@@ -247,7 +247,7 @@ int main(void) {
 		unsigned char dir[KRMQ_MAX_DEPTH]; \
 		int i, d = 0, cmp; \
 		unsigned cnt = 0; \
-		fake.__head.p[0] = *root_, fake.__head.p[1] = 0; \
+		fake = **root_, fake.__head.p[0] = *root_, fake.__head.p[1] = 0; \
 		if (cnt_) *cnt_ = 0; \
 		if (x) { \
 			for (cmp = -1, p = &fake; cmp; cmp = __cmp(x, p)) { \

ksw2.h

@@ -16,6 +16,13 @@
 #define KSW_EZ_SPLICE_REV  0x200
 #define KSW_EZ_SPLICE_FLANK 0x400
 
+// The subset of CIGAR operators used by ksw code.
+// Use MM_CIGAR_* from minimap.h if you need the full list.
+#define KSW_CIGAR_MATCH  0
+#define KSW_CIGAR_INS    1
+#define KSW_CIGAR_DEL    2
+#define KSW_CIGAR_N_SKIP 3
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -137,13 +144,13 @@ static inline void ksw_backtrack(void *km, int is_rot, int is_rev, int min_intro
 		else if (!(tmp >> (state + 2) & 1)) state = 0; // if requesting other states, _state_ stays the same if it is a continuation; otherwise, set to H
 		if (state == 0) state = tmp & 7; // TODO: probably this line can be merged into the "else if" line right above; not 100% sure
 		if (force_state >= 0) state = force_state;
-		if (state == 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 0, 1), --i, --j; // match
-		else if (state == 1 || (state == 3 && min_intron_len <= 0)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 2, 1), --i; // deletion
-		else if (state == 3 && min_intron_len > 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 3, 1), --i; // intron
-		else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, 1), --j; // insertion
+		if (state == 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_MATCH, 1), --i, --j;
+		else if (state == 1 || (state == 3 && min_intron_len <= 0)) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_DEL, 1), --i;
+		else if (state == 3 && min_intron_len > 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_N_SKIP, 1), --i;
+		else cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_INS, 1), --j;
 	}
-	if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, min_intron_len > 0 && i >= min_intron_len? 3 : 2, i + 1); // first deletion
-	if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, 1, j + 1); // first insertion
+	if (i >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, min_intron_len > 0 && i >= min_intron_len? KSW_CIGAR_N_SKIP : KSW_CIGAR_DEL, i + 1); // first deletion
+	if (j >= 0) cigar = ksw_push_cigar(km, &n_cigar, &m_cigar, cigar, KSW_CIGAR_INS, j + 1); // first insertion
 	if (!is_rev)
 		for (i = 0; i < n_cigar>>1; ++i) // reverse CIGAR
 			tmp = cigar[i], cigar[i] = cigar[n_cigar-1-i], cigar[n_cigar-1-i] = tmp;

lchain.c

@@ -114,8 +114,9 @@ static inline int32_t comput_sc(const mm128_t *ai, const mm128_t *aj, int32_t ma
 		float lin_pen, log_pen;
 		lin_pen = chn_pen_gap * (float)dd + chn_pen_skip * (float)dg;
 		log_pen = dd >= 1? mg_log2(dd + 1) : 0.0f; // mg_log2() only works for dd>=2
-		if (is_cdna) {
-			if (dr > dq) sc -= (int)(lin_pen < log_pen? lin_pen : log_pen); // deletion or jump between paired ends
+		if (is_cdna || sidi != sidj) {
+			if (sidi != sidj && dr == 0) ++sc; // possibly due to overlapping paired ends; give a minor bonus
+			else if (dr > dq || sidi != sidj) sc -= (int)(lin_pen < log_pen? lin_pen : log_pen); // deletion or jump between paired ends
 			else sc -= (int)(lin_pen + .5f * log_pen);
 		} else sc -= (int)(lin_pen + .5f * log_pen);
 	}

main.c

@@ -7,7 +7,7 @@
 #include "mmpriv.h"
 #include "ketopt.h"
 
-#define MM_VERSION "2.19-r1057"
+#define MM_VERSION "2.21-r1071"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -317,7 +317,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -g NUM       stop chain enlongation if there are no minimizers in INT-bp [%d]\n", opt.max_gap);
 		fprintf(fp_help, "    -G NUM       max intron length (effective with -xsplice; changing -r) [200k]\n");
 		fprintf(fp_help, "    -F NUM       max fragment length (effective with -xsr or in the fragment mode) [800]\n");
-		fprintf(fp_help, "    -r NUM       bandwidth used in chaining and DP-based alignment [%d]\n", opt.bw);
+		fprintf(fp_help, "    -r NUM[,NUM] chaining/alignment bandwidth and long-join bandwidth [%d,%d]\n", opt.bw, opt.bw_long);
 		fprintf(fp_help, "    -n INT       minimal number of minimizers on a chain [%d]\n", opt.min_cnt);
 		fprintf(fp_help, "    -m INT       minimal chaining score (matching bases minus log gap penalty) [%d]\n", opt.min_chain_score);
 //		fprintf(fp_help, "    -T INT       SDUST threshold; 0 to disable SDUST [%d]\n", opt.sdust_thres); // TODO: this option is never used; might be buggy

map.c

@@ -277,7 +277,17 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 						 opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 	}
 
-	if (opt->max_occ > opt->mid_occ && rep_len > 0 && !(opt->flag & MM_F_RMQ)) {
+	if (opt->bw_long > opt->bw && (opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN)) == 0 && n_segs == 1 && n_regs0 > 1) { // re-chain/long-join for long sequences
+		int32_t st = (int32_t)a[0].y, en = (int32_t)a[(int32_t)u[0] - 1].y;
+		if (qlen_sum - (en - st) > opt->rmq_rescue_size || en - st > qlen_sum * opt->rmq_rescue_ratio) {
+			int32_t i;
+			for (i = 0, n_a = 0; i < n_regs0; ++i) n_a += (int32_t)u[i];
+			kfree(b->km, u);
+			radix_sort_128x(a, a + n_a);
+			a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
+							  opt->chain_gap_scale * 0.01 * mi->k, 0.0f, n_a, a, &n_regs0, &u, b->km);
+		}
+	} else if (opt->max_occ > opt->mid_occ && rep_len > 0 && !(opt->flag & MM_F_RMQ)) { // re-chain, mostly for short reads
 		int rechain = 0;
 		if (n_regs0 > 0) { // test if the best chain has all the segments
 			int n_chained_segs = 1, max = 0, max_i = -1, max_off = -1, off = 0;
@@ -300,16 +310,6 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			a = mg_lchain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->max_chain_iter, opt->min_cnt, opt->min_chain_score,
 							 opt->chain_gap_scale * 0.01 * mi->k, 0.0f, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
-	} else if (opt->bw_long > opt->bw && (opt->flag & (MM_F_RMQ|MM_F_NO_LJOIN)) == 0 && n_segs == 1 && n_regs0 > 1) {
-		int32_t st = (int32_t)a[0].y, en = (int32_t)a[(int32_t)u[0] - 1].y;
-		if (qlen_sum - (en - st) > opt->rmq_rescue_size || en - st > qlen_sum * opt->rmq_rescue_ratio) {
-			int32_t i;
-			for (i = 0, n_a = 0; i < n_regs0; ++i) n_a += (int32_t)u[i];
-			kfree(b->km, u);
-			radix_sort_128x(a, a + n_a);
-			a = mg_lchain_rmq(opt->max_gap, opt->rmq_inner_dist, opt->bw_long, opt->max_chain_skip, opt->rmq_size_cap, opt->min_cnt, opt->min_chain_score,
-							  opt->chain_gap_scale * 0.01 * mi->k, 0.0f, n_a, a, &n_regs0, &u, b->km);
-		}
 	}
 	b->frag_gap = max_chain_gap_ref;
 	b->rep_len = rep_len;

minimap.h

@@ -46,6 +46,18 @@
 
 #define MM_MAX_SEG       255
 
+#define MM_CIGAR_MATCH      0
+#define MM_CIGAR_INS        1
+#define MM_CIGAR_DEL        2
+#define MM_CIGAR_N_SKIP     3
+#define MM_CIGAR_SOFTCLIP   4
+#define MM_CIGAR_HARDCLIP   5
+#define MM_CIGAR_PADDING    6
+#define MM_CIGAR_EQ_MATCH   7
+#define MM_CIGAR_X_MISMATCH 8
+
+#define MM_CIGAR_STR  "MIDNSHP=XB"
+
 #ifdef __cplusplus
 extern "C" {
 #endif

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "26 May 2021" "minimap2-2.19 (r1057)" "Bioinformatics tools"
+.TH minimap2 1 "6 July 2021" "minimap2-2.21 (r1071)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -165,9 +165,12 @@ Stop chain enlongation if there are no minimizers within
 .IR NUM -bp
 [10k].
 .TP
-.BI -r \ NUM
-Bandwidth used in chaining and DP-based alignment [500,20k]. This option
-approximately controls the maximum gap size.
+.BI -r \ NUM1 [, NUM2 ]
+Bandwidth for chaining and base alignment [500,20k].
+.I NUM1
+is used for initial chaining and alignment extension;
+.I NUM2
+for RMQ-based re-chaining and closing gaps in alignments.
 .TP
 .BI -n \ INT
 Discard chains consisting of
@@ -589,8 +592,8 @@ Long-read splice alignment for PacBio CCS reads
 .B sr
 Short single-end reads without splicing
 .RB ( -k21
-.B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -r50 -p.5 -N20 -f1000,5000 -n2 -m20
-.B -s40 -g200 -2K50m --heap-sort=yes
+.B -w11 --sr --frag=yes -A2 -B8 -O12,32 -E2,1 -r100 -p.5 -N20 -f1000,5000 -n2 -m20
+.B -s40 -g100 -2K50m --heap-sort=yes
 .BR --secondary=no ).
 .TP
 .B ava-pb

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.19-r1057';
+var paftools_version = '2.21-r1071';
 
 /*****************************
  ***** Library functions *****
@@ -1419,7 +1419,7 @@ function paf_view(args)
 
 	var s_ref = new Bytes(), s_qry = new Bytes(), s_mid = new Bytes(); // these are used to show padded alignment
 	var re_cs = /([:=\-\+\*])(\d+|[A-Za-z]+)/g;
-	var re_cg = /(\d+)([MIDNSH])/g;
+	var re_cg = /(\d+)([MIDNSHP=X])/g;
 
 	var buf = new Bytes();
 	var file = args[getopt.ind] == "-"? new File() : new File(args[getopt.ind]);
@@ -1899,7 +1899,7 @@ function paf_splice2bed(args)
 		a.length = 0;
 	}
 
-	var re = /(\d+)([MIDNSH])/g;
+	var re = /(\d+)([MIDNSHP=X])/g;
 	var c, fmt = "bed", fn_name_conv = null, keep_multi = false;
 	while ((c = getopt(args, "f:n:m")) != null) {
 		if (c == 'f') fmt = getopt.arg;
@@ -2369,7 +2369,7 @@ function paf_junceval(args)
 
 	file = getopt.ind+1 >= args.length || args[getopt.ind+1] == '-'? new File() : new File(args[getopt.ind+1]);
 	var last_qname = null;
-	var re_cigar = /(\d+)([MIDNSHX=])/g;
+	var re_cigar = /(\d+)([MIDNSHP=X])/g;
 	while (file.readline(buf) >= 0) {
 		var m, t = buf.toString().split("\t");
 		var ctg_name = null, cigar = null, pos = null, qname = t[0];

options.c

@@ -75,7 +75,7 @@ void mm_mapopt_update(mm_mapopt_t *opt, const mm_idx_t *mi)
 void mm_mapopt_max_intron_len(mm_mapopt_t *opt, int max_intron_len)
 {
 	if ((opt->flag & MM_F_SPLICE) && max_intron_len > 0)
-		opt->max_gap_ref = opt->bw = max_intron_len;
+		opt->max_gap_ref = opt->bw = opt->bw_long = max_intron_len;
 }
 
 int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
@@ -109,7 +109,7 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		io->flag = 0, io->k = 19, io->w = 19;
 		mo->bw = mo->bw_long = 100000;
 		mo->max_gap = 10000;
-		mo->flag |= MM_F_RMQ | MM_F_NO_LJOIN;
+		mo->flag |= MM_F_RMQ;
 		mo->min_mid_occ = 50, mo->max_mid_occ = 500;
 		mo->min_dp_max = 200;
 		mo->best_n = 50;
@@ -156,6 +156,11 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 
 int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
 {
+	if (mo->bw > mo->bw_long) {
+		if (mm_verbose >= 1)
+			fprintf(stderr, "[ERROR]\033[1;31m with '-rNUM1,NUM2', NUM1 (%d) can't be larger than NUM2 (%d)\033[0m\n", mo->bw, mo->bw_long);
+		return -8;
+	}
 	if ((mo->flag & MM_F_RMQ) && (mo->flag & (MM_F_SR|MM_F_SPLICE))) {
 		if (mm_verbose >= 1)
 			fprintf(stderr, "[ERROR]\033[1;31m --rmq doesn't work with --sr or --splice\033[0m\n");

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.19'
+__version__ = '2.21'
 
 cmappy.mm_reset_timer()
 
@@ -82,7 +82,7 @@ cdef class Alignment:
 
 	@property
 	def cigar_str(self):
-		return "".join(map(lambda x: str(x[0]) + 'MIDNSH'[x[1]], self._cigar))
+		return "".join(map(lambda x: str(x[0]) + 'MIDNSHP=XB'[x[1]], self._cigar))
 
 	def __str__(self):
 		if self._strand > 0: strand = '+'

seed.c

@@ -46,19 +46,20 @@ void mm_seed_select(int32_t n, mm_seed_t *a, int len, int max_occ, int max_max_o
 				int32_t pe = i == n? len : (uint32_t)a[i].q_pos>>1;
 				int32_t j, k, st = last0 + 1, en = i;
 				int32_t max_high_occ = (int32_t)((double)(pe - ps) / dist + .499);
-				//fprintf(stderr, "Y\t%d\t%d\n", ps, pe);
-				if (max_high_occ > MAX_MAX_HIGH_OCC)
-					max_high_occ = MAX_MAX_HIGH_OCC;
-				for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)
-					b[k] = (uint64_t)a[j].n<<32 | j;
-				ks_heapmake_uint64_t(k, b); // initialize the binomial heap
-				for (; j < en; ++j) { // if there are more, choose top max_high_occ
-					if (a[j].n < (int32_t)(b[0]>>32)) { // then update the heap
-						b[0] = (uint64_t)a[j].n<<32 | j;
-						ks_heapdown_uint64_t(0, k, b);
+				if (max_high_occ > 0) {
+					if (max_high_occ > MAX_MAX_HIGH_OCC)
+						max_high_occ = MAX_MAX_HIGH_OCC;
+					for (j = st, k = 0; j < en && k < max_high_occ; ++j, ++k)
+						b[k] = (uint64_t)a[j].n<<32 | j;
+					ks_heapmake_uint64_t(k, b); // initialize the binomial heap
+					for (; j < en; ++j) { // if there are more, choose top max_high_occ
+						if (a[j].n < (int32_t)(b[0]>>32)) { // then update the heap
+							b[0] = (uint64_t)a[j].n<<32 | j;
+							ks_heapdown_uint64_t(0, k, b);
+						}
 					}
+					for (j = 0; j < k; ++j) a[(uint32_t)b[j]].flt = 1;
 				}
-				for (j = 0; j < k; ++j) a[(uint32_t)b[j]].flt = 1;
 				for (j = st; j < en; ++j) a[j].flt ^= 1;
 				for (j = st; j < en; ++j)
 					if (a[j].n > max_max_occ)

setup.py

@@ -23,7 +23,7 @@ def readme():
 
 setup(
 	name = 'mappy',
-	version = '2.19',
+	version = '2.21',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),
@@ -33,7 +33,7 @@ setup(
 	keywords = 'sequence-alignment',
 	scripts = ['python/minimap2.py'],
 	ext_modules = [Extension('mappy',
-		sources = ['python/mappy.pyx', 'align.c', 'bseq.c', 'lchain.c', 'format.c', 'hit.c', 'index.c', 'pe.c', 'options.c',
+		sources = ['python/mappy.pyx', 'align.c', 'bseq.c', 'lchain.c', 'seed.c', 'format.c', 'hit.c', 'index.c', 'pe.c', 'options.c',
 				   'ksw2_extd2_sse.c', 'ksw2_exts2_sse.c', 'ksw2_extz2_sse.c', 'ksw2_ll_sse.c',
 				   'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c', 'esterr.c', 'splitidx.c'],
 		depends = ['minimap.h', 'bseq.h', 'kalloc.h', 'kdq.h', 'khash.h', 'kseq.h', 'ksort.h',