From cf2bae6e9b5ff709c640e8ea56f9b22dbbdf311f Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 10 Jan 2019 09:57:05 -0500 Subject: [PATCH 01/20] r904: fixed a corner-case segfault. Resolves #307. --- align.c | 4 +++- main.c | 2 +- 2 files changed, 4 insertions(+), 2 deletions(-) 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/main.c b/main.c index a7a7fea..9f52bb7 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.14-r894-dirty" +#define MM_VERSION "2.14-r904-dirty" #ifdef __linux__ #include From c404f49569fa2d606b652418ffa4b9743bcaf641 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 10 Jan 2019 12:34:45 -0500 Subject: [PATCH 02/20] Release minimap2-2.15 (r905) --- NEWS.md | 36 ++++++++++++++++++++++++++++++++++++ README.md | 4 ++-- cookbook.md | 4 ++-- main.c | 2 +- minimap2.1 | 2 +- misc/paftools.js | 2 +- python/mappy.pyx | 2 +- setup.py | 2 +- 8 files changed, 45 insertions(+), 9 deletions(-) diff --git a/NEWS.md b/NEWS.md index e5060bc..e540e50 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,39 @@ +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..3fe3387 100644 --- a/README.md +++ b/README.md @@ -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.15/minimap2-2.15_x64-linux.tar.bz2 | tar -jxvf - +./minimap2-2.15_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 diff --git a/cookbook.md b/cookbook.md index d43d190..6345f0e 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.15/minimap2-2.15_x64-linux.tar.bz2 | tar jxf - +cp minimap2-2.15_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/main.c b/main.c index 9f52bb7..9481685 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.14-r904-dirty" +#define MM_VERSION "2.15-r905" #ifdef __linux__ #include diff --git a/minimap2.1 b/minimap2.1 index b9badbb..ffb8158 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 "10 January 2019" "minimap2-2.15 (r905)" "Bioinformatics tools" .SH NAME .PP minimap2 - mapping and alignment between collections of DNA sequences diff --git a/misc/paftools.js b/misc/paftools.js index c84ffa8..77e3c94 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.15-r905'; /***************************** ***** Library functions ***** diff --git a/python/mappy.pyx b/python/mappy.pyx index b9e703f..8d55fca 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.15' cmappy.mm_reset_timer() diff --git a/setup.py b/setup.py index 1760057..bf689f4 100644 --- a/setup.py +++ b/setup.py @@ -33,7 +33,7 @@ def readme(): setup( name = 'mappy', - version = '2.14', + version = '2.15', url = 'https://github.com/lh3/minimap2', description = 'Minimap2 python binding', long_description = readme(), From 48e230f40d121a0c40fc5a49480be9223dc4bc7a Mon Sep 17 00:00:00 2001 From: Heng Li Date: Fri, 11 Jan 2019 19:39:09 -0500 Subject: [PATCH 03/20] r906: de tag is wrongly calculated given "N" Resolves #309 --- format.c | 2 +- main.c | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/format.c b/format.c index e81ef32..69e25e0 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) diff --git a/main.c b/main.c index 9481685..02373a7 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.15-r905" +#define MM_VERSION "2.15-r906-dirty" #ifdef __linux__ #include From 30abcf3cf95c4e23345b90f12ad2502d4ba657d2 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sun, 13 Jan 2019 17:52:31 -0500 Subject: [PATCH 04/20] r907: copy tag "cs" in sam2paf Resolves #310 --- misc/paftools.js | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/misc/paftools.js b/misc/paftools.js index 77e3c94..468111d 100755 --- a/misc/paftools.js +++ b/misc/paftools.js @@ -1,6 +1,6 @@ #!/usr/bin/env k8 -var paftools_version = '2.15-r905'; +var paftools_version = '2.15-r907-dirty'; /***************************** ***** Library functions ***** @@ -1623,13 +1623,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,7 +1666,7 @@ function paf_sam2paf(args) } // parse MD var cs = []; - if (MD != null) { + if (MD != null && cs_str == null) { 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 @@ -1731,7 +1732,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")); From 597212b9f37ea311a67d5e57495658fc1b752af7 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 23 Jan 2019 11:18:50 -0500 Subject: [PATCH 05/20] r908: added an assertion to detect a potential bug as in #311 --- format.c | 1 + main.c | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/format.c b/format.c index 69e25e0..e300fe7 100644 --- a/format.c +++ b/format.c @@ -368,6 +368,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]); diff --git a/main.c b/main.c index 02373a7..19359f6 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.15-r906-dirty" +#define MM_VERSION "2.15-r908-dirty" #ifdef __linux__ #include From a8e1e3cbb879f7c73dce7f550af7745dfed5e911 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sat, 26 Jan 2019 17:59:36 -0500 Subject: [PATCH 06/20] updated citation with page numbers --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 3fe3387..8430951 100644 --- a/README.md +++ b/README.md @@ -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 From eba237f39d549bd4e842e7d2c3f2c3797b2bef0c Mon Sep 17 00:00:00 2001 From: Heng Li Date: Tue, 29 Jan 2019 09:20:23 -0500 Subject: [PATCH 07/20] r910: meaningful error message (#320) when minimap2 fails to create temporary files --- main.c | 2 +- splitidx.c | 7 +++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/main.c b/main.c index 19359f6..1421959 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.15-r908-dirty" +#define MM_VERSION "2.15-r910-dirty" #ifdef __linux__ #include diff --git a/splitidx.c b/splitidx.c index 6d12071..478a76e 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, "[E::%s] failed to create file '%s'\n", __func__, 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) { From 8b05880f73c76ccc00d2d1aa62eb3c999ca1c130 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Tue, 29 Jan 2019 10:42:20 -0500 Subject: [PATCH 08/20] r911: option -o to output to file (#319) --- main.c | 14 +++++++++++--- minimap2.1 | 7 ++++++- splitidx.c | 2 +- 3 files changed, 18 insertions(+), 5 deletions(-) diff --git a/main.c b/main.c index 1421959..69c45b8 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.15-r910-dirty" +#define MM_VERSION "2.15-r911-dirty" #ifdef __linux__ #include @@ -99,7 +99,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,6 +165,14 @@ 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 @@ -293,7 +301,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/minimap2.1 b/minimap2.1 index ffb8158..45ff89c 100644 --- a/minimap2.1 +++ b/minimap2.1 @@ -1,4 +1,4 @@ -.TH minimap2 1 "10 January 2019" "minimap2-2.15 (r905)" "Bioinformatics tools" +.TH minimap2 1 "10 January 2019" "minimap2-2.15-dirty (r911)" "Bioinformatics tools" .SH NAME .PP minimap2 - mapping and alignment between collections of DNA sequences @@ -384,6 +384,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 diff --git a/splitidx.c b/splitidx.c index 478a76e..1054a64 100644 --- a/splitidx.c +++ b/splitidx.c @@ -13,7 +13,7 @@ FILE *mm_split_init(const char *prefix, const mm_idx_t *mi) sprintf(fn, "%s.%.4d.tmp", prefix, mi->index); if ((fp = fopen(fn, "wb")) == NULL) { if (mm_verbose >= 1) - fprintf(stderr, "[E::%s] failed to create file '%s'\n", __func__, fn); + 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); From aaf32338184f71deb50a6f78da7a79dc06583feb Mon Sep 17 00:00:00 2001 From: Heng Li Date: Tue, 29 Jan 2019 12:53:20 -0500 Subject: [PATCH 09/20] added mappy.Aligner.seq_names to return seq names Resolves #312 --- python/README.rst | 6 ++++++ python/mappy.pyx | 10 ++++++++++ 2 files changed, 16 insertions(+) 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/mappy.pyx b/python/mappy.pyx index 8d55fca..7b0667f 100644 --- a/python/mappy.pyx +++ b/python/mappy.pyx @@ -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)) From 1739a260fb044eded97b3271cd5204c79363fa76 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Tue, 5 Feb 2019 14:19:17 -0500 Subject: [PATCH 10/20] r913: output tag "rl", length of unseedable regs --- format.c | 17 +++++++++++++++-- main.c | 2 +- map.c | 8 ++++---- minimap2.1 | 3 ++- mmpriv.h | 2 ++ 5 files changed, 24 insertions(+), 8 deletions(-) diff --git a/format.c b/format.c index e300fe7..2379a00 100644 --- a/format.c +++ b/format.c @@ -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); + 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]; @@ -377,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; @@ -528,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); @@ -535,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/main.c b/main.c index 69c45b8..798f6d7 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.15-r911-dirty" +#define MM_VERSION "2.15-r913-dirty" #ifdef __linux__ #include diff --git a/map.c b/map.c index 16fbe06..0f1bac2 100644 --- a/map.c +++ b/map.c @@ -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/minimap2.1 b/minimap2.1 index 45ff89c..016af5e 100644 --- a/minimap2.1 +++ b/minimap2.1 @@ -1,4 +1,4 @@ -.TH minimap2 1 "10 January 2019" "minimap2-2.15-dirty (r911)" "Bioinformatics tools" +.TH minimap2 1 "5 Feburary 2019" "minimap2-2.15-dirty (r913)" "Bioinformatics tools" .SH NAME .PP minimap2 - mapping and alignment between collections of DNA sequences @@ -617,6 +617,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/mmpriv.h b/mmpriv.h index 2688bca..a059701 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -61,8 +61,10 @@ 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); From 32ab6ce15b6d94dd66d7ad97e3b866c751e71d5b Mon Sep 17 00:00:00 2001 From: Heng Li Date: Tue, 12 Feb 2019 19:30:49 -0500 Subject: [PATCH 11/20] r914: fixed two harmless division by 0 Resolves #326 --- chain.c | 1 + hit.c | 1 + main.c | 2 +- 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/chain.c b/chain.c index 4947a3e..d8c1eae 100644 --- a/chain.c +++ b/chain.c @@ -28,6 +28,7 @@ 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) return 0; f = (int32_t*)kmalloc(km, n * 4); p = (int32_t*)kmalloc(km, n * 4); t = (int32_t*)kmalloc(km, n * 4); 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/main.c b/main.c index 798f6d7..970884e 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.15-r913-dirty" +#define MM_VERSION "2.15-r914-dirty" #ifdef __linux__ #include From e02bebd96d703707e30a85520d5e2f4bf84835b1 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 14 Feb 2019 10:04:04 -0500 Subject: [PATCH 12/20] r915: fixed a bug caused by the latest change --- chain.c | 5 ++++- main.c | 2 +- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/chain.c b/chain.c index d8c1eae..f7c05cd 100644 --- a/chain.c +++ b/chain.c @@ -28,7 +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) return 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); diff --git a/main.c b/main.c index 970884e..fac05b1 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.15-r914-dirty" +#define MM_VERSION "2.15-r915-dirty" #ifdef __linux__ #include From 19208fb06b288ad43cf51f4072bb7cea886030c8 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Sun, 17 Feb 2019 09:35:23 -0500 Subject: [PATCH 13/20] r916: support long cs in sam-to-paf conversion --- misc/paftools.js | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/misc/paftools.js b/misc/paftools.js index 468111d..75ed5bd 100755 --- a/misc/paftools.js +++ b/misc/paftools.js @@ -1,6 +1,6 @@ #!/usr/bin/env k8 -var paftools_version = '2.15-r907-dirty'; +var paftools_version = '2.15-r916-dirty'; /***************************** ***** Library functions ***** @@ -1590,11 +1590,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); } @@ -1666,8 +1671,8 @@ function paf_sam2paf(args) } // parse MD var cs = []; - if (MD != null && cs_str == 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; @@ -1681,13 +1686,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; } From ea84fc0a530aec95f8bed699d0cedc0a6ec3b2e1 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 27 Feb 2019 11:21:57 -0500 Subject: [PATCH 14/20] r917: fixed a bug in command-line parsing Resolves #344 --- ketopt.h | 14 +++++++++----- main.c | 2 +- 2 files changed, 10 insertions(+), 6 deletions(-) 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 fac05b1..34c4696 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.15-r915-dirty" +#define MM_VERSION "2.15-r917-dirty" #ifdef __linux__ #include From ccf1680aaf927efa9baa37fb0f5176212484b87b Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 27 Feb 2019 12:43:33 -0500 Subject: [PATCH 15/20] make it explicit that -x is preferred for prebuilt --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 8430951..5857cdf 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 From d431dc0181ea63974ae16eb737eb0c0dee8a768a Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 27 Feb 2019 14:41:01 -0500 Subject: [PATCH 16/20] r917: added --max-chain-iter to avoid worst case Resolves #324 --- chain.c | 3 ++- main.c | 2 ++ map.c | 4 ++-- minimap.h | 2 +- minimap2.1 | 10 ++++++++-- mmpriv.h | 2 +- options.c | 1 + python/cmappy.pxd | 5 +++-- 8 files changed, 20 insertions(+), 9 deletions(-) diff --git a/chain.c b/chain.c index f7c05cd..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; @@ -49,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/main.c b/main.c index 34c4696..6ff3a4b 100644 --- a/main.c +++ b/main.c @@ -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' }, @@ -179,6 +180,7 @@ int main(int argc, char *argv[]) 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 diff --git a/map.c b/map.c index 0f1bac2..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; 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 016af5e..ca5f4ba 100644 --- a/minimap2.1 +++ b/minimap2.1 @@ -1,4 +1,4 @@ -.TH minimap2 1 "5 Feburary 2019" "minimap2-2.15-dirty (r913)" "Bioinformatics tools" +.TH minimap2 1 "27 Feburary 2019" "minimap2-2.15-dirty (r917)" "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 diff --git a/mmpriv.h b/mmpriv.h index a059701..c9b91cd 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -69,7 +69,7 @@ void mm_write_sam3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se 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/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 From 90b7b83ec7ec35c648cf7403aa232396fa43e110 Mon Sep 17 00:00:00 2001 From: Kevin Chan Date: Tue, 18 Dec 2018 17:52:20 -0800 Subject: [PATCH 17/20] fix typo in command line help --- main.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/main.c b/main.c index 6ff3a4b..e1a113c 100644 --- a/main.c +++ b/main.c @@ -278,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"); From 55e39c2d3027eae926feadded302258f442f1dca Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 27 Feb 2019 15:03:19 -0500 Subject: [PATCH 18/20] r921: output unmapped reads in full PAF --- format.c | 2 +- main.c | 2 +- minimap2.1 | 4 +++- misc/paftools.js | 20 +++++++++++--------- 4 files changed, 16 insertions(+), 12 deletions(-) diff --git a/format.c b/format.c index 2379a00..29138d1 100644 --- a/format.c +++ b/format.c @@ -305,7 +305,7 @@ void mm_write_paf3(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const { 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; } diff --git a/main.c b/main.c index e1a113c..ced4c0b 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.15-r917-dirty" +#define MM_VERSION "2.15-r921-dirty" #ifdef __linux__ #include diff --git a/minimap2.1 b/minimap2.1 index ca5f4ba..cfcfb1d 100644 --- a/minimap2.1 +++ b/minimap2.1 @@ -474,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 diff --git a/misc/paftools.js b/misc/paftools.js index 75ed5bd..4fe692e 100755 --- a/misc/paftools.js +++ b/misc/paftools.js @@ -1,6 +1,6 @@ #!/usr/bin/env k8 -var paftools_version = '2.15-r916-dirty'; +var paftools_version = '2.15-r921-dirty'; /***************************** ***** 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; From 6b391e3373fa9410f2bd30a31914c454353fdfe9 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 28 Feb 2019 15:49:24 -0500 Subject: [PATCH 19/20] Release minimap2-2.16 (r922) --- NEWS.md | 49 ++++++++++++++++++++++++++++++++++++++++++++++++ README.md | 4 ++-- cookbook.md | 4 ++-- main.c | 2 +- minimap2.1 | 2 +- misc/paftools.js | 2 +- python/mappy.pyx | 2 +- setup.py | 2 +- 8 files changed, 58 insertions(+), 9 deletions(-) diff --git a/NEWS.md b/NEWS.md index e540e50..7cb5401 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,52 @@ +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) ----------------------------------- diff --git a/README.md b/README.md index 5857cdf..11088d9 100644 --- a/README.md +++ b/README.md @@ -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.15/minimap2-2.15_x64-linux.tar.bz2 | tar -jxvf - -./minimap2-2.15_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 diff --git a/cookbook.md b/cookbook.md index 6345f0e..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.15/minimap2-2.15_x64-linux.tar.bz2 | tar jxf - -cp minimap2-2.15_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/main.c b/main.c index ced4c0b..aa0cf91 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.15-r921-dirty" +#define MM_VERSION "2.16-r922" #ifdef __linux__ #include diff --git a/minimap2.1 b/minimap2.1 index cfcfb1d..c605508 100644 --- a/minimap2.1 +++ b/minimap2.1 @@ -1,4 +1,4 @@ -.TH minimap2 1 "27 Feburary 2019" "minimap2-2.15-dirty (r917)" "Bioinformatics tools" +.TH minimap2 1 "28 Feburary 2019" "minimap2-2.16-dirty (r922)" "Bioinformatics tools" .SH NAME .PP minimap2 - mapping and alignment between collections of DNA sequences diff --git a/misc/paftools.js b/misc/paftools.js index 4fe692e..459eebd 100755 --- a/misc/paftools.js +++ b/misc/paftools.js @@ -1,6 +1,6 @@ #!/usr/bin/env k8 -var paftools_version = '2.15-r921-dirty'; +var paftools_version = '2.16-r922'; /***************************** ***** Library functions ***** diff --git a/python/mappy.pyx b/python/mappy.pyx index 7b0667f..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.15' +__version__ = '2.16' cmappy.mm_reset_timer() diff --git a/setup.py b/setup.py index bf689f4..b61a5e2 100644 --- a/setup.py +++ b/setup.py @@ -33,7 +33,7 @@ def readme(): setup( name = 'mappy', - version = '2.15', + version = '2.16', url = 'https://github.com/lh3/minimap2', description = 'Minimap2 python binding', long_description = readme(), From 169216bfff5721616763a468126aa1432d897ce9 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Thu, 28 Feb 2019 15:58:12 -0500 Subject: [PATCH 20/20] manpage was wrongly marked as "dirty" --- minimap2.1 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/minimap2.1 b/minimap2.1 index c605508..be1f9ed 100644 --- a/minimap2.1 +++ b/minimap2.1 @@ -1,4 +1,4 @@ -.TH minimap2 1 "28 Feburary 2019" "minimap2-2.16-dirty (r922)" "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