debugging code

NEWS.md

@@ -1,142 +1,3 @@
-Release 2.17-r941 (4 May 2019)
-------------------------------
-
-Changes since the last release:
-
- * Fixed flawed CIGARs like `5I6D7I` (#392).
-
- * Bugfix: TLEN should be 0 when either end is unmapped (#373 and #365).
-
- * Bugfix: mappy is unable to write index (#372).
-
- * Added option `--junc-bed` to load known gene annotations in the BED12
-   format. Minimap2 prefers annotated junctions over novel junctions (#197 and
-   #348). GTF can be converted to BED12 with `paftools.js gff2bed`.
-
- * Added option `--sam-hit-only` to suppress unmapped hits in SAM (#377).
-
- * Added preset `splice:hq` for high-quality CCS or mRNA sequences. It applies
-   better scoring and improves the sensitivity to small exons. This preset may
-   introduce false small introns, but the overall accuracy should be higher.
-
-This version produces nearly identical alignments to v2.16, except for CIGARs
-affected by the bug mentioned above.
-
-(2.17: 5 May 2019, r941)
-
-
-
-Release 2.16-r922 (28 February 2019)
-------------------------------------
-
-This release is 50% faster for mapping ultra-long nanopore reads at comparable
-accuracy. For short-read mapping, long-read overlapping and ordinary long-read
-mapping, the performance and accuracy remain similar. This speedup is achieved
-with a new heuristic to limit the number of chaining iterations (#324). Users
-can disable the heuristic by increasing a new option `--max-chain-iter` to a
-huge number.
-
-Other changes to minimap2:
-
- * Implemented option `--paf-no-hit` to output unmapped query sequences in PAF.
-   The strand and reference name columns are both `*` at an unmapped line. The
-   hidden option is available in earlier minimap2 but had a different 2-column
-   output format instead of PAF.
-
- * Fixed a bug that leads to wrongly calculated `de` tags when ambiguous bases
-   are involved (#309). This bug only affects v2.15.
-
- * Fixed a bug when parsing command-line option `--splice` (#344). This bug was
-   introduced in v2.13.
-
- * Fixed two division-by-zero cases (#326). They don't affect final alignments
-   because the results of the divisions are not used in both case.
-
- * Added an option `-o` to output alignments to a specified file. It is still
-   recommended to use UNIX pipes for on-the-fly conversion or compression.
-
- * Output a new `rl` tag to give the length of query regions harboring
-   repetitive seeds.
-
-Changes to paftool.js:
-
- * Added a new option to convert the MD tag to the long form of the cs tag.
-
-Changes to mappy:
-
- * Added the `mappy.Aligner.seq_names` method to return sequence names (#312).
-
-For NA12878 ultra-long reads, this release changes the alignments of <0.1% of
-reads in comparison to v2.15. All these reads have highly fragmented alignments
-and are likely to be problematic anyway. For shorter or well aligned reads,
-this release should produce mostly identical alignments to v2.15.
-
-(2.16: 28 February 2019, r922)
-
-
-
-Release 2.15-r905 (10 January 2019)
------------------------------------
-
-Changes to minimap2:
-
- * Fixed a rare segmentation fault when option -H is in use (#307). This may
-   happen when there are very long homopolymers towards the 5'-end of a read.
-
- * Fixed wrong CIGARs when option --eqx is used (#266).
-
- * Fixed a typo in the base encoding table (#264). This should have no
-   practical effect.
-
- * Fixed a typo in the example code (#265).
-
- * Improved the C++ compatibility by removing "register" (#261). However,
-   minimap2 still can't be compiled in the pedantic C++ mode (#306).
-
- * Output a new "de" tag for gap-compressed sequence divergence.
-
-Changes to paftools.js:
-
- * Added "asmgene" to evaluate the completeness of an assembly by measuring the
-   uniquely mapped single-copy genes. This command learns the idea of BUSCO.
-
- * Added "vcfpair" to call a phased VCF from phased whole-genome assemblies. An
-   earlier version of this script is used to produce the ground truth for the
-   syndip benchmark [PMID:30013044].
-
-This release produces identical alignment coordinates and CIGARs in comparison
-to v2.14. Users are advised to upgrade due to the several bug fixes.
-
-(2.15: 10 Janurary 2019, r905)
-
-
-
-Release 2.14-r883 (5 November 2018)
------------------------------------
-
-Notable changes:
-
- * Fixed two minor bugs caused by typos (#254 and #266).
-
- * Fixed a bug that made minimap2 abort when --eqx was used together with --MD
-   or --cs (#257).
-
- * Added --cap-sw-mem to cap the size of DP matrices (#259). Base alignment may
-   take a lot of memory in the splicing mode. This may lead to issues when we
-   run minimap2 on a cluster with a hard memory limit. The new option avoids
-   unlimited memory usage at the cost of missing a few long introns.
-
- * Conforming to C99 and C11 when possible (#261).
-
- * Warn about malformatted FASTA or FASTQ (#252 and #255).
-
-This release occasionally produces base alignments different from v2.13. The
-overall alignment accuracy remain similar.
-
-(2.14: 5 November 2018, r883)
-
-
-
 Release 2.13-r850 (11 October 2018)
 -----------------------------------
 

README.md

@@ -9,8 +9,8 @@ cd minimap2 && make
 # long sequences against a reference genome
 ./minimap2 -a test/MT-human.fa test/MT-orang.fa > test.sam
 # create an index first and then map
-./minimap2 -x map-ont -d MT-human-ont.mmi test/MT-human.fa
+./minimap2 -d MT-human.mmi test/MT-human.fa
-./minimap2 -a MT-human-ont.mmi test/MT-orang.fa > test.sam
+./minimap2 -a MT-human.mmi test/MT-orang.fa > test.sam
 # use presets (no test data)
 ./minimap2 -ax map-pb ref.fa pacbio.fq.gz > aln.sam       # PacBio genomic reads
 ./minimap2 -ax map-ont ref.fa ont.fq.gz > aln.sam         # Oxford Nanopore genomic reads
@@ -18,7 +18,7 @@ cd minimap2 && make
 ./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)
 ./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam  # noisy Nanopore Direct RNA-seq
-./minimap2 -ax splice:hq -uf ref.fa query.fa > aln.sam    # Final PacBio Iso-seq or traditional cDNA
+./minimap2 -ax splice -uf -C5 ref.fa query.fa > aln.sam   # Final PacBio Iso-seq or traditional cDNA
 ./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf             # intra-species asm-to-asm alignment
 ./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf     # PacBio read overlap
 ./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf    # Nanopore read overlap
@@ -71,8 +71,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.17/minimap2-2.17_x64-linux.tar.bz2 | tar -jxvf -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.13/minimap2-2.13_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.17_x64-linux/minimap2
+./minimap2-2.13_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
@@ -139,7 +139,7 @@ Nanopore reads.
 #### <a name="map-long-splice"></a>Map long mRNA/cDNA reads
 
 ```sh
-minimap2 -ax splice:hq -uf ref.fa iso-seq.fq > aln.sam       # PacBio Iso-seq/traditional cDNA
+minimap2 -ax splice -uf -C5 ref.fa iso-seq.fq > aln.sam      # PacBio Iso-seq/traditional cDNA
 minimap2 -ax splice ref.fa nanopore-cdna.fa > aln.sam        # Nanopore 2D cDNA-seq
 minimap2 -ax splice -uf -k14 ref.fa direct-rna.fq > aln.sam  # Nanopore Direct RNA-seq
 minimap2 -ax splice --splice-flank=no SIRV.fa SIRV-seq.fa    # mapping against SIRV control
@@ -324,7 +324,7 @@ There is not a specific mailing list for the time being.
 If you use minimap2 in your work, please cite:
 
 > Li, H. (2018). Minimap2: pairwise alignment for nucleotide sequences.
-> *Bioinformatics*, **34**:3094-3100. [doi:10.1093/bioinformatics/bty191][doi]
+> Bioinformatics. [doi:10.1093/bioinformatics/bty191][doi]
 
 ## <a name="dguide"></a>Developers' Guide
 
@@ -355,8 +355,6 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
   billion bases or longer (2,147,483,647 to be exact). The total length of all
   sequences can well exceed this threshold.
 
-* Minimap2 often misses small exons.
-
 
 
 [paf]: https://github.com/lh3/miniasm/blob/master/PAF.md

align.c

@@ -88,13 +88,14 @@ static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq,
 	return max_zdrop > opt->zdrop? 1 : 0;
 }
 
-static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift)
+static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift, int left_aln)
 {
 	mm_extra_t *p = r->p;
 	int32_t toff = 0, qoff = 0, to_shrink = 0;
 	uint32_t k;
 	*qshift = *tshift = 0;
 	if (p->n_cigar <= 1) return;
+	if (!left_aln) goto end_left_aln;
 	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;
@@ -123,25 +124,7 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 		}
 	}
 	assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
-	for (k = 0; k < p->n_cigar - 2; ++k) { // fix CIGAR like 5I6D7I
+end_left_aln:
-		if ((p->cigar[k]&0xf) > 0 && (p->cigar[k]&0xf) + (p->cigar[k+1]&0xf) == 3) {
-			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)
-					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;
-				for (k += 2; k < l; ++k)
-					p->cigar[k] &= 0xf;
-				to_shrink = 1;
-			}
-			k = l;
-		}
-	}
 	if (to_shrink) { // squeeze out zero-length operations
 		int32_t l = 0;
 		for (k = 0; k < p->n_cigar; ++k) // squeeze out zero-length operations
@@ -166,6 +149,78 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq,
 	}
 }
 
+static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e, int left_aln)
+{
+	uint32_t k, l;
+	int32_t s = 0, max = 0, qshift, tshift, toff = 0, qoff = 0;
+	mm_extra_t *p = r->p;
+	if (p == 0) return;
+	mm_fix_cigar(r, qseq, tseq, &qshift, &tshift, left_aln);
+	qseq += qshift, tseq += tshift; // qseq and tseq may be shifted due to the removal of leading I/D
+	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
+			int n_ambi = 0, n_diff = 0;
+			for (l = 0; l < len; ++l) {
+				int cq = qseq[qoff + l], ct = tseq[toff + l];
+				if (ct > 3 || cq > 3) ++n_ambi;
+				else if (ct != cq) ++n_diff;
+				s += mat[ct * 5 + cq];
+				if (s < 0) s = 0;
+				else max = max > s? max : s;
+			}
+			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
+			int n_ambi = 0;
+			for (l = 0; l < len; ++l)
+				if (qseq[qoff + l] > 3) ++n_ambi;
+			r->blen += len - n_ambi, p->n_ambi += n_ambi;
+			s -= q + e * len;
+			if (s < 0) s = 0;
+			qoff += len;
+		} else if (op == 2) { // deletion
+			int n_ambi = 0;
+			for (l = 0; l < len; ++l)
+				if (tseq[toff + l] > 3) ++n_ambi;
+			r->blen += len - n_ambi, p->n_ambi += n_ambi;
+			s -= q + e * len;
+			if (s < 0) s = 0;
+			toff += len;
+		} else if (op == 3) { // intron
+			toff += len;
+		}
+	}
+	p->dp_max = max;
+	assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
+}
+
+static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) // TODO: this calls the libc realloc()
+{
+	mm_extra_t *p;
+	if (n_cigar == 0) return;
+	if (r->p == 0) {
+		uint32_t capacity = n_cigar + sizeof(mm_extra_t)/4;
+		kroundup32(capacity);
+		r->p = (mm_extra_t*)calloc(capacity, 4);
+		r->p->capacity = capacity;
+	} else if (r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4 > r->p->capacity) {
+		r->p->capacity = r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4;
+		kroundup32(r->p->capacity);
+		r->p = (mm_extra_t*)realloc(r->p, r->p->capacity * 4);
+	}
+	p = r->p;
+	if (p->n_cigar > 0 && (p->cigar[p->n_cigar-1]&0xf) == (cigar[0]&0xf)) { // same CIGAR op at the boundary
+		p->cigar[p->n_cigar-1] += cigar[0]>>4<<4;
+		if (n_cigar > 1) memcpy(p->cigar + p->n_cigar, cigar + 1, (n_cigar - 1) * 4);
+		p->n_cigar += n_cigar - 1;
+	} else {
+		memcpy(p->cigar + p->n_cigar, cigar, n_cigar * 4);
+		p->n_cigar += n_cigar;
+	}
+}
+
 static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq) // written by @armintoepfer
 {
 	uint32_t n_EQX = 0;
@@ -237,80 +292,8 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
 	r->p = p;
 }
 
-static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e, int is_eqx)
+static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *ob, const int8_t *mat,
-{
+						  int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
-	uint32_t k, l;
-	int32_t s = 0, max = 0, qshift, tshift, toff = 0, qoff = 0;
-	mm_extra_t *p = r->p;
-	if (p == 0) return;
-	mm_fix_cigar(r, qseq, tseq, &qshift, &tshift);
-	qseq += qshift, tseq += tshift; // qseq and tseq may be shifted due to the removal of leading I/D
-	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
-			int n_ambi = 0, n_diff = 0;
-			for (l = 0; l < len; ++l) {
-				int cq = qseq[qoff + l], ct = tseq[toff + l];
-				if (ct > 3 || cq > 3) ++n_ambi;
-				else if (ct != cq) ++n_diff;
-				s += mat[ct * 5 + cq];
-				if (s < 0) s = 0;
-				else max = max > s? max : s;
-			}
-			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
-			int n_ambi = 0;
-			for (l = 0; l < len; ++l)
-				if (qseq[qoff + l] > 3) ++n_ambi;
-			r->blen += len - n_ambi, p->n_ambi += n_ambi;
-			s -= q + e * len;
-			if (s < 0) s = 0;
-			qoff += len;
-		} else if (op == 2) { // deletion
-			int n_ambi = 0;
-			for (l = 0; l < len; ++l)
-				if (tseq[toff + l] > 3) ++n_ambi;
-			r->blen += len - n_ambi, p->n_ambi += n_ambi;
-			s -= q + e * len;
-			if (s < 0) s = 0;
-			toff += len;
-		} else if (op == 3) { // intron
-			toff += len;
-		}
-	}
-	p->dp_max = max;
-	assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
-	if (is_eqx) mm_update_cigar_eqx(r, qseq, tseq); // NB: it has to be called here as changes to qseq and tseq are not returned
-}
-
-static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) // TODO: this calls the libc realloc()
-{
-	mm_extra_t *p;
-	if (n_cigar == 0) return;
-	if (r->p == 0) {
-		uint32_t capacity = n_cigar + sizeof(mm_extra_t)/4;
-		kroundup32(capacity);
-		r->p = (mm_extra_t*)calloc(capacity, 4);
-		r->p->capacity = capacity;
-	} else if (r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4 > r->p->capacity) {
-		r->p->capacity = r->p->n_cigar + n_cigar + sizeof(mm_extra_t)/4;
-		kroundup32(r->p->capacity);
-		r->p = (mm_extra_t*)realloc(r->p, r->p->capacity * 4);
-	}
-	p = r->p;
-	if (p->n_cigar > 0 && (p->cigar[p->n_cigar-1]&0xf) == (cigar[0]&0xf)) { // same CIGAR op at the boundary
-		p->cigar[p->n_cigar-1] += cigar[0]>>4<<4;
-		if (n_cigar > 1) memcpy(p->cigar + p->n_cigar, cigar + 1, (n_cigar - 1) * 4);
-		p->n_cigar += n_cigar - 1;
-	} else {
-		memcpy(p->cigar + p->n_cigar, cigar, n_cigar * 4);
-		p->n_cigar += n_cigar;
-	}
-}
-
-static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const uint8_t *junc, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
 {
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
 		int i;
@@ -320,15 +303,12 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
 		for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr);
 		fputc('\n', stderr);
 	}
-	if (opt->max_sw_mat > 0 && (int64_t)tlen * qlen > opt->max_sw_mat) {
+	if (opt->flag & MM_F_SPLICE)
-		ksw_reset_extz(ez);
+		ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, flag, ez);
-		ez->zdropped = 1;
-	} else if (opt->flag & MM_F_SPLICE)
-		ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, opt->junc_bonus, flag, junc, ez);
 	else if (opt->q == opt->q2 && opt->e == opt->e2)
 		ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, flag, ez);
 	else
-		ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez);
+		ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ob, ez);
 	if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
 		int i;
 		fprintf(stderr, "score=%d, cigar=", ez->score);
@@ -437,7 +417,7 @@ static void mm_filter_bad_seeds_alt(void *km, int as1, int cnt1, mm128_t *a, int
 			gap2 = ((int32_t)a[as1 + j].y - (int32_t)a[as1 + j - 1].y) - (int32_t)(a[as1 + j].x - a[as1 + j - 1].x);
 			q_span_pre = a[as1 + j - 1].y >> 32 & 0xff;
 			rs2 = (int32_t)a[as1 + j - 1].x + q_span_pre;
-			qs2 = (int32_t)a[as1 + j - 1].y + q_span_pre;
+			qs2 = (int32_t)a[as1 + j - 1].x + q_span_pre;
 			m = rs2 - re1 < qs2 - qe1? rs2 - re1 : qs2 - qe1;
 			gap2 = gap2 > 0? gap2 : -gap2;
 			if (m > gap1 + gap2) break;
@@ -566,11 +546,11 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 {
 	int is_sr = !!(opt->flag & MM_F_SR), is_splice = !!(opt->flag & MM_F_SPLICE);
 	int32_t rid = a[r->as].x<<1>>33, rev = a[r->as].x>>63, as1, cnt1;
-	uint8_t *tseq, *qseq, *junc;
+	uint8_t *tseq, *qseq;
 	int32_t i, l, bw, dropped = 0, extra_flag = 0, rs0, re0, qs0, qe0;
 	int32_t rs, re, qs, qe;
 	int32_t rs1, qs1, re1, qe1;
-	int8_t mat[25];
+	int8_t mat[25], *ob;
 
 	if (is_sr) assert(!(mi->flag & MM_I_HPC)); // HPC won't work with SR because with HPC we can't easily tell if there is a gap
 
@@ -614,9 +594,11 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 		qs0 = 0, qe0 = qlen;
 		l = qs;
 		l += l * opt->a + opt->end_bonus > opt->q? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0;
+		l = l < opt->bw? l : opt->bw;
 		rs0 = rs - l > 0? rs - l : 0;
 		l = qlen - qe;
 		l += l * opt->a + opt->end_bonus > opt->q? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0;
+		l = l < opt->bw? l : opt->bw;
 		re0 = re + l < (int32_t)mi->seq[rid].len? re + l : mi->seq[rid].len;
 	} else {
 		// compute rs0 and qs0
@@ -632,7 +614,6 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				if (++l > opt->min_cnt) {
 					l = rs0 - x > qs0 - y? rs0 - x : qs0 - y;
 					rs1 = rs0 - l, qs1 = qs0 - l;
-					if (rs1 < 0) rs1 = 0; // not strictly necessary; better have this guard for explicit
 					break;
 				}
 			}
@@ -646,7 +627,6 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 			l = l < rs? l : rs;
 			rs1 = rs1 > rs - l? rs1 : rs - l;
 			rs0 = rs0 < rs1? rs0 : rs1;
-			rs0 = rs0 < rs? rs0 : rs;
 		} else rs0 = rs, qs0 = qs;
 		// compute re0 and qe0
 		re0 = (int32_t)a[r->as + r->cnt - 1].x + 1;
@@ -685,16 +665,14 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	assert(re0 > rs0);
 	tseq = (uint8_t*)kmalloc(km, re0 - rs0);
-	junc = (uint8_t*)kmalloc(km, re0 - rs0);
+	ob = (int8_t*)kmalloc(km, re0 - rs0);
 
-	if (qs > 0 && rs > 0) { // left extension; probably the condition can be changed to "qs > qs0 && rs > rs0"
+	if (qs > 0 && rs > 0) { // left extension
 		qseq = &qseq0[rev][qs0];
-		mm_idx_getseq(mi, rid, rs0, rs, tseq);
+		mm_idx_getseq2(mi, rid, rs0, rs, tseq, ob);
-		mm_idx_bed_junc(mi, rid, rs0, rs, junc);
 		mm_seq_rev(qs - qs0, qseq);
 		mm_seq_rev(rs - rs0, tseq);
-		mm_seq_rev(rs - rs0, junc);
+		mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, ob, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez);
-		mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, junc, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez);
 		if (ez->n_cigar > 0) {
 			mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			r->p->dp_score += ez->max;
@@ -719,8 +697,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				bw1 = qe - qs > re - rs? qe - qs : re - rs;
 			// perform alignment
 			qseq = &qseq0[rev][qs];
-			mm_idx_getseq(mi, rid, rs, re, tseq);
+			mm_idx_getseq2(mi, rid, rs, re, tseq, ob);
-			mm_idx_bed_junc(mi, rid, rs, re, junc);
 			if (is_sr) { // perform ungapped alignment
 				assert(qe - qs == re - rs);
 				ksw_reset_extz(ez);
@@ -730,11 +707,11 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 				}
 				ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, 0, 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
+				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, ob, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
 			}
 			// test Z-drop and inversion Z-drop
 			if ((zdrop_code = mm_test_zdrop(km, opt, qseq, tseq, ez->n_cigar, ez->cigar, mat)) != 0)
-				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, junc, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate
+				mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, ob, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate
 			// update CIGAR
 			if (ez->n_cigar > 0)
 				mm_append_cigar(r, ez->n_cigar, ez->cigar);
@@ -759,9 +736,8 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	if (!dropped && qe < qe0 && re < re0) { // right extension
 		qseq = &qseq0[rev][qe];
-		mm_idx_getseq(mi, rid, re, re0, tseq);
+		mm_idx_getseq2(mi, rid, re, re0, tseq, ob);
-		mm_idx_bed_junc(mi, rid, re, re0, junc);
+		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, ob, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
-		mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, junc, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
 		if (ez->n_cigar > 0) {
 			mm_append_cigar(r, ez->n_cigar, ez->cigar);
 			r->p->dp_score += ez->max;
@@ -777,14 +753,16 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
 
 	assert(re1 - rs1 <= re0 - rs0);
 	if (r->p) {
-		mm_idx_getseq(mi, rid, rs1, re1, tseq);
+		int left_aln = (opt->flag & MM_F_SPLICE) || mi->n_R == 0? 1 : 0;
-		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX);
+		mm_idx_getseq2(mi, rid, rs1, re1, tseq, ob);
+		mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e, left_aln);
+		if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r, &qseq0[r->rev][qs1], tseq);
 		if (rev && r->p->trans_strand)
 			r->p->trans_strand ^= 3; // flip to the read strand
 	}
 
 	kfree(km, tseq);
-	kfree(km, junc);
+	kfree(km, ob);
 }
 
 static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
@@ -837,7 +815,8 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
 	}
 	r_inv->rs = r1->re + t_off;
 	r_inv->re = r_inv->rs + ez->max_t + 1;
-	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, opt->flag & MM_F_EQX);
+	mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, 0);
+	if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r_inv, &qseq[q_off], &tseq[t_off]);
 	ret = 1;
 end_align1_inv:
 	kfree(km, tseq);

bseq.c

@@ -15,7 +15,7 @@ unsigned char seq_comp_table[256] = {
 	 48,  49,  50,	51,	 52,  53,  54,	55,	 56,  57,  58,	59,	 60,  61,  62,	63,
 	 64, 'T', 'V', 'G', 'H', 'E', 'F', 'C', 'D', 'I', 'J', 'M', 'L', 'K', 'N', 'O',
 	'P', 'Q', 'Y', 'S', 'A', 'A', 'B', 'W', 'X', 'R', 'Z',	91,	 92,  93,  94,	95,
-	 96, 't', 'v', 'g', 'h', 'e', 'f', 'c', 'd', 'i', 'j', 'm', 'l', 'k', 'n', 'o',
+	 64, 't', 'v', 'g', 'h', 'e', 'f', 'c', 'd', 'i', 'j', 'm', 'l', 'k', 'n', 'o',
 	'p', 'q', 'y', 's', 'a', 'a', 'b', 'w', 'x', 'r', 'z', 123, 124, 125, 126, 127,
 	128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
 	144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159,
@@ -39,7 +39,7 @@ mm_bseq_file_t *mm_bseq_open(const char *fn)
 {
 	mm_bseq_file_t *fp;
 	gzFile f;
-	f = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(0, "r");
+	f = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
 	if (f == 0) return 0;
 	fp = (mm_bseq_file_t*)calloc(1, sizeof(mm_bseq_file_t));
 	fp->fp = f;
@@ -65,8 +65,6 @@ static inline char *kstrdup(const kstring_t *s)
 static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual, int with_comment)
 {
 	int i;
-	if (ks->name.l == 0)
-		fprintf(stderr, "[WARNING]\033[1;31m empty sequence name in the input.\033[0m\n");
 	s->name = kstrdup(&ks->name);
 	s->seq = kstrdup(&ks->seq);
 	for (i = 0; i < (int)ks->seq.l; ++i) // convert U to T
@@ -80,7 +78,6 @@ static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual, int with_
 mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int with_comment, int frag_mode, int *n_)
 {
 	int64_t size = 0;
-	int ret;
 	kvec_t(mm_bseq1_t) a = {0,0,0};
 	kseq_t *ks = fp->ks;
 	*n_ = 0;
@@ -90,7 +87,7 @@ mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 		size = fp->s.l_seq;
 		memset(&fp->s, 0, sizeof(mm_bseq1_t));
 	}
-	while ((ret = kseq_read(ks)) >= 0) {
+	while (kseq_read(ks) >= 0) {
 		mm_bseq1_t *s;
 		assert(ks->seq.l <= INT32_MAX);
 		if (a.m == 0) kv_resize(mm_bseq1_t, 0, a, 256);
@@ -110,8 +107,6 @@ mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
 			break;
 		}
 	}
-	if (ret < -1)
-		fprintf(stderr, "[WARNING]\033[1;31m wrong FASTA/FASTQ record. Continue anyway.\033[0m\n");
 	*n_ = a.n;
 	return a.a;
 }

chain.c

@@ -14,12 +14,12 @@ static const char LogTable256[256] = {
 
 static inline int ilog2_32(uint32_t v)
 {
-	uint32_t t, tt;
+	register uint32_t t, tt;
 	if ((tt = v>>16)) return (t = tt>>8) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
 	return (t = v>>8) ? 8 + LogTable256[t] : LogTable256[v];
 }
 
-mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km)
+mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km)
 { // TODO: make sure this works when n has more than 32 bits
 	int32_t k, *f, *p, *t, *v, n_u, n_v;
 	int64_t i, j, st = 0;
@@ -28,10 +28,6 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 	mm128_t *b, *w;
 
 	if (_u) *_u = 0, *n_u_ = 0;
-	if (n == 0 || a == 0) {
-		kfree(km, a);
-		return 0;
-	}
 	f = (int32_t*)kmalloc(km, n * 4);
 	p = (int32_t*)kmalloc(km, n * 4);
 	t = (int32_t*)kmalloc(km, n * 4);
@@ -49,7 +45,6 @@ mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int m
 		int32_t max_f = q_span, n_skip = 0, min_d;
 		int32_t sidi = (a[i].y & MM_SEED_SEG_MASK) >> MM_SEED_SEG_SHIFT;
 		while (st < i && ri > a[st].x + max_dist_x) ++st;
-		if (i - st > max_iter) st = i - max_iter;
 		for (j = i - 1; j >= st; --j) {
 			int64_t dr = ri - a[j].x;
 			int32_t dq = qi - (int32_t)a[j].y, dd, sc, log_dd;

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.17/minimap2-2.17_x64-linux.tar.bz2 | tar jxf -
+curl -L https://github.com/lh3/minimap2/releases/download/v2.13/minimap2-2.13_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.17_x64-linux/{minimap2,k8,paftools.js} .  # copy executables
+cp minimap2-2.13_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

@@ -45,7 +45,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, "MIDSHN"[r->p->cigar[i]&0xf]);
 				putchar('\n');
 				free(r->p);
 			}

format.c

@@ -92,8 +92,7 @@ static void sam_write_rg_line(kstring_t *str, const char *s)
 		if (mm_verbose >= 1) fprintf(stderr, "[ERROR] the read group line contained literal <tab> characters -- replace with escaped tabs: \\t\n");
 		goto err_set_rg;
 	}
-	rg_line = (char*)malloc(strlen(s) + 1);
+	rg_line = strdup(s);
-	strcpy(rg_line, s);
 	mm_escape(rg_line);
 	if ((p = strstr(rg_line, "\tID:")) == 0) {
 		if (mm_verbose >= 1) fprintf(stderr, "[ERROR] no ID within the read group line\n");
@@ -140,8 +139,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);
+		assert(op >= 0 && op <= 3);
-		if (op == 0 || op == 7 || op == 8) { // match
+		if (op == 0) { // match
 			int l_tmp = 0;
 			for (j = 0; j < len; ++j) {
 				if (qseq[q_off + j] != tseq[t_off + j]) {
@@ -188,8 +187,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);
+		assert(op >= 0 && op <= 3);
-		if (op == 0 || op == 7 || op == 8) { // match
+		if (op == 0) { // match
 			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]]);
@@ -261,18 +260,6 @@ int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_r
 	return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0);
 }
 
-double mm_event_identity(const mm_reg1_t *r)
-{
-	int32_t i, n_gapo = 0, n_gap = 0;
-	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)
-			++n_gapo, n_gap += len;
-	}
-	return (double)r->mlen / (r->blen - n_gap + n_gapo);
-}
-
 static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 {
 	int type;
@@ -285,28 +272,20 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
 	}
 	mm_sprintf_lite(s, "\ttp:A:%c\tcm:i:%d\ts1:i:%d", type, r->cnt, r->score);
 	if (r->parent == r->id) mm_sprintf_lite(s, "\ts2:i:%d", r->subsc);
-	if (r->p) {
+	if (r->div >= 0.0f && r->div <= 1.0f) {
-		char buf[16];
+		char buf[8];
-		double div;
-		div = 1.0 - mm_event_identity(r);
-		if (div == 0.0) buf[0] = '0', buf[1] = 0;
-		else snprintf(buf, 16, "%.4f", 1.0 - mm_event_identity(r));
-		mm_sprintf_lite(s, "\tde:f:%s", buf);
-	} else if (r->div >= 0.0f && r->div <= 1.0f) {
-		char buf[16];
 		if (r->div == 0.0f) buf[0] = '0', buf[1] = 0;
-		else snprintf(buf, 16, "%.4f", r->div);
+		else sprintf(buf, "%.4f", r->div);
 		mm_sprintf_lite(s, "\tdv:f:%s", buf);
 	}
 	if (r->split) mm_sprintf_lite(s, "\tzd:i:%d", r->split);
 }
 
-void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag, int rep_len)
+void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag)
 {
 	s->l = 0;
 	if (r == 0) {
-		mm_sprintf_lite(s, "%s\t%d\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0", t->name, t->l_seq);
+		mm_sprintf_lite(s, "%s\t%d", t->name, t->l_seq);
-		if (rep_len >= 0) mm_sprintf_lite(s, "\trl:i:%d", rep_len);
 		return;
 	}
 	mm_sprintf_lite(s, "%s\t%d\t%d\t%d\t%c\t", t->name, t->l_seq, r->qs, r->qe, "+-"[r->rev]);
@@ -316,7 +295,6 @@ void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const
 	mm_sprintf_lite(s, "\t%d\t%d", r->mlen, r->blen);
 	mm_sprintf_lite(s, "\t%d", r->mapq);
 	write_tags(s, r);
-	if (rep_len >= 0) mm_sprintf_lite(s, "\trl:i:%d", rep_len);
 	if (r->p && (opt_flag & MM_F_OUT_CG)) {
 		uint32_t k;
 		mm_sprintf_lite(s, "\tcg:Z:");
@@ -329,11 +307,6 @@ void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const
 		mm_sprintf_lite(s, "\t%s", t->comment);
 }
 
-void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag)
-{
-	mm_write_paf3(s, mi, t, r, km, opt_flag, -1);
-}
-
 static void sam_write_sq(kstring_t *s, char *seq, int l, int rev, int comp)
 {
 	extern unsigned char seq_comp_table[256];
@@ -375,7 +348,6 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 			if (clip_len[1]) mm_sprintf_lite(s, ",%u", clip_len[1]<<4|clip_char);
 		} else {
 			int clip_char = (sam_flag&0x800) && !(opt_flag&MM_F_SOFTCLIP)? 'H' : 'S';
-			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]);
@@ -384,7 +356,7 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
 	}
 }
 
-void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag, int rep_len)
+void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag)
 {
 	const int max_bam_cigar_op = 65535;
 	int flag, n_regs = n_regss[seg_idx], cigar_in_tag = 0;
@@ -460,17 +432,17 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 		int tlen = 0;
 		if (this_rid >= 0 && r_next) {
 			if (this_rid == r_next->rid) {
-				if (r) {
+				int this_pos5 = r && r->rev? r->re - 1 : this_pos;
-					int this_pos5 = r->rev? r->re - 1 : this_pos;
+				int next_pos5 = r_next->rev? r_next->re - 1 : r_next->rs;
-					int next_pos5 = r_next->rev? r_next->re - 1 : r_next->rs;
+				tlen = next_pos5 - this_pos5;
-					tlen = next_pos5 - this_pos5;
-				}
 				mm_sprintf_lite(s, "\t=\t");
 			} else mm_sprintf_lite(s, "\t%s\t", mi->seq[r_next->rid].name);
 			mm_sprintf_lite(s, "%d\t", r_next->rs + 1);
 		} else if (r_next) { // && this_rid < 0
 			mm_sprintf_lite(s, "\t%s\t%d\t", mi->seq[r_next->rid].name, r_next->rs + 1);
 		} else if (this_rid >= 0) { // && r_next == NULL
+			int this_pos5 = this_rev? r->re - 1 : this_pos; // this_rev is only true when r != NULL
+			tlen = this_pos - this_pos5; // next_pos5 will be this_pos
 			mm_sprintf_lite(s, "\t=\t%d\t", this_pos + 1); // next segment will take r's coordinate
 		} else mm_sprintf_lite(s, "\t*\t0\t"); // neither has coordinates
 		if (tlen > 0) ++tlen;
@@ -535,7 +507,6 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 		if (cigar_in_tag)
 			write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag);
 	}
-	if (rep_len >= 0) mm_sprintf_lite(s, "\trl:i:%d", rep_len);
 
 	if ((opt_flag & MM_F_COPY_COMMENT) && t->comment)
 		mm_sprintf_lite(s, "\t%s", t->comment);
@@ -543,11 +514,6 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
 	s->s[s->l] = 0; // we always have room for an extra byte (see str_enlarge)
 }
 
-void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag)
-{
-	mm_write_sam3(s, mi, t, seg_idx, reg_idx, n_seg, n_regss, regss, km, opt_flag, -1);
-}
-
 void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs)
 {
 	int i;

hit.c

@@ -449,7 +449,6 @@ void mm_set_mapq(void *km, int n_regs, mm_reg1_t *regs, int min_chain_sc, int ma
 	int64_t sum_sc = 0;
 	float uniq_ratio;
 	int i;
-	if (n_regs == 0) return;
 	for (i = 0; i < n_regs; ++i)
 		if (regs[i].parent == regs[i].id)
 			sum_sc += regs[i].score;

index.c

@@ -31,16 +31,6 @@ typedef struct mm_idx_bucket_s {
 	void *h;     // hash table indexing _p_ and minimizers appearing once
 } mm_idx_bucket_t;
 
-typedef struct {
-	int32_t st, en, max; // max is not used for now
-	int32_t score:30, strand:2;
-} mm_idx_intv1_t;
-
-typedef struct mm_idx_intv_s {
-	int32_t n, m;
-	mm_idx_intv1_t *a;
-} mm_idx_intv_t;
-
 mm_idx_t *mm_idx_init(int w, int k, int b, int flag)
 {
 	mm_idx_t *mi;
@@ -65,17 +55,12 @@ void mm_idx_destroy(mm_idx_t *mi)
 			kh_destroy(idx, (idxhash_t*)mi->B[i].h);
 		}
 	}
-	if (mi->I) {
-		for (i = 0; i < mi->n_seq; ++i)
-			free(mi->I[i].a);
-		free(mi->I);
-	}
 	if (!mi->km) {
 		for (i = 0; i < mi->n_seq; ++i)
 			free(mi->seq[i].name);
 		free(mi->seq);
 	} else km_destroy(mi->km);
-	free(mi->B); free(mi->S); free(mi);
+	free(mi->R); free(mi->B); free(mi->S); free(mi);
 }
 
 const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n)
@@ -149,7 +134,7 @@ int mm_idx_name2id(const mm_idx_t *mi, const char *name)
 	return k == kh_end(h)? -1 : kh_val(h, k);
 }
 
-int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq)
+int mm_idx_getseq2(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq, int8_t *b)
 {
 	uint64_t i, st1, en1;
 	if (rid >= mi->n_seq || st >= mi->seq[rid].len) return -1;
@@ -158,9 +143,26 @@ int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, ui
 	en1 = mi->seq[rid].offset + en;
 	for (i = st1; i < en1; ++i)
 		seq[i - st1] = mm_seq4_get(mi->S, i);
+	if (b) memset(b, 0, en - st);
+	if (b && mi->R) {
+		uint32_t i, z;
+		memset(b, 0, en - st);
+		z = mm_idx_bed_query(mi, (uint64_t)rid << 32 | st);
+		for (i = z < 0? 0 : z; i < mi->n_R; ++i) {
+			uint32_t j, rr, rs, re;
+			rr = mi->R[i].x >> 32, rs = (uint32_t)mi->R[i].x, re = mi->R[i].end;
+			if (rr > rid || rs >= en) break;
+			if (rr < rid) continue;
+			re = re < en? re : en;
+			for (j = st > rs? st : rs; j < re; ++j)
+				b[j - st] = mi->R[i].score;
+		}
+	}
 	return en - st;
 }
 
+int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq) { return mm_idx_getseq2(mi, rid, st, en, seq, 0); }
+
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f)
 {
 	int i;
@@ -607,25 +609,24 @@ int mm_idx_reader_eof(const mm_idx_reader_t *r) // TODO: in extremely rare cases
 #include "kseq.h"
 KSTREAM_DECLARE(gzFile, gzread)
 
-#define sort_key_bed(a) ((a).st)
+#define sort_key_bed(a) ((a).x)
-KRADIX_SORT_INIT(bed, mm_idx_intv1_t, sort_key_bed, 4)
+KRADIX_SORT_INIT(bed, mm_idx_bed_t, sort_key_bed, 8)
 
-mm_idx_intv_t *mm_idx_read_bed(const mm_idx_t *mi, const char *fn, int read_junc)
+mm_idx_bed_t *mm_idx_bed_read_list(const mm_idx_t *mi, const char *fn, uint32_t *n_)
 {
 	gzFile fp;
 	kstream_t *ks;
 	kstring_t str = {0,0,0};
-	mm_idx_intv_t *I;
+	uint32_t n = 0, m = 0;
+	mm_idx_bed_t *r = 0;
 
 	fp = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
 	if (fp == 0) return 0;
-	I = (mm_idx_intv_t*)calloc(mi->n_seq, sizeof(*I));
 	ks = ks_init(fp);
 	while (ks_getuntil(ks, KS_SEP_LINE, &str, 0) >= 0) {
-		mm_idx_intv_t *r;
+		mm_idx_bed_t t;
-		mm_idx_intv1_t t = {-1,-1,-1,-1,0};
+		char *p, *q;
-		char *p, *q, *bl, *bs;
+		int i, id = -1, st = -1, en = -1, sc = -1;
-		int32_t i, id = -1, n_blk = 0;
 		for (p = q = str.s, i = 0;; ++p) {
 			if (*p == 0 || isspace(*p)) {
 				int32_t c = *p;
@@ -634,93 +635,56 @@ mm_idx_intv_t *mm_idx_read_bed(const mm_idx_t *mi, const char *fn, int read_junc
 					id = mm_idx_name2id(mi, q);
 					if (id < 0) break; // unknown name; TODO: throw a warning
 				} else if (i == 1) { // start
-					t.st = atol(q); // TODO: watch out integer overflow!
+					st = atoi(q);
-					if (t.st < 0) break;
+					if (st < 0) break;
 				} else if (i == 2) { // end
-					t.en = atol(q);
+					en = atoi(q);
-					if (t.en < 0) break;
+					if (en < 0) break;
+				} else if (i == 3) { // name; do nothing
 				} else if (i == 4) { // BED score
-					t.score = atol(q);
+					sc = atoi(q);
-				} else if (i == 5) { // strand
+					assert(sc >= 0 && sc <= 127);
-					t.strand = *q == '+'? 1 : *q == '-'? -1 : 0;
+				} else break;
-				} else if (i == 9) {
-					if (!isdigit(*q)) break;
-					n_blk = atol(q);
-				} else if (i == 10) {
-					bl = q;
-				} else if (i == 11) {
-					bs = q;
-					break;
-				}
 				if (c == 0) break;
 				++i, q = p + 1;
 			}
 		}
-		if (id < 0 || t.st < 0 || t.st >= t.en) continue;
+		if (en < 0) en = st + 1;
-		r = &I[id];
+		if (st < 0 || st >= en) continue;
-		if (i >= 11 && read_junc) { // BED12
+		if (m == n) EXPAND(r, m);
-			int32_t st, sz, en;
+		t.x = (uint64_t)id << 32 | st, t.end = en, t.score = sc >= 0? sc : 0, t.idx = -1;
-			st = strtol(bs, &bs, 10); ++bs;
+		r[n++] = t;
-			sz = strtol(bl, &bl, 10); ++bl;
-			en = t.st + st + sz;
-			for (i = 1; i < n_blk; ++i) {
-				mm_idx_intv1_t s = t;
-				if (r->n == r->m) {
-					r->m = r->m? r->m + (r->m>>1) : 16;
-					r->a = (mm_idx_intv1_t*)realloc(r->a, sizeof(*r->a) * r->m);
-				}
-				st = strtol(bs, &bs, 10); ++bs;
-				sz = strtol(bl, &bl, 10); ++bl;
-				s.st = en, s.en = t.st + st;
-				en = t.st + st + sz;
-				if (s.en > s.st) r->a[r->n++] = s;
-			}
-		} else {
-			if (r->n == r->m) {
-				r->m = r->m? r->m + (r->m>>1) : 16;
-				r->a = (mm_idx_intv1_t*)realloc(r->a, sizeof(*r->a) * r->m);
-			}
-			r->a[r->n++] = t;
-		}
 	}
-	free(str.s);
 	ks_destroy(ks);
 	gzclose(fp);
-	return I;
+	*n_ = n;
+	return r;
 }
 
-int mm_idx_bed_read(mm_idx_t *mi, const char *fn, int read_junc)
+int mm_idx_bed_attach(mm_idx_t *mi, uint32_t n, mm_idx_bed_t *r) // TODO: check errors
 {
-	int32_t i;
+	radix_sort_bed(r, r + n);
-	if (mi->h == 0) mm_idx_index_name(mi);
+	mi->R = r, mi->n_R = n;
-	mi->I = mm_idx_read_bed(mi, fn, read_junc);
-	if (mi->I == 0) return -1;
-	for (i = 0; i < mi->n_seq; ++i) // TODO: eliminate redundant intervals
-		radix_sort_bed(mi->I[i].a, mi->I[i].a + mi->I[i].n);
 	return 0;
 }
 
-int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uint8_t *s)
+int mm_idx_bed_read(mm_idx_t *mi, const char *fn)
 {
-	int32_t i, left, right;
+	mm_idx_bed_t *r;
-	mm_idx_intv_t *r;
+	uint32_t n;
-	memset(s, 0, en - st);
+	if (mi->h == 0) mm_idx_index_name(mi);
-	if (mi->I == 0 || ctg < 0 || ctg >= mi->n_seq) return -1;
+	r = mm_idx_bed_read_list(mi, fn, &n);
-	r = &mi->I[ctg];
+	return mm_idx_bed_attach(mi, n, r);
-	left = 0, right = r->n;
+}
-	while (right > left) {
+
+int mm_idx_bed_query(const mm_idx_t *mi, uint64_t x)
+{
+	int32_t left = -1, right = mi->n_R;
+	while (right - left > 1) {
 		int32_t mid = left + ((right - left) >> 1);
-		if (r->a[mid].st >= st) right = mid;
+		if (mi->R[mid].x > x) right = mid;
-		else left = mid + 1;
+		else if (mi->R[mid].x < x) left = mid;
-	}
+		else return mid;
-	for (i = left; i < r->n; ++i) {
-		if (st <= r->a[i].st && en >= r->a[i].en && r->a[i].strand != 0) {
-			if (r->a[i].strand > 0) {
-				s[r->a[i].st - st] |= 1, s[r->a[i].en - 1 - st] |= 2;
-			} else {
-				s[r->a[i].st - st] |= 8, s[r->a[i].en - 1 - st] |= 4;
-			}
-		}
 	}
 	return left;
 }

ketopt.h

@@ -73,17 +73,13 @@ static int ketopt(ketopt_t *s, int argc, char *argv[], int permute, const char *
 		}
 		s->opt = 0, opt = '?', s->pos = -1;
 		if (longopts) { /* parse long options */
-			int k, n_exact = 0, n_partial = 0;
+			int k, n_matches = 0;
-			const ko_longopt_t *o = 0, *o_exact = 0, *o_partial = 0;
+			const ko_longopt_t *o = 0;
 			for (j = 2; argv[s->i][j] != '\0' && argv[s->i][j] != '='; ++j) {} /* find the end of the option name */
 			for (k = 0; longopts[k].name != 0; ++k)
-				if (strncmp(&argv[s->i][2], longopts[k].name, j - 2) == 0) {
+				if (strncmp(&argv[s->i][2], longopts[k].name, j - 2) == 0)
-					if (longopts[k].name[j - 2] == 0) ++n_exact, o_exact = &longopts[k];
+					++n_matches, o = &longopts[k];
-					else ++n_partial, o_partial = &longopts[k];
+			if (n_matches == 1) {
-				}
-			if (n_exact > 1 || (n_exact == 0 && n_partial > 1)) return '?';
-			o = n_exact == 1? o_exact : n_partial == 1? o_partial : 0;
-			if (o) {
 				s->opt = opt = o->val, s->longidx = o - longopts;
 				if (argv[s->i][j] == '=') s->arg = &argv[s->i][j + 1];
 				if (o->has_arg == 1 && argv[s->i][j] == '\0') {
@@ -96,7 +92,7 @@ static int ketopt(ketopt_t *s, int argc, char *argv[], int permute, const char *
 		char *p;
 		if (s->pos == 0) s->pos = 1;
 		opt = s->opt = argv[s->i][s->pos++];
-		p = strchr((char*)ostr, opt);
+		p = strchr(ostr, opt);
 		if (p == 0) {
 			opt = '?'; /* unknown option */
 		} else if (p[1] == ':') {

ksort.h

@@ -37,7 +37,7 @@ typedef struct {
 	int depth;
 } ks_isort_stack_t;
 
-#define KSORT_SWAP(type_t, a, b) { type_t t=(a); (a)=(b); (b)=t; }
+#define KSORT_SWAP(type_t, a, b) { register type_t t=(a); (a)=(b); (b)=t; }
 
 #define KSORT_INIT(name, type_t, __sort_lt) \
 	void ks_heapdown_##name(size_t i, size_t n, type_t l[]) \

ksw2.h

@@ -58,10 +58,10 @@ void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t
 			  int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int flag, ksw_extz_t *ez);
 
 void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
+				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez);
 
 void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez);
+				   int8_t gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
 
 void ksw_extf2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t mch, int8_t mis, int8_t e, int w, int xdrop, ksw_extz_t *ez);
 

ksw2_dispatch.c

@@ -17,20 +17,18 @@
 void __cpuidex(int cpuid[4], int func_id, int subfunc_id)
 {
 #if defined(__x86_64__)
-	__asm__ volatile ("cpuid"
+	asm volatile ("cpuid"
 			: "=a" (cpuid[0]), "=b" (cpuid[1]), "=c" (cpuid[2]), "=d" (cpuid[3])
 			: "0" (func_id), "2" (subfunc_id));
 #else // on 32bit, ebx can NOT be used as PIC code
-	__asm__ volatile ("xchgl %%ebx, %1; cpuid; xchgl %%ebx, %1"
+	asm volatile ("xchgl %%ebx, %1; cpuid; xchgl %%ebx, %1"
 			: "=a" (cpuid[0]), "=r" (cpuid[1]), "=c" (cpuid[2]), "=d" (cpuid[3])
 			: "0" (func_id), "2" (subfunc_id));
 #endif
 }
 #endif
 
-static int ksw_simd = -1;
+int x86_simd(void)
-
-static int x86_simd(void)
 {
 	int flag = 0, cpuid[4], max_id;
 	__cpuidex(cpuid, 0, 0);
@@ -56,41 +54,44 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 {
 	extern void ksw_extz2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
 	extern void ksw_extz2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
-	if (ksw_simd < 0) ksw_simd = x86_simd();
+	unsigned simd;
-	if (ksw_simd & SIMD_SSE4_1)
+	simd = x86_simd();
+	if (simd & SIMD_SSE4_1)
 		ksw_extz2_sse41(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, end_bonus, flag, ez);
-	else if (ksw_simd & SIMD_SSE2)
+	else if (simd & SIMD_SSE2)
 		ksw_extz2_sse2(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, end_bonus, flag, ez);
 	else abort();
 }
 
 void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez)
 {
 	extern void ksw_extd2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez);
 	extern void ksw_extd2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez);
-	if (ksw_simd < 0) ksw_simd = x86_simd();
+	unsigned simd;
-	if (ksw_simd & SIMD_SSE4_1)
+	simd = x86_simd();
-		ksw_extd2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, ez);
+	if (simd & SIMD_SSE4_1)
-	else if (ksw_simd & SIMD_SSE2)
+		ksw_extd2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, qd, ez);
-		ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, ez);
+	else if (simd & SIMD_SSE2)
+		ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, qd, ez);
 	else abort();
 }
 
 void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
 {
 	extern void ksw_exts2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez);
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
 	extern void ksw_exts2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez);
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
-	if (ksw_simd < 0) ksw_simd = x86_simd();
+	unsigned simd;
-	if (ksw_simd & SIMD_SSE4_1)
+	simd = x86_simd();
-		ksw_exts2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, junc_bonus, flag, junc, ez);
+	if (simd & SIMD_SSE4_1)
-	else if (ksw_simd & SIMD_SSE2)
+		ksw_exts2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
-		ksw_exts2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, junc_bonus, flag, junc, ez);
+	else if (simd & SIMD_SSE2)
+		ksw_exts2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
 	else abort();
 }
 #endif

ksw2_extd2_sse.c

@@ -17,17 +17,18 @@
 #ifdef KSW_CPU_DISPATCH
 #ifdef __SSE4_1__
 void ksw_extd2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez)
 #else
 void ksw_extd2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez)
 #endif
 #else
 void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez)
 #endif // ~KSW_CPU_DISPATCH
 {
 #define __dp_code_block1 \
+	dt = _mm_load_si128(&dv[t]); \
 	z = _mm_load_si128(&s[t]); \
 	xt1 = _mm_load_si128(&x[t]);                     /* xt1 <- x[r-1][t..t+15] */ \
 	tmp = _mm_srli_si128(xt1, 15);                   /* tmp <- x[r-1][t+15] */ \
@@ -50,10 +51,10 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 #define __dp_code_block2 \
 	_mm_store_si128(&u[t], _mm_sub_epi8(z, vt1));    /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \
 	_mm_store_si128(&v[t], _mm_sub_epi8(z, ut));     /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \
-	tmp = _mm_sub_epi8(z, q_); \
+	tmp = _mm_sub_epi8(z, _mm_add_epi8(dt, q_)); \
 	a = _mm_sub_epi8(a, tmp); \
 	b = _mm_sub_epi8(b, tmp); \
-	tmp = _mm_sub_epi8(z, q2_); \
+	tmp = _mm_sub_epi8(z, _mm_add_epi8(dt, q2_)); \
 	a2= _mm_sub_epi8(a2, tmp); \
 	b2= _mm_sub_epi8(b2, tmp);
 
@@ -61,8 +62,8 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX);
 	int32_t *H = 0, H0 = 0, last_H0_t = 0;
 	uint8_t *qr, *sf, *mem, *mem2 = 0;
-	__m128i q_, q2_, qe_, qe2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_;
+	__m128i q_, q2_, qe_, qe2_, e_, e2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_;
-	__m128i *u, *v, *x, *y, *x2, *y2, *s, *p = 0;
+	__m128i *u, *v, *x, *y, *x2, *y2, *s, *p = 0, *dv;
 
 	ksw_reset_extz(ez);
 	if (m <= 1 || qlen <= 0 || tlen <= 0) return;
@@ -72,6 +73,8 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	zero_   = _mm_set1_epi8(0);
 	q_      = _mm_set1_epi8(q);
 	q2_     = _mm_set1_epi8(q2);
+	e_      = _mm_set1_epi8(e);
+	e2_     = _mm_set1_epi8(e2);
 	qe_     = _mm_set1_epi8(q + e);
 	qe2_    = _mm_set1_epi8(q2 + e2);
 	sc_mch_ = _mm_set1_epi8(mat[0]);
@@ -96,16 +99,23 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		++long_thres;
 	long_diff = long_thres * (e - e2) - (q2 - q) - e2;
 
-	mem = (uint8_t*)kcalloc(km, tlen_ * 8 + qlen_ + 1, 16);
+	mem = (uint8_t*)kcalloc(km, tlen_ * 9 + qlen_ + 1, 16);
 	u = (__m128i*)(((size_t)mem + 15) >> 4 << 4); // 16-byte aligned
-	v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_;
+	v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_, dv = y2 + tlen_;
-	s = y2 + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16;
+	s = dv + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16;
 	memset(u,  -q  - e,  tlen_ * 16);
 	memset(v,  -q  - e,  tlen_ * 16);
 	memset(x,  -q  - e,  tlen_ * 16);
 	memset(y,  -q  - e,  tlen_ * 16);
 	memset(x2, -q2 - e2, tlen_ * 16);
 	memset(y2, -q2 - e2, tlen_ * 16);
+	if (qd) {
+		int8_t *tmp = (int8_t*)dv;
+		for (t = 0; t < tlen; ++t) tmp[t] = -qd[t];
+		fprintf(stderr, "%d\t%d\t%x\n", tlen, qlen, flag&KSW_EZ_RIGHT);
+		for (t = 0; t < tlen; ++t) fputc("ACGTN"[target[t]], stderr); fputc('\n', stderr);
+		for (t = 0; t < tlen; ++t) fputc('0' + qd[t], stderr); fputc('\n', stderr);
+	}
 	if (!approx_max) {
 		H = (int32_t*)kmalloc(km, tlen_ * 16 * 4);
 		for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF;
@@ -182,7 +192,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		assert(en_ - st_ + 1 <= n_col_);
 		if (!with_cigar) { // score only
 			for (t = st_; t <= en_; ++t) {
-				__m128i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp;
+				__m128i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt;
 				__dp_code_block1;
 #ifdef __SSE4_1__
 				z = _mm_max_epi8(z, a);
@@ -191,10 +201,10 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				z = _mm_max_epi8(z, b2);
 				z = _mm_min_epi8(z, sc_mch_);
 				__dp_code_block2; // save u[] and v[]; update a, b, a2 and b2
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_max_epi8(a,  zero_), qe_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_max_epi8(a,  zero_), _mm_add_epi8(dt, qe_)));
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_max_epi8(b,  zero_), qe_));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_max_epi8(b,  zero_), _mm_add_epi8(dt, qe_)));
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_max_epi8(a2, zero_), qe2_));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_max_epi8(a2, zero_), _mm_add_epi8(dt, qe2_)));
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_max_epi8(b2, zero_), qe2_));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_max_epi8(b2, zero_), _mm_add_epi8(dt, qe2_)));
 #else
 				tmp = _mm_cmpgt_epi8(a,  z);
 				z = _mm_or_si128(_mm_andnot_si128(tmp, z), _mm_and_si128(tmp, a));
@@ -208,20 +218,20 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 				z = _mm_or_si128(_mm_and_si128(tmp, sc_mch_), _mm_andnot_si128(tmp, z));
 				__dp_code_block2;
 				tmp = _mm_cmpgt_epi8(a, zero_);
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  qe_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  _mm_add_epi8(dt, qe_)));
 				tmp = _mm_cmpgt_epi8(b, zero_);
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  qe_));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  _mm_add_epi8(dt, qe_)));
 				tmp = _mm_cmpgt_epi8(a2, zero_);
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), qe2_));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), _mm_add_epi8(dt, qe2_)));
 				tmp = _mm_cmpgt_epi8(b2, zero_);
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), qe2_));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), _mm_add_epi8(dt, qe2_)));
 #endif
 			}
 		} else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment
 			__m128i *pr = p + (size_t)r * n_col_ - st_;
 			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
-				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp;
+				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt;
 				__dp_code_block1;
 #ifdef __SSE4_1__
 				d = _mm_and_si128(_mm_cmpgt_epi8(a, z), _mm_set1_epi8(1));       // d = a  > z? 1 : 0
@@ -251,16 +261,16 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 #endif
 				__dp_code_block2;
 				tmp = _mm_cmpgt_epi8(a, zero_);
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  qe_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_and_si128(tmp, a),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x08))); // d = a > 0? 1<<3 : 0
 				tmp = _mm_cmpgt_epi8(b, zero_);
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  qe_));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_and_si128(tmp, b),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x10))); // d = b > 0? 1<<4 : 0
 				tmp = _mm_cmpgt_epi8(a2, zero_);
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), qe2_));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x20))); // d = a > 0? 1<<5 : 0
 				tmp = _mm_cmpgt_epi8(b2, zero_);
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), qe2_));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x40))); // d = b > 0? 1<<6 : 0
 				_mm_store_si128(&pr[t], d);
 			}
@@ -268,7 +278,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 			__m128i *pr = p + (size_t)r * n_col_ - st_;
 			off[r] = st, off_end[r] = en;
 			for (t = st_; t <= en_; ++t) {
-				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp;
+				__m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt;
 				__dp_code_block1;
 #ifdef __SSE4_1__
 				d = _mm_andnot_si128(_mm_cmpgt_epi8(z, a), _mm_set1_epi8(1));    // d = z > a?  0 : 1
@@ -298,16 +308,16 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 #endif
 				__dp_code_block2;
 				tmp = _mm_cmpgt_epi8(zero_, a);
-				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, a),  qe_));
+				_mm_store_si128(&x[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, a),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x08))); // d = a > 0? 1<<3 : 0
 				tmp = _mm_cmpgt_epi8(zero_, b);
-				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, b),  qe_));
+				_mm_store_si128(&y[t],  _mm_sub_epi8(_mm_andnot_si128(tmp, b),  _mm_add_epi8(dt, qe_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x10))); // d = b > 0? 1<<4 : 0
 				tmp = _mm_cmpgt_epi8(zero_, a2);
-				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, a2), qe2_));
+				_mm_store_si128(&x2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, a2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x20))); // d = a > 0? 1<<5 : 0
 				tmp = _mm_cmpgt_epi8(zero_, b2);
-				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, b2), qe2_));
+				_mm_store_si128(&y2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, b2), _mm_add_epi8(dt, qe2_)));
 				d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x40))); // d = b > 0? 1<<6 : 0
 				_mm_store_si128(&pr[t], d);
 			}
@@ -376,6 +386,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 		last_st = st, last_en = en;
 		//for (t = st0; t <= en0; ++t) printf("(%d,%d)\t(%d,%d,%d,%d)\t%d\n", r, t, ((int8_t*)u)[t], ((int8_t*)v)[t], ((int8_t*)x)[t], ((int8_t*)y)[t], H[t]); // for debugging
 	}
+	fprintf(stderr, "score: %d\n", ez->score);
 	kfree(km, mem);
 	if (!approx_max) kfree(km, H);
 	if (with_cigar) { // backtrack

ksw2_exts2_sse.c

@@ -17,14 +17,14 @@
 #ifdef KSW_CPU_DISPATCH
 #ifdef __SSE4_1__
 void ksw_exts2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
 #else
 void ksw_exts2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
 #endif
 #else
 void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
-				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int8_t junc_bonus, int flag, const uint8_t *junc, ksw_extz_t *ez)
+				   int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez)
 #endif // ~KSW_CPU_DISPATCH
 {
 #define __dp_code_block1 \
@@ -113,53 +113,20 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
 	if (flag & (KSW_EZ_SPLICE_FOR|KSW_EZ_SPLICE_REV)) {
 		int semi_cost = flag&KSW_EZ_SPLICE_FLANK? -noncan/2 : 0; // GTr or yAG is worth 0.5 bit; see PMID:18688272
 		memset(donor, -noncan, tlen_ * 16);
+		for (t = 0; t < tlen - 4; ++t) {
+			int can_type = 0; // type of canonical site: 0=none, 1=GT/AG only, 2=GTr/yAG
+			if ((flag & KSW_EZ_SPLICE_FOR) && target[t+1] == 2 && target[t+2] == 3) can_type = 1; // GTr...
+			if ((flag & KSW_EZ_SPLICE_REV) && target[t+1] == 1 && target[t+2] == 3) can_type = 1; // CTr...
+			if (can_type && (target[t+3] == 0 || target[t+3] == 2)) can_type = 2;
+			if (can_type) ((int8_t*)donor)[t] = can_type == 2? 0 : semi_cost;
+		}
 		memset(acceptor, -noncan, tlen_ * 16);
-		if (!(flag & KSW_EZ_REV_CIGAR)) {
+		for (t = 2; t < tlen; ++t) {
-			for (t = 0; t < tlen - 4; ++t) {
+			int can_type = 0;
-				int can_type = 0; // type of canonical site: 0=none, 1=GT/AG only, 2=GTr/yAG
+			if ((flag & KSW_EZ_SPLICE_FOR) && target[t-1] == 0 && target[t] == 2) can_type = 1; // ...yAG
-				if ((flag & KSW_EZ_SPLICE_FOR) && target[t+1] == 2 && target[t+2] == 3) can_type = 1; // GTr...
+			if ((flag & KSW_EZ_SPLICE_REV) && target[t-1] == 0 && target[t] == 1) can_type = 1; // ...yAC
-				if ((flag & KSW_EZ_SPLICE_REV) && target[t+1] == 1 && target[t+2] == 3) can_type = 1; // CTr...
+			if (can_type && (target[t-2] == 1 || target[t-2] == 3)) can_type = 2;
-				if (can_type && (target[t+3] == 0 || target[t+3] == 2)) can_type = 2;
+			if (can_type) ((int8_t*)acceptor)[t] = can_type == 2? 0 : semi_cost;
-				if (can_type) ((int8_t*)donor)[t] = can_type == 2? 0 : semi_cost;
-			}
-			if (junc)
-				for (t = 0; t < tlen - 1; ++t)
-					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t+1]&1)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t+1]&8)))
-						((int8_t*)donor)[t] += junc_bonus;
-			for (t = 2; t < tlen; ++t) {
-				int can_type = 0;
-				if ((flag & KSW_EZ_SPLICE_FOR) && target[t-1] == 0 && target[t] == 2) can_type = 1; // ...yAG
-				if ((flag & KSW_EZ_SPLICE_REV) && target[t-1] == 0 && target[t] == 1) can_type = 1; // ...yAC
-				if (can_type && (target[t-2] == 1 || target[t-2] == 3)) can_type = 2;
-				if (can_type) ((int8_t*)acceptor)[t] = can_type == 2? 0 : semi_cost;
-			}
-			if (junc)
-				for (t = 0; t < tlen; ++t)
-					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t]&2)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t]&4)))
-						((int8_t*)acceptor)[t] += junc_bonus;
-		} else {
-			for (t = 0; t < tlen - 4; ++t) {
-				int can_type = 0; // type of canonical site: 0=none, 1=GT/AG only, 2=GTr/yAG
-				if ((flag & KSW_EZ_SPLICE_FOR) && target[t+1] == 2 && target[t+2] == 0) can_type = 1; // GAy...
-				if ((flag & KSW_EZ_SPLICE_REV) && target[t+1] == 1 && target[t+2] == 0) can_type = 1; // CAy...
-				if (can_type && (target[t+3] == 1 || target[t+3] == 3)) can_type = 2;
-				if (can_type) ((int8_t*)donor)[t] = can_type == 2? 0 : semi_cost;
-			}
-			if (junc)
-				for (t = 0; t < tlen - 1; ++t)
-					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t+1]&2)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t+1]&4)))
-						((int8_t*)donor)[t] += junc_bonus;
-			for (t = 2; t < tlen; ++t) {
-				int can_type = 0;
-				if ((flag & KSW_EZ_SPLICE_FOR) && target[t-1] == 3 && target[t] == 2) can_type = 1; // ...rTG
-				if ((flag & KSW_EZ_SPLICE_REV) && target[t-1] == 3 && target[t] == 1) can_type = 1; // ...rTC
-				if (can_type && (target[t-2] == 0 || target[t-2] == 2)) can_type = 2;
-				if (can_type) ((int8_t*)acceptor)[t] = can_type == 2? 0 : semi_cost;
-			}
-			if (junc)
-				for (t = 0; t < tlen; ++t)
-					if (((flag & KSW_EZ_SPLICE_FOR) && (junc[t]&1)) || ((flag & KSW_EZ_SPLICE_REV) && (junc[t]&8)))
-						((int8_t*)acceptor)[t] += junc_bonus;
 		}
 	}
 

main.c

@@ -6,7 +6,7 @@
 #include "mmpriv.h"
 #include "ketopt.h"
 
-#define MM_VERSION "2.17-r941"
+#define MM_VERSION "2.13-r852-dirty"
 
 #ifdef __linux__
 #include <sys/resource.h>
@@ -60,12 +60,7 @@ static ko_longopt_t long_options[] = {
 	{ "split-prefix",   ko_required_argument, 334 },
 	{ "no-end-flt",     ko_no_argument,       335 },
 	{ "hard-mask-level",ko_no_argument,       336 },
-	{ "cap-sw-mem",     ko_required_argument, 337 },
+	{ "bed",            ko_required_argument, 337 },
-	{ "max-qlen",       ko_required_argument, 338 },
-	{ "max-chain-iter", ko_required_argument, 339 },
-	{ "junc-bed",       ko_required_argument, 340 },
-	{ "junc-bonus",     ko_required_argument, 341 },
-	{ "sam-hit-only",   ko_no_argument,       342 },
 	{ "help",           ko_no_argument,       'h' },
 	{ "max-intron-len", ko_required_argument, 'G' },
 	{ "version",        ko_no_argument,       'V' },
@@ -103,12 +98,12 @@ static inline void yes_or_no(mm_mapopt_t *opt, int flag, int long_idx, const cha
 
 int main(int argc, char *argv[])
 {
-	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:yYPo:";
+	const char *opt_str = "2aSDw:k:K:t:r:f:Vv:g:G:I:d:XT:s:x:Hcp:M:n:z:A:B:O:E:m:N:Qu:R:hF:LC:yYP";
 	ketopt_t o = KETOPT_INIT;
 	mm_mapopt_t opt;
 	mm_idxopt_t ipt;
 	int i, c, n_threads = 3, n_parts, old_best_n = -1;
-	char *fnw = 0, *rg = 0, *junc_bed = 0, *s;
+	char *fnw = 0, *fn_bed = 0, *rg = 0, *s;
 	FILE *fp_help = stderr;
 	mm_idx_reader_t *idx_rdr;
 	mm_idx_t *mi;
@@ -128,7 +123,7 @@ int main(int argc, char *argv[])
 			fprintf(stderr, "[ERROR] missing option argument\n");
 			return 1;
 		} else if (c == '?') {
-			fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[o.i - 1]);
+			fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[o.i]);
 			return 1;
 		}
 	}
@@ -169,21 +164,12 @@ int main(int argc, char *argv[])
 		else if (c == 'R') rg = o.arg;
 		else if (c == 'h') fp_help = stdout;
 		else if (c == '2') opt.flag |= MM_F_2_IO_THREADS;
-		else if (c == 'o') {
-			if (strcmp(o.arg, "-") != 0) {
-				if (freopen(o.arg, "wb", stdout) == NULL) {
-					fprintf(stderr, "[ERROR]\033[1;31m failed to write the output to file '%s'\033[0m\n", o.arg);
-					exit(1);
-				}
-			}
-		}
 		else if (c == 300) ipt.bucket_bits = atoi(o.arg); // --bucket-bits
 		else if (c == 302) opt.seed = atoi(o.arg); // --seed
 		else if (c == 303) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
 		else if (c == 304) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
 		else if (c == 306) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
 		else if (c == 307) opt.max_chain_skip = atoi(o.arg); // --max-chain-skip
-		else if (c == 339) opt.max_chain_iter = atoi(o.arg); // --max-chain-iter
 		else if (c == 308) opt.min_ksw_len = atoi(o.arg); // --min-dp-len
 		else if (c == 309) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
 		else if (c == 310) opt.flag |= MM_F_SPLICE; // --splice
@@ -205,10 +191,7 @@ int main(int argc, char *argv[])
 		else if (c == 334) opt.split_prefix = o.arg; // --split-prefix
 		else if (c == 335) opt.flag |= MM_F_NO_END_FLT; // --no-end-flt
 		else if (c == 336) opt.flag |= MM_F_HARD_MLEVEL; // --hard-mask-level
-		else if (c == 337) opt.max_sw_mat = mm_parse_num(o.arg); // --cap-sw-mat
+		else if (c == 337) fn_bed = o.arg; // --bed-prefer
-		else if (c == 338) opt.max_qlen = mm_parse_num(o.arg); // --max-qlen
-		else if (c == 340) junc_bed = o.arg; // --junc-bed
-		else if (c == 342) opt.flag |= MM_F_SAM_HIT_ONLY; // --sam-hit-only
 		else if (c == 314) { // --frag
 			yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
 		} else if (c == 315) { // --secondary
@@ -283,7 +266,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "  Indexing:\n");
 		fprintf(fp_help, "    -H           use homopolymer-compressed k-mer (preferrable for PacBio)\n");
 		fprintf(fp_help, "    -k INT       k-mer size (no larger than 28) [%d]\n", ipt.k);
-		fprintf(fp_help, "    -w INT       minimizer window size [%d]\n", ipt.w);
+		fprintf(fp_help, "    -w INT       minizer window size [%d]\n", ipt.w);
 		fprintf(fp_help, "    -I NUM       split index for every ~NUM input bases [4G]\n");
 		fprintf(fp_help, "    -d FILE      dump index to FILE []\n");
 		fprintf(fp_help, "  Mapping:\n");
@@ -308,7 +291,7 @@ int main(int argc, char *argv[])
 		fprintf(fp_help, "    -u CHAR      how to find GT-AG. f:transcript strand, b:both strands, n:don't match GT-AG [n]\n");
 		fprintf(fp_help, "  Input/Output:\n");
 		fprintf(fp_help, "    -a           output in the SAM format (PAF by default)\n");
-		fprintf(fp_help, "    -o FILE      output alignments to FILE [stdout]\n");
+		fprintf(fp_help, "    -Q           don't output base quality in SAM\n");
 		fprintf(fp_help, "    -L           write CIGAR with >65535 ops at the CG tag\n");
 		fprintf(fp_help, "    -R STR       SAM read group line in a format like '@RG\\tID:foo\\tSM:bar' []\n");
 		fprintf(fp_help, "    -c           output CIGAR in PAF\n");
@@ -367,8 +350,8 @@ int main(int argc, char *argv[])
 			fprintf(stderr, "[M::%s::%.3f*%.2f] loaded/built the index for %d target sequence(s)\n",
 					__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), mi->n_seq);
 		if (argc != o.ind + 1) mm_mapopt_update(&opt, mi);
+		if (fn_bed) mm_idx_bed_read(mi, fn_bed);
 		if (mm_verbose >= 3) mm_idx_stat(mi);
-		if (junc_bed) mm_idx_bed_read(mi, junc_bed, 1);
 		if (!(opt.flag & MM_F_FRAG_MODE)) {
 			for (i = o.ind + 1; i < argc; ++i)
 				mm_map_file(mi, argv[i], &opt, n_threads);

map.c

@@ -284,7 +284,6 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		qlen_sum += qlens[i], n_regs[i] = 0, regs[i] = 0;
 
 	if (qlen_sum == 0 || n_segs <= 0 || n_segs > MM_MAX_SEG) return;
-	if (opt->max_qlen > 0 && qlen_sum > opt->max_qlen) return;
 
 	hash  = qname? __ac_X31_hash_string(qname) : 0;
 	hash ^= __ac_Wang_hash(qlen_sum) + __ac_Wang_hash(opt->seed);
@@ -312,7 +311,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 		if (max_chain_gap_ref < opt->max_gap) max_chain_gap_ref = opt->max_gap;
 	} else max_chain_gap_ref = opt->max_gap;
 
-	a = mm_chain_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, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
+	a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 
 	if (opt->max_occ > opt->mid_occ && rep_len > 0) {
 		int rechain = 0;
@@ -334,7 +333,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char **
 			kfree(b->km, mini_pos);
 			if (opt->flag & MM_F_HEAP_SORT) a = collect_seed_hits_heap(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
 			else a = collect_seed_hits(b->km, opt, opt->max_occ, mi, qname, &mv, qlen_sum, &n_a, &rep_len, &n_mini_pos, &mini_pos);
-			a = mm_chain_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, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
+			a = mm_chain_dp(max_chain_gap_ref, max_chain_gap_qry, opt->bw, opt->max_chain_skip, opt->min_cnt, opt->min_chain_score, is_splice, n_segs, n_a, a, &n_regs0, &u, b->km);
 		}
 	}
 	b->frag_gap = max_chain_gap_ref;
@@ -584,16 +583,16 @@ static void *worker_pipeline(void *shared, int step, void *in)
 						if ((p->opt->flag & MM_F_NO_PRINT_2ND) && r->id != r->parent)
 							continue;
 						if (p->opt->flag & MM_F_OUT_SAM)
-							mm_write_sam3(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag, s->rep_len[i]);
+							mm_write_sam2(&p->str, mi, t, i - seg_st, j, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
 						else
-							mm_write_paf3(&p->str, mi, t, r, km, p->opt->flag, s->rep_len[i]);
+							mm_write_paf(&p->str, mi, t, r, km, p->opt->flag);
 						mm_err_puts(p->str.s);
 					}
-				} else if ((p->opt->flag & MM_F_PAF_NO_HIT) || ((p->opt->flag & MM_F_OUT_SAM) && !(p->opt->flag & MM_F_SAM_HIT_ONLY))) { // output an empty hit, if requested
+				} else if (p->opt->flag & (MM_F_OUT_SAM|MM_F_PAF_NO_HIT)) { // output an empty hit, if requested
 					if (p->opt->flag & MM_F_OUT_SAM)
-						mm_write_sam3(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag, s->rep_len[i]);
+						mm_write_sam2(&p->str, mi, t, i - seg_st, -1, s->n_seg[k], &s->n_reg[seg_st], (const mm_reg1_t*const*)&s->reg[seg_st], km, p->opt->flag);
 					else
-						mm_write_paf3(&p->str, mi, t, 0, 0, p->opt->flag, s->rep_len[i]);
+						mm_write_paf(&p->str, mi, t, 0, 0, p->opt->flag);
 					mm_err_puts(p->str.s);
 				}
 			}

minimap.h

@@ -35,7 +35,6 @@
 #define MM_F_PAF_NO_HIT    0x8000000 // output unmapped reads to PAF
 #define MM_F_NO_END_FLT    0x10000000
 #define MM_F_HARD_MLEVEL   0x20000000
-#define MM_F_SAM_HIT_ONLY  0x40000000
 
 #define MM_I_HPC          0x1
 #define MM_I_NO_SEQ       0x2
@@ -60,14 +59,21 @@ typedef struct {
 	uint32_t len;    // length
 } mm_idx_seq_t;
 
+typedef struct {
+	uint64_t x;
+	int32_t end, idx;
+	int32_t score; // NB: wasting 4 bytes due to memory alignment
+} mm_idx_bed_t;
+
 typedef struct {
 	int32_t b, w, k, flag;
 	uint32_t n_seq;            // number of reference sequences
 	int32_t index;
+	uint32_t n_R;
 	mm_idx_seq_t *seq;         // sequence name, length and offset
 	uint32_t *S;               // 4-bit packed sequence
 	struct mm_idx_bucket_s *B; // index (hidden)
-	struct mm_idx_intv_s *I;   // intervals (hidden)
+	mm_idx_bed_t *R;
 	void *km, *h;
 } mm_idx_t;
 
@@ -105,16 +111,14 @@ typedef struct {
 } mm_idxopt_t;
 
 typedef struct {
-	int64_t flag;    // see MM_F_* macros
 	int seed;
 	int sdust_thres; // score threshold for SDUST; 0 to disable
-
+	int flag;        // see MM_F_* macros
-	int max_qlen;    // max query length
 
 	int bw;          // bandwidth
 	int max_gap, max_gap_ref; // break a chain if there are no minimizers in a max_gap window
 	int max_frag_len;
-	int max_chain_skip, max_chain_iter;
+	int max_chain_skip;
 	int min_cnt;         // min number of minimizers on each chain
 	int min_chain_score; // min chaining score
 
@@ -129,7 +133,6 @@ typedef struct {
 	int a, b, q, e, q2, e2; // matching score, mismatch, gap-open and gap-ext penalties
 	int sc_ambi; // score when one or both bases are "N"
 	int noncan;      // cost of non-canonical splicing sites
-	int junc_bonus;
 	int zdrop, zdrop_inv;   // break alignment if alignment score drops too fast along the diagonal
 	int end_bonus;
 	int min_dp_max;  // drop an alignment if the score of the max scoring segment is below this threshold
@@ -144,7 +147,6 @@ typedef struct {
 	int32_t mid_occ;     // ignore seeds with occurrences above this threshold
 	int32_t max_occ;
 	int mini_batch_size; // size of a batch of query bases to process in parallel
-	int64_t max_sw_mat;
 
 	const char *split_prefix;
 } mm_mapopt_t;
@@ -368,8 +370,10 @@ int mm_idx_index_name(mm_idx_t *mi);
 int mm_idx_name2id(const mm_idx_t *mi, const char *name);
 int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq);
 
-int mm_idx_bed_read(mm_idx_t *mi, const char *fn, int read_junc);
+// BED operations
-int mm_idx_bed_junc(const mm_idx_t *mi, int32_t ctg, int32_t st, int32_t en, uint8_t *s);
+int mm_idx_bed_read(mm_idx_t *mi, const char *fn);
+int mm_idx_bed_attach(mm_idx_t *mi, uint32_t n, mm_idx_bed_t *r);
+int mm_idx_bed_query(const mm_idx_t *mi, uint64_t x);
 
 // deprecated APIs for backward compatibility
 void mm_mapopt_init(mm_mapopt_t *opt);

minimap2.1

@@ -1,4 +1,4 @@
-.TH minimap2 1 "4 May 2019" "minimap2-2.17 (r941)" "Bioinformatics tools"
+.TH minimap2 1 "11 October 2018" "minimap2-2.13 (r850)" "Bioinformatics tools"
 .SH NAME
 .PP
 minimap2 - mapping and alignment between collections of DNA sequences
@@ -232,19 +232,13 @@ Honor option
 and disable a heurstic to save unmapped subsequences.
 .TP
 .BI --max-chain-skip \ INT
-A heuristics that stops chaining early [25]. Minimap2 uses dynamic programming
+A heuristics that stops chaining early [50]. Minimap2 uses dynamic programming
 for chaining. The time complexity is quadratic in the number of seeds. This
 option makes minimap2 exits the inner loop if it repeatedly sees seeds already
 on chains. Set
 .I INT
 to a large number to switch off this heurstics.
 .TP
-.BI --max-chain-iter \ INT
-Check up to
-.I INT
-partial chains during chaining [5000]. This is a heuristic to avoid quadratic
-time complexity in the worst case.
-.TP
 .B --no-long-join
 Disable the long gap patching heuristic. When this option is applied, the
 maximum alignment gap is mostly controlled by
@@ -280,10 +274,6 @@ Only map to the reverse complement strand of the reference sequences.
 .BR --heap-sort = no | yes
 If yes, sort anchors with heap merge, instead of radix sort. Heap merge is
 faster for short reads, but slower for long reads. [no]
-.TP
-.B --no-pairing
-Treat two reads in a pair as independent reads. The mate related fields in SAM
-are still properly populated.
 .SS Alignment options
 .TP 10
 .BI -A \ INT
@@ -364,17 +354,6 @@ on SIRV data, please add
 .B --splice-flank=no
 to the command line.
 .TP
-.BR --junc-bed \ FILE
-Gene annotations in the BED12 format (aka 12-column BED), or intron positions
-in 5-column BED. With this option, minimap2 prefers splicing in annotations.
-BED12 file can be converted from GTF/GFF3 with `paftools.js gff2bed anno.gtf'
-[].
-.TP
-.BR --junc-bonus \ INT
-Score bonus for a splice donor or acceptor found in annotation (effective with
-.BR --junc-bed )
-[0].
-.TP
 .BI --end-seed-pen \ INT
 Drop a terminal anchor if
 .IR s <log( g )+ INT ,
@@ -390,22 +369,12 @@ It helps to avoid tiny terminal exons. [6]
 .B --no-end-flt
 Don't filter seeds towards the ends of chains before performing base-level
 alignment.
-.TP
-.BI --cap-sw-mem \ NUM
-Skip alignment if the DP matrix size is above
-.IR NUM .
-Set 0 to disable [0].
 .SS Input/output options
 .TP 10
 .B -a
 Generate CIGAR and output alignments in the SAM format. Minimap2 outputs in PAF
 by default.
 .TP
-.BI -o \ FILE
-Output alignments to
-.I FILE
-[stdout].
-.TP
 .B -Q
 Ignore base quality in the input file.
 .TP
@@ -480,18 +449,6 @@ memory.
 .BR --secondary = yes | no
 Whether to output secondary alignments [yes]
 .TP
-.BI --max-qlen \ NUM
-Filter out query sequences longer than
-.IR NUM .
-.TP
-.B --paf-no-hit
-In PAF, output unmapped queries; the strand and the reference name fields are
-set to `*'. Warning: some paftools.js commands may not work with such output
-for the moment.
-.TP
-.B --sam-hit-only
-In SAM, don't output unmapped reads.
-.TP
 .B --version
 Print version number to stdout
 .SS Preset options
@@ -521,7 +478,7 @@ is determined by the sequencing error mode.
 .B asm5
 Long assembly to reference mapping
 .RB ( -k19
-.B -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200 -N50
+.B -w19 -A1 -B19 -O39,81 -E3,1 -s200 -z200
 .BR --min-occ-floor=100 ).
 Typically, the alignment will not extend to regions with 5% or higher sequence
 divergence. Only use this preset if the average divergence is far below 5%.
@@ -529,14 +486,14 @@ divergence. Only use this preset if the average divergence is far below 5%.
 .B asm10
 Long assembly to reference mapping
 .RB ( -k19
-.B -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200 -N50
+.B -w19 -A1 -B9 -O16,41 -E2,1 -s200 -z200
 .BR --min-occ-floor=100 ).
 Up to 10% sequence divergence.
 .TP
 .B asm20
 Long assembly to reference mapping
 .RB ( -k19
-.B -w10 -A1 -B4 -O6,26 -E2,1 -s200 -z200 -N50
+.B -w10 -A1 -B6 -O6,26 -E2,1 -s200 -z200
 .BR --min-occ-floor=100 ).
 Up to 20% sequence divergence.
 .TP
@@ -558,7 +515,7 @@ is that this preset is not using HPC minimizers.
 .B splice
 Long-read spliced alignment
 .RB ( -k15
-.B -w5 --splice -g2000 -G200k -A1 -B2 -O2,32 -E1,0 -C9 -z200 -ub --junc-bonus=9
+.B -w5 --splice -g2000 -G200k -A1 -B2 -O2,32 -E1,0 -C9 -z200 -ub
 .BR --splice-flank=yes ).
 In the splice mode, 1) long deletions are taken as introns and represented as
 the
@@ -568,12 +525,6 @@ costs are different during chaining; 4) the computation of the
 .RB ` ms '
 tag ignores introns to demote hits to pseudogenes.
 .TP
-.B splice:hq
-Long-read splice alignment for PacBio CCS reads
-.RB ( -xsplice
-.B -C5 -O6,24
-.BR -B4 ).
-.TP
 .B sr
 Short single-end reads without splicing
 .RB ( -k21
@@ -644,8 +595,6 @@ ts	A	Transcript strand (splice mode only)
 cg	Z	CIGAR string (only in PAF)
 cs	Z	Difference string
 dv	f	Approximate per-base sequence divergence
-de	f	Gap-compressed per-base sequence divergence
-rl	i	Length of query regions harboring repetitive seeds
 .TE
 
 .PP

misc.c

@@ -116,7 +116,8 @@ long peakrss(void)
 double realtime(void)
 {
 	struct timeval tp;
-	gettimeofday(&tp, NULL);
+	struct timezone tzp;
+	gettimeofday(&tp, &tzp);
 	return tp.tv_sec + tp.tv_usec * 1e-6;
 }
 

misc/paftools.js

@@ -1,6 +1,6 @@
 #!/usr/bin/env k8
 
-var paftools_version = '2.17-r941';
+var paftools_version = '2.13-r850';
 
 /*****************************
  ***** Library functions *****
@@ -433,7 +433,6 @@ function paf_call(args)
 	while (file.readline(buf) >= 0) {
 		var line = buf.toString();
 		var m, t = line.split("\t", 12);
-		if (t.length < 12 || t[5] == '*') continue; // unmapped
 		for (var i = 6; i <= 11; ++i)
 			t[i] = parseInt(t[i]);
 		if (t[10] < min_cov_len || t[11] < min_mapq) continue;
@@ -565,18 +564,14 @@ function paf_call(args)
 
 function paf_asmstat(args)
 {
-	var c, min_query_len = 0, min_seg_len = 10000, max_diff = 0.01, bp_flank_len = 0, bp_gap_len = 0;
+	var c, min_seg_len = 10000, max_diff = 0.01;
-	while ((c = getopt(args, "l:d:b:g:q:")) != null) {
+	while ((c = getopt(args, "l:d:")) != null) {
 		if (c == 'l') min_seg_len = parseInt(getopt.arg);
 		else if (c == 'd') max_diff = parseFloat(getopt.arg);
-		else if (c == 'b') bp_flank_len = parseInt(getopt.arg);
-		else if (c == 'g') bp_gap_len = parseInt(getopt.arg);
-		else if (c == 'q') min_query_len = parseInt(getopt.arg);
 	}
 	if (getopt.ind == args.length) {
 		print("Usage: paftools.js asmstat [options] <ref.fa.fai> <asm1.paf> [...]");
 		print("Options:");
-		print("  -q INT     ignore query shorter than INT [0]");
 		print("  -l INT     min alignment block length [" + min_seg_len + "]");
 		print("  -d FLOAT   max gap-compressed sequence divergence [" + max_diff + "]");
 		exit(1);
@@ -592,7 +587,7 @@ function paf_asmstat(args)
 	}
 	file.close();
 
-	function process_query(qblocks, qblock_len, bp, qi) {
+	function process_query(qblocks, qblock_len, bp) {
 		qblocks.sort(function(a,b) { return a[0]-b[0]; });
 		var last_k = null, last_blen = null, st = -1, en = -1, qcov = 0;
 		for (var k = 0; k < qblocks.length; ++k) {
@@ -617,7 +612,6 @@ function paf_asmstat(args)
 				var min = blen < last_blen? blen : last_blen;
 				var flank = k == 0? min : blen;
 				bp.push([flank, gap]);
-				qi.bp.push([flank, gap]);
 			}
 			last_k = k, last_blen = blen;
 		}
@@ -660,7 +654,7 @@ function paf_asmstat(args)
 		return (NM - n_gaps + n_gapo) / (n_M + n_gapo);
 	}
 
-	var labels = ['Length', 'l_cov', 'Rcov', 'Rdup', 'Qcov', 'NG75', 'NG50', 'NGA50', '#breaks', 'bp(' + min_seg_len + ',0)', 'bp(' + min_seg_len + ',10k)'];
+	var labels = ['Length', 'NG50', 'Coverage', 'Qcov', 'NGA50', '#breaks', 'bp(' + min_seg_len + ',0)', 'bp(' + min_seg_len + ',10k)'];
 	var rst = [];
 	for (var i = 0; i < labels.length; ++i)
 		rst[i] = [];
@@ -670,23 +664,17 @@ function paf_asmstat(args)
 	for (var i = 0; i < n_asm; ++i) {
 		var n_breaks = 0, qcov = 0;
 		var fn = args[getopt.ind + 1 + i];
-		var label = fn.replace(/.paf(.gz)?$/, "");
+		header.push(fn.replace(/.paf(.gz)?$/, ""));
-		header.push(label);
 		var ref_blocks = [], qblock_len = [], qblocks = [], bp = [];
-		var query = {}, qinfo = {};
+		var query = {};
 		var last_qname = null;
 		file = new File(fn);
 		while (file.readline(buf) >= 0) {
 			var m, line = buf.toString();
 			var t = line.split("\t");
 			t[1] = parseInt(t[1]);
-			if (t[1] < min_query_len) continue;
+			if (t.length >= 2) query[t[0]] = t[1];
-			if (t.length < 2) continue;
+			if (t.length < 9) continue;
-			query[t[0]] = t[1];
-			if (qinfo[t[0]] == null) qinfo[t[0]] = {};
-			qinfo[t[0]].len = t[1];
-			qinfo[t[0]].bp = [];
-			if (t.length < 9 || t[5] == "*") continue;
 			if (!/\ttp:A:[PI]/.test(line)) continue;
 			if ((m = /\tcg:Z:(\S+)/.exec(line)) == null) continue;
 			var cigar = m[1];
@@ -702,7 +690,7 @@ function paf_asmstat(args)
 			if (t[3] - t[2] < min_seg_len) continue;
 			if (t[0] != last_qname) {
 				if (last_qname != null)
-					qcov += process_query(qblocks, qblock_len, bp, qinfo[last_qname]);
+					qcov += process_query(qblocks, qblock_len, bp);
 				qblocks = [];
 				last_qname = t[0];
 			}
@@ -710,7 +698,7 @@ function paf_asmstat(args)
 			qblocks.push([t[2], t[3], t[4], t[5], t[7], t[8]]);
 		}
 		if (last_qname != null)
-			qcov += process_query(qblocks, qblock_len, bp, qinfo[last_qname]);
+			qcov += process_query(qblocks, qblock_len, bp);
 		file.close();
 
 		// compute NG50
@@ -720,8 +708,7 @@ function paf_asmstat(args)
 			asm_lens.push(query[ctg]);
 		}
 		rst[0][i] = asm_len;
-		rst[5][i] = N50(asm_lens, ref_len, 0.75);
+		rst[1][i] = N50(asm_lens, ref_len, 0.5);
-		rst[6][i] = N50(asm_lens, ref_len, 0.5);
 
 		// compute coverage
 		var l_cov = 0;
@@ -736,195 +723,22 @@ function paf_asmstat(args)
 			} else en = en > ref_blocks[j][2]? en : ref_blocks[j][2];
 		}
 		l_cov += en - st;
-		rst[1][i] = l_cov;
 		rst[2][i] = (100.0 * (l_cov / ref_len)).toFixed(2) + '%';
-		rst[4][i] = (100.0 * (qcov  / asm_len)).toFixed(2) + '%';
+		rst[3][i] = (100.0 * (qcov  / asm_len)).toFixed(2) + '%';
-
-		// compute cov1 and cov2+ lengths; see paf_call() for details
-		var c1_ctg = null, c1_start = 0, c1_end = 0, c1_len = 0;
-		for (var j = 0; j < ref_blocks.length; ++j) {
-			if (ref_blocks[j][0] != c1_ctg || ref_blocks[j][1] >= c1_end) {
-				if (c1_end > c1_start)
-					c1_len += c1_end - c1_start;
-				c1_ctg = ref_blocks[j][0], c1_start = ref_blocks[j][1], c1_end = ref_blocks[j][2];
-			} else if (ref_blocks[j][2] > c1_end) { // overlap
-				if (ref_blocks[j][1] > c1_start)
-					c1_len += ref_blocks[j][1] - c1_start;
-				c1_start = c1_end, c1_end = ref_blocks[j][2];
-			} else if (ref_blocks[j][2] > c1_start) { // contained
-				if (ref_blocks[j][1] > c1_start)
-					c1_len += ref_blocks[j][1] - c1_start;
-				c1_start = ref_blocks[j][2];
-			}
-			//print(ref_blocks[j][0], ref_blocks[j][1], ref_blocks[j][2], c1_start, c1_end, c1_len);
-		}
-		if (c1_end > c1_start)
-			c1_len += c1_end - c1_start;
-		rst[3][i] = (100 * (l_cov - c1_len) / l_cov).toFixed(2) + '%';
 
 		// compute NGA50
-		rst[7][i] = N50(qblock_len, ref_len, 0.5);
+		rst[4][i] = N50(qblock_len, ref_len, 0.5);
 
 		// compute break points
-		rst[8][i] = n_breaks;
+		rst[5][i] = n_breaks;
-		rst[9][i] = count_bp(bp, 500, 0);
+		rst[6][i] = count_bp(bp, 500, 0);
-		rst[10][i] = count_bp(bp, 500, 10000);
+		rst[7][i] = count_bp(bp, 500, 10000);
-
-		// nb-plot; NOT USED
-		/*
-		var qa = [];
-		for (var qn in qinfo)
-			qa.push([qinfo[qn].len, qinfo[qn].bp]);
-		qa = qa.sort(function(a, b) { return b[0] - a[0] });
-		var sum = 0, n_bp = 0, next_quantile = 0.1;
-		for (var j = 0; j < qa.length; ++j) {
-			sum += qa[j][0];
-			for (var k = 0; k < qa[j][1].length; ++k)
-				if (qa[j][1][k][0] >= bp_flank_len && qa[j][1][k][1] >= bp_gap_len)
-					++n_bp;
-			if (sum >= ref_len * next_quantile) {
-				print(label, Math.floor(next_quantile * 100 + .5), qa[j][0], (sum / n_bp).toFixed(0), n_bp);
-				next_quantile += 0.1;
-				if (next_quantile >= 1.0) break;
-			}
-		}
-		*/
-	}
-	buf.destroy();
-
-	if (bp_flank_len <= 0) {
-		print(header.join("\t"));
-		for (var i = 0; i < labels.length; ++i)
-			print(labels[i], rst[i].join("\t"));
-	}
-}
-
-function paf_asmgene(args)
-{
-	var c, opt = { min_cov:0.99, min_iden:0.99 }, print_err = false, auto_only = false;
-	while ((c = getopt(args, "i:c:ea")) != null)
-		if (c == 'i') opt.min_iden = parseFloat(getopt.arg);
-		else if (c == 'c') opt.min_cov = parseFloat(getopt.arg);
-		else if (c == 'e') print_err = true;
-		else if (c == 'a') auto_only = true;
-
-	var n_fn = args.length - getopt.ind;
-	if (n_fn < 2) {
-		print("Usage: paftools.js asmgene [options] <ref-splice.paf> <asm-splice.paf> [...]");
-		print("Options:");
-		print("  -i FLOAT     min identity [" + opt.min_iden + "]");
-		print("  -c FLOAT     min coverage [" + opt.min_cov + "]");
-		print("  -a           only evaluate genes mapped to the autosomes");
-		print("  -e           print fragmented/missing genes");
-		exit(1);
 	}
 
-	function process_query(opt, a) {
+	print(header.join("\t"));
-		var b = [], cnt = [0, 0, 0];
+	for (var i = 0; i < labels.length; ++i)
-		for (var j = 0; j < a.length; ++j) {
+		print(labels[i], rst[i].join("\t"));
-			if (a[j][4] < a[j][5] * opt.min_iden)
-				continue;
-			b.push(a[j].slice(0));
-		}
-		if (b.length == 0) return cnt;
-		// count full
-		var n_full = 0;
-		for (var j = 0; j < b.length; ++j)
-			if (b[j][3] - b[j][2] >= b[j][1] * opt.min_cov)
-				++n_full;
-		cnt[0] = n_full;
-		// compute coverage
-		b = b.sort(function(x, y) { return x[2] - y[2] });
-		var l_cov = 0, st = b[0][2], en = b[0][3];
-		for (var j = 1; j < b.length; ++j) {
-			if (b[j][2] <= en)
-				en = b[j][3] > en? b[j][3] : en;
-			else l_cov += en - st;
-		}
-		l_cov += en - st;
-		cnt[1] = l_cov / b[0][1];
-		cnt[2] = b.length;
-		return cnt;
-	}
 
-	var buf = new Bytes();
-	var gene = {}, header = [], refpos = {};
-	for (var i = getopt.ind; i < args.length; ++i) {
-		var fn = args[i];
-		var label = fn.replace(/.paf(.gz)?$/, "");
-		header.push(label);
-		var file = new File(fn), a = [];
-		while (file.readline(buf) >= 0) {
-			var t = buf.toString().split("\t");
-			var ql = parseInt(t[1]), qs = parseInt(t[2]), qe = parseInt(t[3]), mlen = parseInt(t[9]), blen = parseInt(t[10]), mapq = parseInt(t[11]);
-			if (i == getopt.ind) refpos[t[0]] = [t[0], t[1], t[5], t[7], t[8]];
-			if (gene[t[0]] == null) gene[t[0]] = [];
-			if (a.length && t[0] != a[0][0]) {
-				gene[a[0][0]][i - getopt.ind] = process_query(opt, a);
-				a = [];
-			}
-			a.push([t[0], ql, qs, qe, mlen, blen]);
-		}
-		if (a.length)
-			gene[t[0]][i - getopt.ind] = process_query(opt, a);
-		file.close();
-	}
-
-	// select the longest genes (not optimal, but should be good enough)
-	var gene_list = [], gene_nr = {};
-	for (var g in refpos)
-		gene_list.push(refpos[g]);
-	gene_list = gene_list.sort(function(a, b) { return a[2] < b[2]? -1 : a[2] > b[2]? 1 : a[3] - b[3] });
-	var last = 0;
-	for (var j = 1; j < gene_list.length; ++j) {
-		if (gene_list[j][2] != gene_list[last][2] || gene_list[j][3] >= gene_list[last][4]) {
-			gene_nr[gene_list[last][0]] = 1;
-			last = j;
-		} else if (gene_list[j][1] > gene_list[last][1]) {
-			last = j;
-		}
-	}
-	gene_nr[gene_list[last][0]] = 1;
-
-	// count and print
-	var col1 = ["full_sgl", "full_dup", "frag", "part50+", "part10+", "part10-"];
-	var rst = [];
-	for (var k = 0; k < col1.length; ++k) {
-		rst[k] = [];
-		for (var i = 0; i < n_fn; ++i)
-			rst[k][i] = 0;
-	}
-	for (var g in gene) {
-		if (gene[g][0] == null || gene[g][0][0] != 1) continue;
-		if (gene_nr[g] == null) continue;
-		if (auto_only && /^(chr)?[XY]$/.test(refpos[g][2])) continue;
-		for (var i = 0; i < n_fn; ++i) {
-			if (gene[g][i] == null) {
-				rst[4][i]++;
-				if (print_err) print('M', header[i], refpos[g].join("\t"));
-			} else if (gene[g][i][0] == 1) rst[0][i]++;
-			else if (gene[g][i][0] > 1) {
-				rst[1][i]++;
-				if (print_err) print('D', header[i], refpos[g].join("\t"));
-			} else if (gene[g][i][1] >= opt.min_cov) {
-				rst[2][i]++;
-				if (print_err) print('F', header[i], refpos[g].join("\t"));
-			} else if (gene[g][i][1] >= 0.5) {
-				rst[3][i]++;
-				if (print_err) print('5', header[i], refpos[g].join("\t"));
-			} else if (gene[g][i][1] >= 0.1) {
-				rst[4][i]++;
-				if (print_err) print('1', header[i], refpos[g].join("\t"));
-			} else {
-				rst[5][i]++;
-				if (print_err) print('0', header[i], refpos[g].join("\t")); // TODO: reduce code duplicates...
-			}
-		}
-	}
-	print('H', 'Metric', header.join("\t"));
-	for (var k = 0; k < rst.length; ++k) {
-		print('X', col1[k], rst[k].join("\t"));
-	}
 	buf.destroy();
 }
 
@@ -967,9 +781,7 @@ function paf_stat(args)
 			var t = line.split("\t", 12);
 			var m, rs, cigar = null, is_pri = false, is_sam = false, is_rev = false, tname = null;
 			var atlen = null, aqlen, qs, qe, mapq, ori_qlen;
-			if (t.length < 2) continue;
+			if (t[4] == '+' || t[4] == '-') { // PAF
-			if (t[4] == '+' || t[4] == '-' || t[4] == '*') { // PAF
-				if (t[4] == '*') continue; // unmapped
 				if (!/\ts2:i:\d+/.test(line)) {
 					++n_2nd;
 					continue;
@@ -1202,105 +1014,6 @@ function paf_bedcov(args)
 	warn("# target bases overlapping regions: " + hit_len + ' (' + (100.0 * hit_len / tot_len).toFixed(2) + '%)');
 }
 
-function paf_vcfpair(args)
-{
-	var c, is_male = false, sample = 'syndip', hgver = null;
-	var PAR = { '37':[[0, 2699520], [154931043, 155260560]] };
-	while ((c = getopt(args, "ms:g:")) != null) {
-		if (c == 'm') is_male = true;
-		else if (c == 's') sample = getopt.arg;
-		else if (c == 'g') hgver = getopt.arg;
-	}
-	if (is_male && (hgver == null || PAR[hgver] == null))
-		throw("for a male, -g must be specified to properly handle PARs on chrX");
-
-	if (getopt.ind == args.length) {
-		print("Usage: paftools.js vcfpair [options] <in.pair.vcf>");
-		print("Options:");
-		print("  -m       the sample is male");
-		print("  -g STR   human genome version '37' []");
-		print("  -s STR   sample name [" + sample + "]");
-		exit(1);
-	}
-
-	var re_ctg = is_male? /^(chr)?([0-9]+|X|Y)$/ : /^(chr)?([0-9]+|X)$/;
-	var label = ['1', '2'];
-	var buf = new Bytes();
-	var file = args[getopt.ind] == '-'? new File() : new File(args[getopt.ind]);
-	while (file.readline(buf) >= 0) {
-		var m, line = buf.toString();
-		if (line.charAt(0) == '#') {
-			if (/^##(source|reference)=/.test(line)) continue;
-			if ((m = /^##contig=.*ID=([^\s,]+)/.exec(line)) != null) {
-				if (!re_ctg.test(m[1])) continue;
-			} else if (/^#CHROM/.test(line)) {
-				var t = line.split("\t");
-				--t.length;
-				t[t.length-1] = sample;
-				line = t.join("\t");
-				print('##FILTER=<ID=HET1,Description="Heterozygous in the first haplotype">');
-				print('##FILTER=<ID=HET2,Description="Heterozygous in the second haplotype">');
-				print('##FILTER=<ID=GAP1,Description="Uncalled in the first haplotype">');
-				print('##FILTER=<ID=GAP2,Description="Uncalled in the second haplotype">');
-			}
-			print(line);
-			continue;
-		}
-		var t = line.split("\t");
-		if (!re_ctg.test(t[0])) continue;
-		var GT = null, AD = null, FILTER = [], HT = [null, null];
-		for (var i = 0; i < 2; ++i) {
-			if ((m = /^(\.|[0-9]+)\/(\.|[0-9]+):(\S+)/.exec(t[9+i])) == null) {
-				warn(line);
-				throw Error("malformatted VCF");
-			}
-			var s = m[3].split(",");
-			if (AD == null) {
-				AD = [];
-				for (var j = 0; j < s.length; ++j)
-					AD[j] = 0;
-			}
-			for (var j = 0; j < s.length; ++j)
-				AD[j] += parseInt(s[j]);
-			if (m[1] == '.') {
-				FILTER.push('GAP' + label[i]);
-				HT[i] = '.';
-			} else if (m[1] != m[2]) {
-				FILTER.push('HET' + label[i]);
-				HT[i] = '.';
-			} else HT[i] = m[1];
-		}
-		--t.length;
-		// test if this is in a haploid region
-		var hap = 0, st = parseInt(t[1]), en = st + t[3].length;
-		if (is_male) {
-			if (/^(chr)?X/.test(t[0])) {
-				if (hgver != null && PAR[hgver] != null) {
-					var r = PAR[hgver], in_par = false;
-					for (var i = 0; i < r.length; ++i)
-						if (r[i][0] <= st && en <= r[i][1])
-							in_par = true;
-					hap = in_par? 0 : 2;
-				}
-			} else if (/^(chr)?Y/.test(t[0])) {
-				hap = 1;
-			}
-		}
-		// special treatment for haploid regions
-		if (hap > 0 && FILTER.length == 1) {
-			if ((hap == 2 && FILTER[0] == "GAP1") || (hap == 1 && FILTER[0] == "GAP2"))
-				FILTER.length = 0;
-		}
-		// update VCF
-		t[5] = 30; // fake QUAL
-		t[6] = FILTER.length? FILTER.join(";") : ".";
-		t[9] = HT.join("|") + ":" + AD.join(",");
-		print(t.join("\t"));
-	}
-	file.close();
-	buf.destroy();
-}
-
 /**********************
  * Conversion related *
  **********************/
@@ -1469,21 +1182,15 @@ function paf_view(args)
 
 function paf_gff2bed(args)
 {
-	var c, fn_ucsc_fai = null, is_short = false, keep_gff = false, print_junc = false;
+	var c, fn_ucsc_fai = null, is_short = false, keep_gff = false;
-	while ((c = getopt(args, "u:sgj")) != null) {
+	while ((c = getopt(args, "u:sg")) != null) {
 		if (c == 'u') fn_ucsc_fai = getopt.arg;
 		else if (c == 's') is_short = true;
 		else if (c == 'g') keep_gff = true;
-		else if (c == 'j') print_junc = true;
 	}
 
 	if (getopt.ind == args.length) {
-		print("Usage: paftools.js gff2bed [options] <in.gff>");
+		print("Usage: paftools.js gff2bed [-g] [-u ucsc-genome.fa.fai] <in.gff>");
-		print("Options:");
-		print("  -j       Output junction BED");
-		print("  -s       Print names in the short form");
-		print("  -u FILE  hg38.fa.fai for chr name conversion");
-		print("  -g       Output GFF (used with -u)");
 		exit(1);
 	}
 
@@ -1515,16 +1222,11 @@ function paf_gff2bed(args)
 		'misc_RNA':'0,192,0'
 	};
 
-	function print_bed12(exons, cds_st, cds_en, is_short, print_junc)
+	function print_bed12(exons, cds_st, cds_en, is_short)
 	{
 		if (exons.length == 0) return;
 		var name = is_short? exons[0][7] + "|" + exons[0][5] : exons[0].slice(4, 7).join("|");
 		var a = exons.sort(function(a,b) {return a[1]-b[1]});
-		if (print_junc) {
-		for (var i = 1; i < a.length; ++i)
-			print(a[i][0], a[i-1][2], a[i][1], name, 1000, a[i][3]);
-			return;
-		}
 		var sizes = [], starts = [], st, en;
 		st = a[0][1];
 		en = a[a.length - 1][2];
@@ -1577,7 +1279,7 @@ function paf_gff2bed(args)
 		if (type == "" && biotype != "") type = biotype;
 		if (id == null) throw Error("No transcript_id");
 		if (id != last_id) {
-			print_bed12(exons, cds_st, cds_en, is_short, print_junc);
+			print_bed12(exons, cds_st, cds_en, is_short);
 			exons = [], cds_st = 1<<30, cds_en = 0;
 			last_id = id;
 		}
@@ -1595,7 +1297,7 @@ function paf_gff2bed(args)
 		}
 	}
 	if (last_id != null)
-		print_bed12(exons, cds_st, cds_en, is_short, print_junc);
+		print_bed12(exons, cds_st, cds_en, is_short);
 
 	file.close();
 	buf.destroy();
@@ -1603,16 +1305,11 @@ function paf_gff2bed(args)
 
 function paf_sam2paf(args)
 {
-	var c, pri_only = false, long_cs = false;
+	var c, pri_only = false, use_eq = false;
-	while ((c = getopt(args, "pL")) != null) {
+	while ((c = getopt(args, "p")) != null)
 		if (c == 'p') pri_only = true;
-		else if (c == 'L') long_cs = true;
-	}
 	if (args.length == getopt.ind) {
-		print("Usage: paftools.js sam2paf [options] <in.sam>");
+		print("Usage: paftools.js sam2paf [-p] <in.sam>");
-		print("Options:");
-		print("  -p      convert primary or supplementary alignments only");
-		print("  -L      output the cs tag in the long form");
 		exit(1);
 	}
 
@@ -1641,14 +1338,13 @@ function paf_sam2paf(args)
 		var tlen = ctg_len[t[2]];
 		if (tlen == null) throw Error("at line " + lineno + ": can't find the length of contig " + t[2]);
 		// find tags
-		var nn = 0, NM = null, MD = null, cs_str = null, md_list = [];
+		var nn = 0, NM = null, MD = null, md_list = [];
 		while ((m = re_tag.exec(line)) != null) {
 			if (m[1] == "NM:i") NM = parseInt(m[2]);
 			else if (m[1] == "nn:i") nn = parseInt(m[2]);
 			else if (m[1] == "MD:Z") MD = m[2];
-			else if (m[1] == "cs:Z") cs_str = m[2];
 		}
-		if (t[9] == '*') MD = cs_str = null;
+		if (t[9] == '*') MD = null;
 		// infer various lengths from CIGAR
 		var clip = [0, 0], soft_clip = 0, I = [0, 0], D = [0, 0], M = 0, N = 0, mm = 0, have_M = false, have_ext = false, cigar = [];
 		while ((m = re.exec(t[5])) != null) {
@@ -1684,8 +1380,8 @@ function paf_sam2paf(args)
 		}
 		// parse MD
 		var cs = [];
-		if (MD != null && cs_str == null && t[9] != "*") {
+		if (MD != null) {
-			var k = 0, cx = 0, cy = 0, mx = 0, my = 0; // cx: cigar ref position; cy: cigar query; mx: MD ref; my: MD query
+			var k = 0, cx = 0, cy = 0, mx = 0, my = 0;
 			while ((m = re_MD.exec(MD)) != null) {
 				if (m[2] != null) { // deletion from the reference
 					var len = m[2].length - 1;
@@ -1699,15 +1395,13 @@ function paf_sam2paf(args)
 							if (my + ml < cy + cl) {
 								if (ml > 0) {
 									if (m[3] != null) cs.push('*', m[3], t[9][my]);
-									else if (long_cs) cs.push('=', t[9].substr(my, ml));
 									else cs.push(':', ml);
 								}
 								mx += ml, my += ml, ml = 0;
 								break;
 							} else {
 								var dl = cy + cl - my;
-								if (long_cs) cs.push('=', t[9].substr(my, dl));
+								cs.push(':', dl);
-								else cs.push(':', dl);
 								cx += cl, cy += cl, ++k;
 								mx += dl, my += dl, ml -= dl;
 							}
@@ -1752,8 +1446,7 @@ function paf_sam2paf(args)
 		var tags = ["tp:A:" + type];
 		if (NM != null) tags.push("mm:i:"+mm);
 		tags.push("gn:i:"+(I[1]+D[1]), "go:i:"+(I[0]+D[0]), "cg:Z:" + t[5].replace(/\d+[SH]/g, ''));
-		if (cs_str != null) tags.push("cs:Z:" + cs_str);
+		if (cs.length > 0) tags.push("cs:Z:" + cs.join(""));
-		else if (cs.length > 0) tags.push("cs:Z:" + cs.join(""));
 		// print out
 		var a = [qname, qlen, qs, qe, flag&16? '-' : '+', t[2], tlen, ts, te, mlen, blen, t[4]];
 		print(a.join("\t"), tags.join("\t"));
@@ -2496,7 +2189,6 @@ function main(args)
 		print("");
 		print("  stat       collect basic mapping information in PAF/SAM");
 		print("  asmstat    collect basic assembly information");
-		print("  asmgene    evaluate gene completeness (EXPERIMENTAL)");
 		print("  liftover   simplistic liftOver");
 		print("  call       call variants from asm-to-ref alignment with the cs tag");
 		print("  bedcov     compute the number of bases covered");
@@ -2518,9 +2210,7 @@ function main(args)
 	else if (cmd == 'gff2bed') paf_gff2bed(args);
 	else if (cmd == 'stat') paf_stat(args);
 	else if (cmd == 'asmstat') paf_asmstat(args);
-	else if (cmd == 'asmgene') paf_asmgene(args);
 	else if (cmd == 'liftover' || cmd == 'liftOver') paf_liftover(args);
-	else if (cmd == 'vcfpair') paf_vcfpair(args);
 	else if (cmd == 'call') paf_call(args);
 	else if (cmd == 'mapeval') paf_mapeval(args);
 	else if (cmd == 'bedcov') paf_bedcov(args);

mmpriv.h

@@ -31,6 +31,13 @@
 #define MALLOC(type, len) ((type*)malloc((len) * sizeof(type)))
 #define CALLOC(type, len) ((type*)calloc((len), sizeof(type)))
 
+#define REALLOC(ptr, len) ((ptr) = (__typeof__(ptr))realloc((ptr), (len) * sizeof(*(ptr))))
+
+#define EXPAND(a, m) do { \
+		(m) = (m)? (m) + ((m)>>1) : 16; \
+		REALLOC((a), (m)); \
+	} while (0)
+
 #ifdef __cplusplus
 extern "C" {
 #endif
@@ -61,15 +68,14 @@ void mm_sketch(void *km, const char *str, int len, int w, int k, uint32_t rid, i
 
 void mm_write_sam_hdr(const mm_idx_t *mi, const char *rg, const char *ver, int argc, char *argv[]);
 void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag);
-void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, void *km, int opt_flag, int rep_len);
 void mm_write_sam(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int n_regs, const mm_reg1_t *regs);
 void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regs, const mm_reg1_t *const* regs, void *km, int opt_flag);
-void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int seg_idx, int reg_idx, int n_seg, const int *n_regss, const mm_reg1_t *const* regss, void *km, int opt_flag, int rep_len);
 
 void mm_idxopt_init(mm_idxopt_t *opt);
 const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n);
 int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f);
-mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int max_iter, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
+int mm_idx_getseq2(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq, int8_t *b);
+mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km);
 mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a);
 
 mm_reg1_t *mm_gen_regs(void *km, uint32_t hash, int qlen, int n_u, uint64_t *u, mm128_t *a);

options.c

@@ -23,7 +23,6 @@ void mm_mapopt_init(mm_mapopt_t *opt)
 	opt->max_gap = 5000;
 	opt->max_gap_ref = -1;
 	opt->max_chain_skip = 25;
-	opt->max_chain_iter = 5000;
 
 	opt->mask_level = 0.5f;
 	opt->pri_ratio = 0.8f;
@@ -120,16 +119,13 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
 		mo->mid_occ = 1000;
 		mo->max_occ = 5000;
 		mo->mini_batch_size = 50000000;
-	} else if (strncmp(preset, "splice", 6) == 0 || strcmp(preset, "cdna") == 0) {
+	} else if (strcmp(preset, "splice") == 0 || strcmp(preset, "cdna") == 0) {
 		io->flag = 0, io->k = 15, io->w = 5;
 		mo->flag |= MM_F_SPLICE | MM_F_SPLICE_FOR | MM_F_SPLICE_REV | MM_F_SPLICE_FLANK;
 		mo->max_gap = 2000, mo->max_gap_ref = mo->bw = 200000;
 		mo->a = 1, mo->b = 2, mo->q = 2, mo->e = 1, mo->q2 = 32, mo->e2 = 0;
 		mo->noncan = 9;
-		mo->junc_bonus = 9;
 		mo->zdrop = 200, mo->zdrop_inv = 100; // because mo->a is halved
-		if (strcmp(preset, "splice:hq") == 0)
-			mo->junc_bonus = 5, mo->b = 4, mo->q = 6, mo->q2 = 24;
 	} else return -1;
 	return 0;
 }
@@ -163,11 +159,6 @@ int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
 			fprintf(stderr, "[ERROR]\033[1;31m --for-only and --rev-only can't be applied at the same time\033[0m\n");
 		return -3;
 	}
-	if (mo->e <= 0 || mo->q <= 0) {
-		if (mm_verbose >= 1)
-			fprintf(stderr, "[ERROR]\033[1;31m -O and -E must be positive\033[0m\n");
-		return -1;
-	}
 	if ((mo->q != mo->q2 || mo->e != mo->e2) && !(mo->e > mo->e2 && mo->q + mo->e < mo->q2 + mo->e2)) {
 		if (mm_verbose >= 1)
 			fprintf(stderr, "[ERROR]\033[1;31m dual gap penalties violating E1>E2 and O1+E1<O2+E2\033[0m\n");

pe.c

@@ -54,7 +54,7 @@ void mm_set_pe_thru(const int *qlens, int *n_regs, mm_reg1_t **regs)
 	if (n_pri[0] == 1 && n_pri[1] == 1) {
 		mm_reg1_t *p = &regs[0][pri[0]];
 		mm_reg1_t *q = &regs[1][pri[1]];
-		if (p->rid == q->rid && p->rev == q->rev && abs(p->rs - q->rs) < 3 && abs(p->re - q->re) < 3
+		if (p->rid == q->rid && p->rev == q->rev && abs(p->rs - q->rs) < 3 && abs(p->re - p->re) < 3
 			&& ((p->qs == 0 && qlens[1] - q->qe == 0) || (q->qs == 0 && qlens[0] - p->qe == 0)))
 		{
 			p->pe_thru = q->pe_thru = 1;

python/README.rst

@@ -114,12 +114,6 @@ This method retrieves a (sub)sequence from the index and returns it as a Python
 string. :code:`None` is returned if :code:`name` is not present in the index or
 the start/end coordinates are invalid.
 
-.. code:: python
-
-	mappy.Aligner.seq_names
-
-This property gives the array of sequence names in the index.
-
 Class mappy.Alignment
 ~~~~~~~~~~~~~~~~~~~~~
 

python/cmappy.h

@@ -20,6 +20,11 @@ typedef struct {
 	uint32_t *cigar32;
 } mm_hitpy_t;
 
+typedef struct {
+	int32_t n, m;
+	mm_idx_bed_t *r;
+} mm_bedpy_t;
+
 static inline void mm_reg2hitpy(const mm_idx_t *mi, mm_reg1_t *r, mm_hitpy_t *h)
 {
 	h->ctg = mi->seq[r->rid].name;
@@ -149,4 +154,33 @@ static mm_idx_t *mappy_idx_seq(int w, int k, int is_hpc, int bucket_bits, const
 	return mi;
 }
 
+static mm_bedpy_t *mappy_bed_new(void)
+{
+	return (mm_bedpy_t*)calloc(1, sizeof(mm_bedpy_t));
+}
+
+static int mappy_bed_add(mm_bedpy_t *bed, mm_idx_t *mi, const char *name, uint32_t st, uint32_t en)
+{
+	mm_idx_bed_t *b;
+	int id;
+	if (mi->h == 0) mm_idx_index_name(mi);
+	if (bed->n == bed->m) {
+		bed->m = bed->m? bed->m + (bed->m>>1) : 16;
+		bed->r = (mm_idx_bed_t*)realloc(bed->r, sizeof(mm_idx_bed_t) * bed->m);
+	}
+	id = mm_idx_name2id(mi, name);
+	if (id < 0 || st >= en) return -1;
+	if (en > mi->seq[id].len) en = mi->seq[id].len;
+	b = &bed->r[bed->n++];
+	b->x = (uint64_t)id << 32 | st;
+	b->end = en, b->idx = -1;
+	return 0;
+}
+
+static void mappy_bed_finalize(mm_bedpy_t *bed, mm_idx_t *mi)
+{
+	mm_idx_bed_attach(mi, bed->n, bed->r); // bed->r is now owned by mi and will be deallocated with it
+	free(bed);
+}
+
 #endif

python/cmappy.pxd

@@ -10,14 +10,13 @@ cdef extern from "minimap.h":
 		uint64_t batch_size
 
 	ctypedef struct mm_mapopt_t:
-		int64_t flag
 		int seed
 		int sdust_thres
-		int max_qlen
+		int flag
 		int bw
 		int max_gap, max_gap_ref
 		int max_frag_len
-		int max_chain_skip, max_chain_iter
+		int max_chain_skip
 		int min_cnt
 		int min_chain_score
 		float mask_level
@@ -25,11 +24,10 @@ cdef extern from "minimap.h":
 		int best_n
 		int max_join_long, max_join_short
 		int min_join_flank_sc
-		float min_join_flank_ratio
+		float min_join_flank_ratio;
 		int a, b, q, e, q2, e2
 		int sc_ambi
 		int noncan
-		int junc_bonus
 		int zdrop, zdrop_inv
 		int end_bonus
 		int min_dp_max
@@ -42,7 +40,6 @@ cdef extern from "minimap.h":
 		int32_t mid_occ
 		int32_t max_occ
 		int mini_batch_size
-		int64_t max_sw_mat
 		const char *split_prefix
 
 	int mm_set_opt(char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
@@ -56,6 +53,10 @@ cdef extern from "minimap.h":
 		uint64_t offset
 		uint32_t len
 
+	ctypedef struct mm_idx_bed_t:
+		uint64_t x
+		int32_t end, idx
+
 	ctypedef struct mm_idx_bucket_t:
 		pass
 
@@ -65,6 +66,7 @@ cdef extern from "minimap.h":
 		mm_idx_seq_t *seq
 		uint32_t *S
 		mm_idx_bucket_t *B
+		mm_idx_bed_t *R
 		void *km
 		void *h
 
@@ -115,6 +117,14 @@ cdef extern from "cmappy.h":
 	char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int en, int *l)
 	mm_idx_t *mappy_idx_seq(int w, int k, int is_hpc, int bucket_bits, const char *seq, int l)
 
+	ctypedef struct mm_bedpy_t:
+		int32_t n, m
+		mm_idx_bed_t *r
+
+	mm_bedpy_t *mappy_bed_new()
+	int mappy_bed_add(mm_bedpy_t *bed, mm_idx_t *mi, const char *name, uint32_t st, uint32_t en)
+	void mappy_bed_finalize(mm_bedpy_t *bed, mm_idx_t *mi)
+
 	ctypedef struct kstring_t:
 		unsigned l, m
 		char *s

python/mappy.pyx

@@ -3,7 +3,7 @@ from libc.stdlib cimport free
 cimport cmappy
 import sys
 
-__version__ = '2.17'
+__version__ = '2.13'
 
 cmappy.mm_reset_timer()
 
@@ -112,7 +112,7 @@ cdef class Aligner:
 	cdef cmappy.mm_idxopt_t idx_opt
 	cdef cmappy.mm_mapopt_t map_opt
 
-	def __cinit__(self, fn_idx_in=None, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None, extra_flags=None, seq=None, scoring=None):
+	def __cinit__(self, fn_idx_in=None, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None, extra_flags=None, seq=None, scoring=None, bed=None):
 		cmappy.mm_set_opt(NULL, &self.idx_opt, &self.map_opt) # set the default options
 		if preset is not None:
 			cmappy.mm_set_opt(str.encode(preset), &self.idx_opt, &self.map_opt) # apply preset
@@ -137,12 +137,13 @@ cdef class Aligner:
 					self.map_opt.sc_ambi = scoring[6]
 
 		cdef cmappy.mm_idx_reader_t *r;
+		cdef cmappy.mm_bedpy_t *bed_agg;
 
 		if seq is None:
 			if fn_idx_out is None:
 				r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, NULL)
 			else:
-				r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, str.encode(fn_idx_out))
+				r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, fn_idx_out)
 			if r is not NULL:
 				self._idx = cmappy.mm_idx_reader_read(r, n_threads) # NB: ONLY read the first part
 				cmappy.mm_idx_reader_close(r)
@@ -153,6 +154,14 @@ cdef class Aligner:
 			cmappy.mm_mapopt_update(&self.map_opt, self._idx)
 			self.map_opt.mid_occ = 1000 # don't filter high-occ seeds
 
+		if bed is not None:
+			bed_agg = cmappy.mappy_bed_new()
+			for b in bed:
+				if len(b) < 3: en = int(b[1]) + 1
+				else: en = int(b[2])
+				cmappy.mappy_bed_add(bed_agg, self._idx, str.encode(b[0]), int(b[1]), en)
+			cmappy.mappy_bed_finalize(bed_agg, self._idx)
+
 	def __dealloc__(self):
 		if self._idx is not NULL:
 			cmappy.mm_idx_destroy(self._idx)
@@ -221,16 +230,6 @@ cdef class Aligner:
 	@property
 	def n_seq(self): return self._idx.n_seq
 
-	@property
-	def seq_names(self):
-		cdef char *p
-		sn = []
-		for i in range(self._idx.n_seq):
-			p = self._idx.seq[i].name
-			s = p if isinstance(p, str) else p.decode()
-			sn.append(s)
-		return sn
-
 def fastx_read(fn, read_comment=False):
 	cdef cmappy.kseq_t *ks
 	ks = cmappy.mm_fastx_open(str.encode(fn))

setup.py

@@ -33,7 +33,7 @@ def readme():
 
 setup(
 	name = 'mappy',
-	version = '2.17',
+	version = '2.13',
 	url = 'https://github.com/lh3/minimap2',
 	description = 'Minimap2 python binding',
 	long_description = readme(),

splitidx.c

@@ -11,11 +11,8 @@ FILE *mm_split_init(const char *prefix, const mm_idx_t *mi)
 	uint32_t i, k = mi->k;
 	fn = (char*)calloc(strlen(prefix) + 10, 1);
 	sprintf(fn, "%s.%.4d.tmp", prefix, mi->index);
-	if ((fp = fopen(fn, "wb")) == NULL) {
+	fp = fopen(fn, "wb");
-		if (mm_verbose >= 1)
+	assert(fp);
-			fprintf(stderr, "[ERROR]\033[1;31m failed to write to temporary file '%s'\033[0m\n", fn);
-		exit(1);
-	}
 	mm_err_fwrite(&k, 4, 1, fp);
 	mm_err_fwrite(&mi->n_seq, 4, 1, fp);
 	for (i = 0; i < mi->n_seq; ++i) {