diff --git a/NEWS.md b/NEWS.md index 3a8a06f..e36db68 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,26 @@ +Release 2.13-r850 (11 October 2018) +----------------------------------- + +Changes to minimap2: + + * Fixed wrongly formatted SAM when -L is in use (#231 and #233). + + * Fixed an integer overflow in rare cases. + + * Added --hard-mask-level to fine control split alignments (#244). + + * Made --MD work with spliced alignment (#139). + + * Replaced musl's getopt with ketopt for portability. + + * Log peak memory usage on exit. + +This release should produce alignments identical to v2.12 and v2.11. + +(2.13: 11 October 2018, r850) + + + Release 2.12-r827 (6 August 2018) --------------------------------- diff --git a/README.md b/README.md index b13136d..7ee0922 100644 --- a/README.md +++ b/README.md @@ -14,6 +14,7 @@ cd minimap2 && make # 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 +./minimap2 -ax asm20 ref.fa pacbio-ccs.fq.gz > aln.sam # PacBio CCS genomic reads ./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 @@ -70,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.12/minimap2-2.12_x64-linux.tar.bz2 | tar -jxvf - -./minimap2-2.12_x64-linux/minimap2 +curl -L https://github.com/lh3/minimap2/releases/download/v2.13/minimap2-2.13_x64-linux.tar.bz2 | tar -jxvf - +./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 diff --git a/align.c b/align.c index cea3618..84ef1ef 100644 --- a/align.c +++ b/align.c @@ -593,9 +593,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 diff --git a/cookbook.md b/cookbook.md index fa5f14e..ea77817 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.12/minimap2-2.12_x64-linux.tar.bz2 | tar jxf - -cp minimap2-2.12_x64-linux/{minimap2,k8,paftools.js} . # copy executables +curl -L https://github.com/lh3/minimap2/releases/download/v2.13/minimap2-2.13_x64-linux.tar.bz2 | tar jxf - +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 - diff --git a/hit.c b/hit.c index 9616bf5..92aa6b7 100644 --- a/hit.c +++ b/hit.c @@ -106,7 +106,7 @@ void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a) r->split |= 1, r2->split |= 2; } -void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff) // and compute mm_reg1_t::subsc +void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff, int hard_mask_level) // and compute mm_reg1_t::subsc { int i, j, k, *w; uint64_t *cov; @@ -118,6 +118,7 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff for (i = 1, k = 1; i < n; ++i) { mm_reg1_t *ri = &r[i]; int si = ri->qs, ei = ri->qe, n_cov = 0, uncov_len = 0; + if (hard_mask_level) goto skip_uncov; for (j = 0; j < k; ++j) { // traverse existing primary hits to find overlapping hits mm_reg1_t *rp = &r[w[j]]; int sj = rp->qs, ej = rp->qe; @@ -137,6 +138,7 @@ void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff } if (ei > x) uncov_len += ei - x; } +skip_uncov: for (j = 0; j < k; ++j) { // traverse existing primary hits again mm_reg1_t *rp = &r[w[j]]; int sj = rp->qs, ej = rp->qe, min, max, ol; diff --git a/main.c b/main.c index 4eb24b4..c3af456 100644 --- a/main.c +++ b/main.c @@ -6,7 +6,7 @@ #include "mmpriv.h" #include "ketopt.h" -#define MM_VERSION "2.12-r836-dirty" +#define MM_VERSION "2.13-r852-dirty" #ifdef __linux__ #include @@ -59,7 +59,8 @@ static ko_longopt_t long_options[] = { { "paf-no-hit", ko_no_argument, 333 }, { "split-prefix", ko_required_argument, 334 }, { "no-end-flt", ko_no_argument, 335 }, - { "bed", ko_required_argument, 336 }, + { "hard-mask-level",ko_no_argument, 336 }, + { "bed", ko_required_argument, 337 }, { "help", ko_no_argument, 'h' }, { "max-intron-len", ko_required_argument, 'G' }, { "version", ko_no_argument, 'V' }, @@ -97,11 +98,11 @@ 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:yY"; + 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; + int i, c, n_threads = 3, n_parts, old_best_n = -1; char *fnw = 0, *fn_bed = 0, *rg = 0, *s; FILE *fp_help = stderr; mm_idx_reader_t *idx_rdr; @@ -139,7 +140,7 @@ int main(int argc, char *argv[]) else if (c == 'g') opt.max_gap = (int)mm_parse_num(o.arg); else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(o.arg)); else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(o.arg); - else if (c == 'N') opt.best_n = atoi(o.arg); + else if (c == 'N') old_best_n = opt.best_n, opt.best_n = atoi(o.arg); else if (c == 'p') opt.pri_ratio = atof(o.arg); else if (c == 'M') opt.mask_level = atof(o.arg); else if (c == 'c') opt.flag |= MM_F_OUT_CG | MM_F_CIGAR; @@ -189,7 +190,8 @@ int main(int argc, char *argv[]) else if (c == 333) opt.flag |= MM_F_PAF_NO_HIT; // --paf-no-hit 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) fn_bed = o.arg; // --bed-prefer + else if (c == 336) opt.flag |= MM_F_HARD_MLEVEL; // --hard-mask-level + else if (c == 337) fn_bed = o.arg; // --bed-prefer else if (c == 314) { // --frag yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1); } else if (c == 315) { // --secondary @@ -253,6 +255,10 @@ int main(int argc, char *argv[]) ipt.flag |= MM_I_NO_SEQ; if (mm_check_opt(&ipt, &opt) < 0) return 1; + if (opt.best_n == 0) { + fprintf(stderr, "[WARNING]\033[1;31m changed '-N 0' to '-N %d --secondary=no'.\033[0m\n", old_best_n); + opt.best_n = old_best_n, opt.flag |= MM_F_NO_PRINT_2ND; + } if (argc == o.ind || fp_help == stdout) { fprintf(fp_help, "Usage: minimap2 [options] | [query.fa] [...]\n"); @@ -370,7 +376,7 @@ int main(int argc, char *argv[]) fprintf(stderr, "[M::%s] CMD:", __func__); for (i = 0; i < argc; ++i) fprintf(stderr, " %s", argv[i]); - fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec\n", __func__, realtime() - mm_realtime0, cputime()); + fprintf(stderr, "\n[M::%s] Real time: %.3f sec; CPU: %.3f sec; Peak RSS: %.3f GB\n", __func__, realtime() - mm_realtime0, cputime(), peakrss() / 1024.0 / 1024.0 / 1024.0); } return 0; } diff --git a/map.c b/map.c index 9c99824..6c21103 100644 --- a/map.c +++ b/map.c @@ -248,7 +248,7 @@ static mm128_t *collect_seed_hits(void *km, const mm_mapopt_t *opt, int max_occ, static void chain_post(const mm_mapopt_t *opt, int max_chain_gap_ref, const mm_idx_t *mi, void *km, int qlen, int n_segs, const int *qlens, int *n_regs, mm_reg1_t *regs, mm128_t *a) { if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s) - mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b); + mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL); if (n_segs <= 1) mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs); else mm_select_sub_multi(km, opt->pri_ratio, 0.2f, 0.7f, max_chain_gap_ref, mi->k*2, opt->best_n, n_segs, qlens, n_regs, regs); if (!(opt->flag & (MM_F_SPLICE|MM_F_SR|MM_F_NO_LJOIN))) // long join not working well without primary chains @@ -261,7 +261,7 @@ static mm_reg1_t *align_regs(const mm_mapopt_t *opt, const mm_idx_t *mi, void *k if (!(opt->flag & MM_F_CIGAR)) return regs; regs = mm_align_skeleton(km, opt, mi, qlen, seq, n_regs, regs, a); // this calls mm_filter_regs() if (!(opt->flag & MM_F_ALL_CHAINS)) { // don't choose primary mapping(s) - mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b); + mm_set_parent(km, opt->mask_level, *n_regs, regs, opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL); mm_select_sub(km, opt->pri_ratio, mi->k*2, opt->best_n, n_regs, regs); mm_set_sam_pri(*n_regs, regs); } @@ -359,7 +359,7 @@ void mm_map_frag(const mm_idx_t *mi, int n_segs, const int *qlens, const char ** seg = mm_seg_gen(b->km, hash, n_segs, qlens, n_regs0, regs0, n_regs, regs, a); // split fragment chain to separate segment chains free(regs0); for (i = 0; i < n_segs; ++i) { - mm_set_parent(b->km, opt->mask_level, n_regs[i], regs[i], opt->a * 2 + opt->b); // update mm_reg1_t::parent + mm_set_parent(b->km, opt->mask_level, n_regs[i], regs[i], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL); // update mm_reg1_t::parent regs[i] = align_regs(opt, mi, b->km, qlens[i], seqs[i], &n_regs[i], regs[i], seg[i].a); mm_set_mapq(b->km, n_regs[i], regs[i], opt->min_chain_score, opt->a, rep_len, is_sr); } @@ -502,7 +502,7 @@ static void merge_hits(step_t *s) } } mm_hit_sort(km, &s->n_reg[k], s->reg[k]); - mm_set_parent(km, opt->mask_level, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b); + mm_set_parent(km, opt->mask_level, s->n_reg[k], s->reg[k], opt->a * 2 + opt->b, opt->flag&MM_F_HARD_MLEVEL); if (!(opt->flag & MM_F_ALL_CHAINS)) { mm_select_sub(km, opt->pri_ratio, s->p->mi->k*2, opt->best_n, &s->n_reg[k], s->reg[k]); mm_set_sam_pri(s->n_reg[k], s->reg[k]); diff --git a/minimap.h b/minimap.h index 50e084f..231a9cf 100644 --- a/minimap.h +++ b/minimap.h @@ -34,6 +34,7 @@ #define MM_F_EQX 0x4000000 // use =/X instead of M #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_I_HPC 0x1 #define MM_I_NO_SEQ 0x2 diff --git a/minimap2.1 b/minimap2.1 index cea2e6e..879ef44 100644 --- a/minimap2.1 +++ b/minimap2.1 @@ -1,4 +1,4 @@ -.TH minimap2 1 "6 August 2018" "minimap2-2.12 (r827)" "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 @@ -193,7 +193,7 @@ Primarily used for all-vs-all read overlapping. .BI -p \ FLOAT Minimal secondary-to-primary score ratio to output secondary mappings [0.8]. Between two chains overlaping over half of the shorter chain (controlled by -.BR --mask-level ), +.BR -M ), the chain with a lower score is secondary to the chain with a higher score. If the ratio of the scores is below .IR FLOAT , @@ -226,6 +226,11 @@ Mark as secondary a chain that overlaps with a better chain by .I FLOAT or more of the shorter chain [0.5] .TP +.B --hard-mask-level +Honor option +.B -M +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 for chaining. The time complexity is quadratic in the number of seeds. This diff --git a/misc/paftools.js b/misc/paftools.js index 8cfe461..e489f37 100755 --- a/misc/paftools.js +++ b/misc/paftools.js @@ -1,6 +1,6 @@ #!/usr/bin/env k8 -var paftools_version = 'r767'; +var paftools_version = '2.13-r850'; /***************************** ***** Library functions ***** @@ -341,14 +341,15 @@ function paf_call(args) { var re_cs = /([:=*+-])(\d+|[A-Za-z]+)/g, re_tag = /\t(\S\S:[AZif]):(\S+)/g; var c, min_cov_len = 10000, min_var_len = 50000, gap_thres = 50, gap_thres_long = 1000, min_mapq = 5; - var fa_tmp = null, fa, fa_lens, is_vcf = false; - while ((c = getopt(args, "l:L:g:q:B:f:")) != null) { + var fa_tmp = null, fa, fa_lens, is_vcf = false, sample_name = "sample"; + while ((c = getopt(args, "l:L:g:q:B:f:s:")) != null) { if (c == 'l') min_cov_len = parseInt(getopt.arg); else if (c == 'L') min_var_len = parseInt(getopt.arg); else if (c == 'g') gap_thres = parseInt(getopt.arg); else if (c == 'G') gap_thres_long = parseInt(getopt.arg); else if (c == 'q') min_mapq = parseInt(getopt.arg); else if (c == 'f') fa_tmp = fasta_read(getopt.arg, fa_lens); + else if (c == 's') sample_name = getopt.arg; } if (fa_tmp != null) fa = fa_tmp[0], fa_lens = fa_tmp[1], is_vcf = true; @@ -360,6 +361,7 @@ function paf_call(args) print(" -q INT min mapping quality ["+min_mapq+"]"); print(" -g INT short/long gap threshold (for statistics only) ["+gap_thres+"]"); print(" -f FILE reference sequences (enabling VCF output) [null]"); + print(" -s NAME sample name in VCF header ["+sample_name+"]"); exit(1); } @@ -423,7 +425,7 @@ function paf_call(args) print('##INFO='); print('##INFO='); print('##FORMAT='); - print('#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT sample'); + print('#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT '+sample_name); } var a = [], out = []; diff --git a/mmpriv.h b/mmpriv.h index 2e56bb0..e68e996 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -83,7 +83,7 @@ void mm_split_reg(mm_reg1_t *r, mm_reg1_t *r2, int n, int qlen, mm128_t *a); void mm_sync_regs(void *km, int n_regs, mm_reg1_t *regs); int mm_squeeze_a(void *km, int n_regs, mm_reg1_t *regs, mm128_t *a); int mm_set_sam_pri(int n, mm_reg1_t *r); -void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff); +void mm_set_parent(void *km, float mask_level, int n, mm_reg1_t *r, int sub_diff, int hard_mask_level); void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int *n_, mm_reg1_t *r); void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int max_gap_ref, int min_diff, int best_n, int n_segs, const int *qlens, int *n_, mm_reg1_t *r); void mm_filter_regs(const mm_mapopt_t *opt, int qlen, int *n_regs, mm_reg1_t *regs); diff --git a/options.c b/options.c index d43ebec..c7f0eeb 100644 --- a/options.c +++ b/options.c @@ -174,5 +174,10 @@ int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo) fprintf(stderr, "[ERROR]\033[1;31m Z-drop should not be less than inversion-Z-drop\033[0m\n"); return -5; } + if ((mo->flag & MM_F_NO_PRINT_2ND) && (mo->flag & MM_F_ALL_CHAINS)) { + if (mm_verbose >= 1) + fprintf(stderr, "[ERROR]\033[1;31m -X/-P and --secondary=no can't be applied at the same time\033[0m\n"); + return -5; + } return 0; } diff --git a/python/mappy.pyx b/python/mappy.pyx index 9e7f842..3918f42 100644 --- a/python/mappy.pyx +++ b/python/mappy.pyx @@ -3,7 +3,7 @@ from libc.stdlib cimport free cimport cmappy import sys -__version__ = '2.11' +__version__ = '2.13' cmappy.mm_reset_timer() diff --git a/setup.py b/setup.py index 4d8a90d..6614571 100644 --- a/setup.py +++ b/setup.py @@ -33,7 +33,7 @@ def readme(): setup( name = 'mappy', - version = '2.12', + version = '2.13', url = 'https://github.com/lh3/minimap2', description = 'Minimap2 python binding', long_description = readme(),