diff --git a/NEWS.md b/NEWS.md
index e5060bc..7cb5401 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,88 @@
+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)
-----------------------------------
diff --git a/README.md b/README.md
index d3ba33f..11088d9 100644
--- a/README.md
+++ b/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 -d MT-human.mmi test/MT-human.fa
-./minimap2 -a MT-human.mmi test/MT-orang.fa > test.sam
+./minimap2 -x map-ont -d MT-human-ont.mmi test/MT-human.fa
+./minimap2 -a MT-human-ont.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
@@ -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.14/minimap2-2.14_x64-linux.tar.bz2 | tar -jxvf -
-./minimap2-2.14_x64-linux/minimap2
+curl -L https://github.com/lh3/minimap2/releases/download/v2.16/minimap2-2.16_x64-linux.tar.bz2 | tar -jxvf -
+./minimap2-2.16_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
@@ -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. [doi:10.1093/bioinformatics/bty191][doi]
+> *Bioinformatics*, **34**:3094-3100. [doi:10.1093/bioinformatics/bty191][doi]
## Developers' Guide
diff --git a/align.c b/align.c
index 21808b9..c3b04d3 100644
--- a/align.c
+++ b/align.c
@@ -613,6 +613,7 @@ 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;
}
}
@@ -626,6 +627,7 @@ 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;
@@ -665,7 +667,7 @@ 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);
- if (qs > 0 && rs > 0) { // left extension
+ if (qs > 0 && rs > 0) { // left extension; probably the condition can be changed to "qs > qs0 && rs > rs0"
qseq = &qseq0[rev][qs0];
mm_idx_getseq(mi, rid, rs0, rs, tseq);
mm_seq_rev(qs - qs0, qseq);
diff --git a/chain.c b/chain.c
index 4947a3e..be9d51b 100644
--- a/chain.c
+++ b/chain.c
@@ -19,7 +19,7 @@ static inline int ilog2_32(uint32_t v)
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 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 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)
{ // 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,6 +28,10 @@ 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);
@@ -45,6 +49,7 @@ 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;
diff --git a/cookbook.md b/cookbook.md
index d43d190..ce581bc 100644
--- a/cookbook.md
+++ b/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.14/minimap2-2.14_x64-linux.tar.bz2 | tar jxf -
-cp minimap2-2.14_x64-linux/{minimap2,k8,paftools.js} . # copy executables
+curl -L https://github.com/lh3/minimap2/releases/download/v2.16/minimap2-2.16_x64-linux.tar.bz2 | tar jxf -
+cp minimap2-2.16_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 -
diff --git a/format.c b/format.c
index e81ef32..29138d1 100644
--- a/format.c
+++ b/format.c
@@ -270,7 +270,7 @@ double mm_event_identity(const mm_reg1_t *r)
if (op == 1 || op == 2)
++n_gapo, n_gap += len;
}
- return (double)r->mlen / (r->blen - r->p->n_ambi - n_gap + n_gapo);
+ return (double)r->mlen / (r->blen - n_gap + n_gapo);
}
static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
@@ -301,11 +301,12 @@ static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
if (r->split) mm_sprintf_lite(s, "\tzd:i:%d", r->split);
}
-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)
{
s->l = 0;
if (r == 0) {
- mm_sprintf_lite(s, "%s\t%d", t->name, t->l_seq);
+ mm_sprintf_lite(s, "%s\t%d\t0\t0\t*\t*\t0\t0\t0\t0\t0\t0", 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]);
@@ -315,6 +316,7 @@ void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m
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:");
@@ -327,6 +329,11 @@ void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m
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];
@@ -368,6 +375,7 @@ 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]);
@@ -376,7 +384,7 @@ static void write_sam_cigar(kstring_t *s, int sam_flag, int in_tag, int qlen, co
}
}
-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)
+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)
{
const int max_bam_cigar_op = 65535;
int flag, n_regs = n_regss[seg_idx], cigar_in_tag = 0;
@@ -527,6 +535,7 @@ void mm_write_sam2(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);
@@ -534,6 +543,11 @@ void mm_write_sam2(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;
diff --git a/hit.c b/hit.c
index 92aa6b7..f43b0d6 100644
--- a/hit.c
+++ b/hit.c
@@ -449,6 +449,7 @@ 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;
diff --git a/ketopt.h b/ketopt.h
index 70193a5..8ae1811 100644
--- a/ketopt.h
+++ b/ketopt.h
@@ -73,13 +73,17 @@ 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_matches = 0;
- const ko_longopt_t *o = 0;
+ int k, n_exact = 0, n_partial = 0;
+ const ko_longopt_t *o = 0, *o_exact = 0, *o_partial = 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)
- ++n_matches, o = &longopts[k];
- if (n_matches == 1) {
+ 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];
+ else ++n_partial, o_partial = &longopts[k];
+ }
+ 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') {
diff --git a/main.c b/main.c
index 12face7..aa0cf91 100644
--- a/main.c
+++ b/main.c
@@ -6,7 +6,7 @@
#include "mmpriv.h"
#include "ketopt.h"
-#define MM_VERSION "2.14-r903-avx2-dirty"
+#define MM_VERSION "2.16-r922"
#ifdef __linux__
#include
@@ -62,6 +62,7 @@ static ko_longopt_t long_options[] = {
{ "hard-mask-level",ko_no_argument, 336 },
{ "cap-sw-mem", ko_required_argument, 337 },
{ "max-qlen", ko_required_argument, 338 },
+ { "max-chain-iter", ko_required_argument, 339 },
{ "help", ko_no_argument, 'h' },
{ "max-intron-len", ko_required_argument, 'G' },
{ "version", ko_no_argument, 'V' },
@@ -99,7 +100,7 @@ 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:yYP";
+ 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:";
ketopt_t o = KETOPT_INIT;
mm_mapopt_t opt;
mm_idxopt_t ipt;
@@ -165,12 +166,21 @@ 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
@@ -268,7 +278,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 minizer window size [%d]\n", ipt.w);
+ fprintf(fp_help, " -w INT minimizer 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");
@@ -293,7 +303,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, " -Q don't output base quality in SAM\n");
+ fprintf(fp_help, " -o FILE output alignments to FILE [stdout]\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");
diff --git a/map.c b/map.c
index 16fbe06..41de91e 100644
--- a/map.c
+++ b/map.c
@@ -312,7 +312,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->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->max_chain_iter, 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 +334,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->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->max_chain_iter, 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 +584,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_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);
+ 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]);
else
- mm_write_paf(&p->str, mi, t, r, km, p->opt->flag);
+ mm_write_paf3(&p->str, mi, t, r, km, p->opt->flag, s->rep_len[i]);
mm_err_puts(p->str.s);
}
} 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_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);
+ 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]);
else
- mm_write_paf(&p->str, mi, t, 0, 0, p->opt->flag);
+ mm_write_paf3(&p->str, mi, t, 0, 0, p->opt->flag, s->rep_len[i]);
mm_err_puts(p->str.s);
}
}
diff --git a/minimap.h b/minimap.h
index 0b44325..176342b 100644
--- a/minimap.h
+++ b/minimap.h
@@ -112,7 +112,7 @@ typedef struct {
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;
+ int max_chain_skip, max_chain_iter;
int min_cnt; // min number of minimizers on each chain
int min_chain_score; // min chaining score
diff --git a/minimap2.1 b/minimap2.1
index b9badbb..be1f9ed 100644
--- a/minimap2.1
+++ b/minimap2.1
@@ -1,4 +1,4 @@
-.TH minimap2 1 "12 December 2018" "minimap2-2.14-dirty (r894)" "Bioinformatics tools"
+.TH minimap2 1 "28 Feburary 2019" "minimap2-2.16 (r922)" "Bioinformatics tools"
.SH NAME
.PP
minimap2 - mapping and alignment between collections of DNA sequences
@@ -232,13 +232,19 @@ Honor option
and disable a heurstic to save unmapped subsequences.
.TP
.BI --max-chain-skip \ INT
-A heuristics that stops chaining early [50]. Minimap2 uses dynamic programming
+A heuristics that stops chaining early [25]. 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
@@ -384,6 +390,11 @@ Set 0 to disable [0].
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
@@ -463,7 +474,9 @@ Filter out query sequences longer than
.IR NUM .
.TP
.B --paf-no-hit
-In PAF, output query name and length for an unmapped sequence.
+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 --version
Print version number to stdout
@@ -612,6 +625,7 @@ 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
diff --git a/misc/paftools.js b/misc/paftools.js
index c84ffa8..459eebd 100755
--- a/misc/paftools.js
+++ b/misc/paftools.js
@@ -1,6 +1,6 @@
#!/usr/bin/env k8
-var paftools_version = '2.14-r895-dirty';
+var paftools_version = '2.16-r922';
/*****************************
***** Library functions *****
@@ -433,6 +433,7 @@ 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;
@@ -680,13 +681,12 @@ function paf_asmstat(args)
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 (qinfo[t[0]] == null) qinfo[t[0]] = {};
- qinfo[t[0]].len = t[1];
- qinfo[t[0]].bp = [];
- }
- if (t.length < 9) continue;
+ if (t.length < 2) 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];
@@ -967,7 +967,9 @@ 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[4] == '+' || t[4] == '-') { // PAF
+ if (t.length < 2) continue;
+ if (t[4] == '+' || t[4] == '-' || t[4] == '*') { // PAF
+ if (t[4] == '*') continue; // unmapped
if (!/\ts2:i:\d+/.test(line)) {
++n_2nd;
continue;
@@ -1590,11 +1592,16 @@ function paf_gff2bed(args)
function paf_sam2paf(args)
{
- var c, pri_only = false, use_eq = false;
- while ((c = getopt(args, "p")) != null)
+ var c, pri_only = false, long_cs = false;
+ while ((c = getopt(args, "pL")) != null) {
if (c == 'p') pri_only = true;
+ else if (c == 'L') long_cs = true;
+ }
if (args.length == getopt.ind) {
- print("Usage: paftools.js sam2paf [-p] ");
+ print("Usage: paftools.js sam2paf [options] ");
+ print("Options:");
+ print(" -p convert primary or supplementary alignments only");
+ print(" -L output the cs tag in the long form");
exit(1);
}
@@ -1623,13 +1630,14 @@ 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, md_list = [];
+ var nn = 0, NM = null, MD = null, cs_str = 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 = null;
+ if (t[9] == '*') MD = cs_str = 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) {
@@ -1665,8 +1673,8 @@ function paf_sam2paf(args)
}
// parse MD
var cs = [];
- if (MD != null) {
- var k = 0, cx = 0, cy = 0, mx = 0, my = 0;
+ if (MD != null && cs_str == null && t[9] != "*") {
+ var k = 0, cx = 0, cy = 0, mx = 0, my = 0; // cx: cigar ref position; cy: cigar query; mx: MD ref; my: MD query
while ((m = re_MD.exec(MD)) != null) {
if (m[2] != null) { // deletion from the reference
var len = m[2].length - 1;
@@ -1680,13 +1688,15 @@ 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;
- cs.push(':', dl);
+ if (long_cs) cs.push('=', t[9].substr(my, dl));
+ else cs.push(':', dl);
cx += cl, cy += cl, ++k;
mx += dl, my += dl, ml -= dl;
}
@@ -1731,7 +1741,8 @@ 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.length > 0) tags.push("cs:Z:" + cs.join(""));
+ if (cs_str != null) tags.push("cs:Z:" + cs_str);
+ 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"));
diff --git a/mmpriv.h b/mmpriv.h
index 2688bca..c9b91cd 100644
--- a/mmpriv.h
+++ b/mmpriv.h
@@ -61,13 +61,15 @@ 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 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 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);
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);
diff --git a/options.c b/options.c
index b6d9004..eac73e3 100644
--- a/options.c
+++ b/options.c
@@ -23,6 +23,7 @@ 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;
diff --git a/python/README.rst b/python/README.rst
index 482c735..3d0fae8 100644
--- a/python/README.rst
+++ b/python/README.rst
@@ -114,6 +114,12 @@ 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
~~~~~~~~~~~~~~~~~~~~~
diff --git a/python/cmappy.pxd b/python/cmappy.pxd
index 5a3927e..1bbe749 100644
--- a/python/cmappy.pxd
+++ b/python/cmappy.pxd
@@ -13,10 +13,11 @@ cdef extern from "minimap.h":
int seed
int sdust_thres
int flag
+ int max_qlen
int bw
int max_gap, max_gap_ref
int max_frag_len
- int max_chain_skip
+ int max_chain_skip, max_chain_iter
int min_cnt
int min_chain_score
float mask_level
@@ -24,7 +25,7 @@ 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
diff --git a/python/mappy.pyx b/python/mappy.pyx
index b9e703f..2528d9d 100644
--- a/python/mappy.pyx
+++ b/python/mappy.pyx
@@ -3,7 +3,7 @@ from libc.stdlib cimport free
cimport cmappy
import sys
-__version__ = '2.14'
+__version__ = '2.16'
cmappy.mm_reset_timer()
@@ -221,6 +221,16 @@ 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))
diff --git a/setup.py b/setup.py
index 1760057..b61a5e2 100644
--- a/setup.py
+++ b/setup.py
@@ -33,7 +33,7 @@ def readme():
setup(
name = 'mappy',
- version = '2.14',
+ version = '2.16',
url = 'https://github.com/lh3/minimap2',
description = 'Minimap2 python binding',
long_description = readme(),
diff --git a/splitidx.c b/splitidx.c
index 6d12071..1054a64 100644
--- a/splitidx.c
+++ b/splitidx.c
@@ -11,8 +11,11 @@ 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);
- fp = fopen(fn, "wb");
- assert(fp);
+ if ((fp = fopen(fn, "wb")) == NULL) {
+ if (mm_verbose >= 1)
+ 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) {