Compare commits

...
54 Commits
Author SHA1 Message Date
Heng Li 59f23f7579 Release minimap2-2.14 (r883) 2018-11-06 00:03:16 -05:00
Heng Li 5e55e397e9 r882: guard against -E0 (#263) 2018-11-05 23:36:12 -05:00
Heng Li 88c421e8de r881: a recent change reduces sr accuracy 2018-11-05 22:03:59 -05:00
Heng Li 3db5bfe6e5 r880: fixed false wrong FASTA/Q alert 2018-11-05 20:52:07 -05:00
Heng Li 83dfdd5f50 draft release note 2018-11-05 20:07:57 -05:00
Heng Li 8a2b1cd4c9 updated mappy for extra option max_sw_mat 2018-11-05 19:28:44 -05:00
Heng Li 1ede8ca170 r877: renamed cap-sw-mat to cap-sw-mem 2018-11-05 11:46:38 -05:00
Heng Li 13981404e2 r876: skip DP if taking too much RAM (#259) 2018-11-05 11:43:10 -05:00
Heng Li fd64dd26f6 r875: warn given incorrect FASTA/Q
resolves #252
resolves #255
2018-11-05 10:02:44 -05:00
Heng Li 24df95e4b8 r874: don't call x86_simd() so often
This takes a few percent of time in profiler.
2018-11-05 09:20:35 -05:00
Heng Li a8ee48c2ce r873: comforming to C99/C11; resolves #261 2018-11-05 08:25:07 -05:00
Heng Li 09e089c3dc r872: choose the longest isoform 2018-11-04 23:48:50 -05:00
Heng Li e46cbb7d84 r871: print erroneous genes 2018-11-04 20:37:06 -05:00
Heng Li 57ec73ec6c r870: separate <50% and <10% 2018-11-04 19:31:25 -05:00
Heng Li 9e27575387 r869: classify incomplete genes 2018-11-04 19:21:55 -05:00
Heng Li b4ad8d8bf0 added asmgene
improvements coming; not made public yet
2018-11-04 17:24:05 -05:00
Heng Li e315b9fada hidden options to control bp calculation 2018-11-04 16:36:04 -05:00
Heng Li 42baf287a4 r866: fixed a typo; resolves #262 2018-10-30 09:11:55 -04:00
Heng Li 2ceba22a7a fixed a typo in manpage 2018-10-28 11:51:02 -04:00
Heng Li 9ed56b4a25 r860: MD/cs not working with --eqx 2018-10-26 23:23:53 -04:00
Heng Li ecb6c5c36c Document --no-pairing (#256) 2018-10-23 10:00:21 -04:00
Heng Li 377c7099a8 r858: fixed a bug; resolves #254 2018-10-22 22:47:11 -04:00
Heng Li 51e2abfa60 clarify that minimap2 may miss small exons 2018-10-22 11:16:16 -04:00
Heng Li 7b0a49732e r856: wrongly reported for an unrecognized option
Resolved #250
2018-10-19 20:07:14 -04:00
Heng Li 20268a6068 updated the copyright holder 2018-10-18 11:11:17 -04:00
Heng Li d04ac068fd r852: a minor when large --end-bonus is in use
We may use a large --end-bonus to mimic end-to-end alignment. In the short-read
mode, the candidate alignment region may be out of the band, which leads to
truncated alignment.
2018-10-15 21:28:27 -04:00
Heng Li 5d5d392c02 Release minimap2-2.13 (r850) 2018-10-11 13:18:31 -04:00
Heng Li 170863e553 r849: option -P doesn't work
I don't know why I haven't found it at the beginning.
2018-10-04 16:11:59 -04:00
Heng Li 97f97306a4 r847: guard against -N0 2018-09-27 15:13:44 -04:00
Heng Li 1077b7ddc8 r846: added --hard-mask-level for #244 2018-09-27 14:46:26 -04:00
Heng Li c57b59f02f r845: log peak memory 2018-09-23 20:27:49 -04:00
Torsten Seemann 34be359e25 Add -s sample option for VCF output 2018-09-19 18:24:47 -04:00
Heng Li 8b12da8b0f removed a useless dependency (not in repo) 2018-09-17 13:59:25 -04:00
Heng Li c63a33904f r836: fixed an integer overflow
Forgot this one.
2018-09-14 23:29:31 -04:00
Heng Li 70b0fede64 r835: improved help message. Resolved #232 2018-09-14 22:29:25 -04:00
Heng Li 0b681e51e7 Merge remote-tracking branch 'remotes/origin/master' 2018-09-14 22:22:20 -04:00
Heng Li 7d80d6de4a r832: fixed outdated -L. Resolved #231 and #233 2018-09-14 22:21:33 -04:00
Chris Rands 791e89ce0f PEP8 and other minor styles changes 2018-09-12 10:58:56 -04:00
Heng Li 98c48a1c45 Fixed wrong version links in the cookbook 2018-09-03 09:07:34 -04:00
Heng Li 63d397120a removed getopt.c from MANIFEST 2018-09-01 21:30:46 -04:00
Heng Li 7998fe9906 r829: replaced musl's getopt with ketopt 2018-09-01 21:18:02 -04:00
Heng Li 3a119d606f r828: --MD to support spliced alignment 2018-08-22 10:47:45 -04:00
Heng Li a5eafb75f9 Release minimap2-2.12 (r827) 2018-08-06 12:44:39 -04:00
Heng Li 9a567e4b37 allow mappy to change scoring 2018-08-06 09:52:13 -04:00
Heng Li 8e606bcc06 added a -C5 example to Getting Started 2018-08-06 09:08:33 -04:00
Heng Li a1b7219b5d explain added parameters 2018-08-05 21:32:05 -04:00
Heng Li b0f39a1a61 r823: mappy to index a single sequence 2018-08-05 20:57:05 -04:00
Heng Li 5ab6538757 r822: added option --no-end-flt 2018-08-05 19:42:12 -04:00
Heng Li b32296e18f r821: fixed memory when -y is used 2018-07-31 15:14:37 -04:00
Heng Li 99ecdf7b5d mappy to support arm64 (#203) 2018-07-24 23:53:02 -04:00
Heng Li ff9917a1c4 r819: mappy to support cs/MD 2018-07-24 23:29:55 -04:00
Mark Bicknell 8c064a5f29 Fixed mm_idx_is_idx to return the correct result on Windows. 2018-07-17 09:16:33 -04:00
Heng Li 0e137670fc Merge branch 'split' 2018-07-15 22:24:15 -04:00
Heng Li f3417da838 bugfix: unmapped records are duplicated in output 2018-07-14 22:54:05 -04:00
32 changed files with 836 additions and 540 deletions
+2 -1
View File
@@ -1,6 +1,7 @@
The MIT License
Copyright (c) 2017 Broad Institute, Inc.
Copyright (c) 2018- Dana-Farber Cancer Institute
2017-2018 Broad Institute, Inc.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
+1 -1
View File
@@ -1,9 +1,9 @@
include *.h
include Makefile
include ksw2_dispatch.c
include getopt.c
include main.c
include README.md
include sse2neon/emmintrin.h
include python/mappy.c
include python/cmappy.h
include python/cmappy.pxd
+7 -7
View File
@@ -32,8 +32,8 @@ all:$(PROG)
extra:all $(PROG_EXTRA)
minimap2:main.o getopt.o libminimap2.a
$(CC) $(CFLAGS) main.o getopt.o -o $@ -L. -lminimap2 $(LIBS)
minimap2:main.o libminimap2.a
$(CC) $(CFLAGS) main.o -o $@ -L. -lminimap2 $(LIBS)
minimap2-lite:example.o libminimap2.a
$(CC) $(CFLAGS) $< -o $@ -L. -lminimap2 $(LIBS)
@@ -41,8 +41,8 @@ minimap2-lite:example.o libminimap2.a
libminimap2.a:$(OBJS)
$(AR) -csru $@ $(OBJS)
sdust:sdust.c getopt.o kalloc.o kalloc.h kdq.h kvec.h kseq.h sdust.h
$(CC) -D_SDUST_MAIN $(CFLAGS) $< getopt.o kalloc.o -o $@ -lz
sdust:sdust.c kalloc.o kalloc.h kdq.h kvec.h kseq.h ketopt.h sdust.h
$(CC) -D_SDUST_MAIN $(CFLAGS) $< kalloc.o -o $@ -lz
# SSE-specific targets on x86/x86_64
@@ -99,7 +99,6 @@ chain.o: minimap.h mmpriv.h bseq.h kalloc.h
esterr.o: mmpriv.h minimap.h bseq.h
example.o: minimap.h kseq.h
format.o: kalloc.h mmpriv.h minimap.h bseq.h
getopt.o: getopt.h
hit.o: mmpriv.h minimap.h bseq.h kalloc.h khash.h
index.o: kthread.h bseq.h minimap.h mmpriv.h kvec.h kalloc.h khash.h
kalloc.o: kalloc.h
@@ -108,11 +107,12 @@ ksw2_exts2_sse.o: ksw2.h kalloc.h
ksw2_extz2_sse.o: ksw2.h kalloc.h
ksw2_ll_sse.o: ksw2.h kalloc.h
kthread.o: kthread.h
main.o: bseq.h minimap.h mmpriv.h getopt.h
main.o: bseq.h minimap.h mmpriv.h ketopt.h
map.o: kthread.h kvec.h kalloc.h sdust.h mmpriv.h minimap.h bseq.h khash.h
map.o: ksort.h
misc.o: mmpriv.h minimap.h bseq.h ksort.h
options.o: mmpriv.h minimap.h bseq.h
pe.o: mmpriv.h minimap.h bseq.h kvec.h kalloc.h ksort.h
sdust.o: kalloc.h kdq.h kvec.h sdust.h
sdust.o: kalloc.h kdq.h kvec.h ketopt.h sdust.h
sketch.o: kvec.h kalloc.h mmpriv.h minimap.h bseq.h
splitidx.o: mmpriv.h minimap.h bseq.h
+73
View File
@@ -1,3 +1,76 @@
Release 2.14-r883 (5 November 2018)
-----------------------------------
Notable changes:
* Fixed two minor bugs caused by typos (#254 and #266).
* Fixed a bug that made minimap2 abort when --eqx was used together with --MD
or --cs (#257).
* Added --cap-sw-mem to cap the size of DP matrices (#259). Base alignment may
take a lot of memory in the splicing mode. This may lead to issues when we
run minimap2 on a cluster with a hard memory limit. The new option avoids
unlimited memory usage at the cost of missing a few long introns.
* Conforming to C99 and C11 when possible (#261).
* Warn about malformatted FASTA or FASTQ (#252 and #255).
This release occasionally produces base alignments different from v2.13. The
overall alignment accuracy remain similar.
(2.14: 5 November 2018, r883)
Release 2.13-r850 (11 October 2018)
-----------------------------------
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)
---------------------------------
Changes to minimap2:
* Added option --split-prefix to write proper alignments (correct mapping
quality and clustered query sequences) given a multi-part index (#141 and
#189; mostly by @hasindu2008).
* Fixed a memory leak when option -y is in use.
Changes to mappy:
* Support the MD/cs tag (#183 and #203).
* Allow mappy to index a single sequence, to add extra flags and to change the
scoring system.
Minimap2 should produce alignments identical to v2.11.
(2.12: 6 August 2018, r827)
Release 2.11-r797 (20 June 2018)
--------------------------------
+8 -4
View File
@@ -14,9 +14,11 @@ 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
./minimap2 -ax splice -k14 -uf ref.fa reads.fa > aln.sam # Nanopore Direct RNA-seq
./minimap2 -ax splice ref.fa rna-reads.fa > aln.sam # spliced long reads (strand unknown)
./minimap2 -ax splice -uf -k14 ref.fa reads.fa > aln.sam # noisy Nanopore Direct RNA-seq
./minimap2 -ax splice -uf -C5 ref.fa query.fa > aln.sam # Final PacBio Iso-seq or traditional cDNA
./minimap2 -cx asm5 asm1.fa asm2.fa > aln.paf # intra-species asm-to-asm alignment
./minimap2 -x ava-pb reads.fa reads.fa > overlaps.paf # PacBio read overlap
./minimap2 -x ava-ont reads.fa reads.fa > overlaps.paf # Nanopore read overlap
@@ -69,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.11/minimap2-2.11_x64-linux.tar.bz2 | tar -jxvf -
./minimap2-2.11_x64-linux/minimap2
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
```
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
@@ -353,6 +355,8 @@ mappy` or [from BioConda][mappyconda] via `conda install -c bioconda mappy`.
billion bases or longer (2,147,483,647 to be exact). The total length of all
sequences can well exceed this threshold.
* Minimap2 often misses small exons.
[paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
+11 -7
View File
@@ -300,7 +300,10 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint
for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr);
fputc('\n', stderr);
}
if (opt->flag & MM_F_SPLICE)
if (opt->max_sw_mat > 0 && (int64_t)tlen * qlen > opt->max_sw_mat) {
ksw_reset_extz(ez);
ez->zdropped = 1;
} else if (opt->flag & MM_F_SPLICE)
ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, flag, ez);
else if (opt->q == opt->q2 && opt->e == opt->e2)
ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, flag, ez);
@@ -414,7 +417,7 @@ static void mm_filter_bad_seeds_alt(void *km, int as1, int cnt1, mm128_t *a, int
gap2 = ((int32_t)a[as1 + j].y - (int32_t)a[as1 + j - 1].y) - (int32_t)(a[as1 + j].x - a[as1 + j - 1].x);
q_span_pre = a[as1 + j - 1].y >> 32 & 0xff;
rs2 = (int32_t)a[as1 + j - 1].x + q_span_pre;
qs2 = (int32_t)a[as1 + j - 1].x + q_span_pre;
qs2 = (int32_t)a[as1 + j - 1].y + q_span_pre;
m = rs2 - re1 < qs2 - qe1? rs2 - re1 : qs2 - qe1;
gap2 = gap2 > 0? gap2 : -gap2;
if (m > gap1 + gap2) break;
@@ -563,11 +566,12 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
re = (int32_t)a[as1+cnt1-1].x + 1;
qe = (int32_t)a[as1+cnt1-1].y + 1;
} else {
if (is_splice) {
mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1);
} else {
mm_fix_bad_ends(r, a, opt->bw, opt->min_chain_score * 2, &as1, &cnt1);
}
if (!(opt->flag & MM_F_NO_END_FLT)) {
if (is_splice)
mm_fix_bad_ends_splice(km, opt, mi, r, mat, qlen, qseq0, a, &as1, &cnt1);
else
mm_fix_bad_ends(r, a, opt->bw, opt->min_chain_score * 2, &as1, &cnt1);
} else as1 = r->as, cnt1 = r->cnt;
mm_filter_bad_seeds(km, as1, cnt1, a, 10, 40, opt->max_gap>>1, 10);
mm_filter_bad_seeds_alt(km, as1, cnt1, a, 30, opt->max_gap>>1);
mm_adjust_minier(mi, qseq0, &a[as1], &rs, &qs);
+7 -2
View File
@@ -39,7 +39,7 @@ mm_bseq_file_t *mm_bseq_open(const char *fn)
{
mm_bseq_file_t *fp;
gzFile f;
f = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(fileno(stdin), "r");
f = fn && strcmp(fn, "-")? gzopen(fn, "r") : gzdopen(0, "r");
if (f == 0) return 0;
fp = (mm_bseq_file_t*)calloc(1, sizeof(mm_bseq_file_t));
fp->fp = f;
@@ -65,6 +65,8 @@ static inline char *kstrdup(const kstring_t *s)
static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual, int with_comment)
{
int i;
if (ks->name.l == 0)
fprintf(stderr, "[WARNING]\033[1;31m empty sequence name in the input.\033[0m\n");
s->name = kstrdup(&ks->name);
s->seq = kstrdup(&ks->seq);
for (i = 0; i < (int)ks->seq.l; ++i) // convert U to T
@@ -78,6 +80,7 @@ static inline void kseq2bseq(kseq_t *ks, mm_bseq1_t *s, int with_qual, int with_
mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int with_comment, int frag_mode, int *n_)
{
int64_t size = 0;
int ret;
kvec_t(mm_bseq1_t) a = {0,0,0};
kseq_t *ks = fp->ks;
*n_ = 0;
@@ -87,7 +90,7 @@ mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
size = fp->s.l_seq;
memset(&fp->s, 0, sizeof(mm_bseq1_t));
}
while (kseq_read(ks) >= 0) {
while ((ret = kseq_read(ks)) >= 0) {
mm_bseq1_t *s;
assert(ks->seq.l <= INT32_MAX);
if (a.m == 0) kv_resize(mm_bseq1_t, 0, a, 256);
@@ -107,6 +110,8 @@ mm_bseq1_t *mm_bseq_read3(mm_bseq_file_t *fp, int chunk_size, int with_qual, int
break;
}
}
if (ret < -1)
fprintf(stderr, "[WARNING]\033[1;31m wrong FASTA/FASTQ record. Continue anyway.\033[0m\n");
*n_ = a.n;
return a.a;
}
+6 -6
View File
@@ -24,18 +24,18 @@
This cookbook walks you through a variety of applications of minimap2 and its
companion script `paftools.js`. All data here are freely available from the
minimap2 release page at version tag [v2.11][v2.11]. Some examples only work
with v2.11 or later.
minimap2 release page at version tag [v2.10][v2.10]. Some examples only work
with v2.10 or later.
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.11/minimap2-2.11_x64-linux.tar.bz2 | tar jxf -
cp minimap2-2.11_x64-linux/{minimap2,k8,paftools.js} . # copy 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
export PATH="$PATH:"`pwd` # put the current directory on PATH
# download example datasets
curl -L https://github.com/lh3/minimap2/releases/download/v2.11/cookbook-data.tgz | tar zxf -
curl -L https://github.com/lh3/minimap2/releases/download/v2.10/cookbook-data.tgz | tar zxf -
```
## <a name="map-reads"></a>Mapping Genomic Reads
@@ -240,4 +240,4 @@ with `-x ava-pb` (99% vs 93% with `-x ava-ont`).
[pbsim]: https://github.com/pfaucon/PBSIM-PacBio-Simulator
[mason2]: https://github.com/seqan/seqan/tree/master/apps/mason2
[paf]: https://github.com/lh3/miniasm/blob/master/PAF.md
[v2.11]: https://github.com/lh3/minimap2/releases/tag/v2.11
[v2.10]: https://github.com/lh3/minimap2/releases/tag/v2.10
+45 -17
View File
@@ -92,7 +92,8 @@ static void sam_write_rg_line(kstring_t *str, const char *s)
if (mm_verbose >= 1) fprintf(stderr, "[ERROR] the read group line contained literal <tab> characters -- replace with escaped tabs: \\t\n");
goto err_set_rg;
}
rg_line = strdup(s);
rg_line = (char*)malloc(strlen(s) + 1);
strcpy(rg_line, s);
mm_escape(rg_line);
if ((p = strstr(rg_line, "\tID:")) == 0) {
if (mm_verbose >= 1) fprintf(stderr, "[ERROR] no ID within the read group line\n");
@@ -133,14 +134,14 @@ void mm_write_sam_hdr(const mm_idx_t *idx, const char *rg, const char *ver, int
free(str.s);
}
static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp, int no_iden)
static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp, int no_iden, int write_tag)
{
int i, q_off, t_off;
mm_sprintf_lite(s, "\tcs:Z:");
if (write_tag) mm_sprintf_lite(s, "\tcs:Z:");
for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
assert(op >= 0 && op <= 3);
if (op == 0) { // match
assert((op >= 0 && op <= 3) || op == 7 || op == 8);
if (op == 0 || op == 7 || op == 8) { // match
int l_tmp = 0;
for (j = 0; j < len; ++j) {
if (qseq[q_off + j] != tseq[t_off + j]) {
@@ -181,14 +182,14 @@ static void write_cs_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
assert(t_off == r->re - r->rs && q_off == r->qe - r->qs);
}
static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp)
static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq, const mm_reg1_t *r, char *tmp, int write_tag)
{
int i, q_off, t_off, l_MD = 0;
mm_sprintf_lite(s, "\tMD:Z:");
if (write_tag) mm_sprintf_lite(s, "\tMD:Z:");
for (i = q_off = t_off = 0; i < (int)r->p->n_cigar; ++i) {
int j, op = r->p->cigar[i]&0xf, len = r->p->cigar[i]>>4;
assert(op >= 0 && op <= 2); // introns (aka reference skips) are not supported
if (op == 0) { // match
assert((op >= 0 && op <= 3) || op == 7 || op == 8);
if (op == 0 || op == 7 || op == 8) { // match
for (j = 0; j < len; ++j) {
if (qseq[q_off + j] != tseq[t_off + j]) {
mm_sprintf_lite(s, "%d%c", l_MD, "ACGTN"[tseq[t_off + j]]);
@@ -204,13 +205,15 @@ static void write_MD_core(kstring_t *s, const uint8_t *tseq, const uint8_t *qseq
mm_sprintf_lite(s, "%d^%s", l_MD, tmp);
l_MD = 0;
t_off += len;
} else if (op == 3) { // reference skip
t_off += len;
}
}
if (l_MD > 0) mm_sprintf_lite(s, "%d", l_MD);
assert(t_off == r->re - r->rs && q_off == r->qe - r->qs);
}
static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int no_iden, int is_MD)
static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const mm_reg1_t *r, int no_iden, int is_MD, int write_tag)
{
extern unsigned char seq_nt4_table[256];
int i;
@@ -230,11 +233,34 @@ static void write_cs_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const mm_
qseq[r->qe - i - 1] = c >= 4? 4 : 3 - c;
}
}
if (is_MD) write_MD_core(s, tseq, qseq, r, tmp);
else write_cs_core(s, tseq, qseq, r, tmp, no_iden);
if (is_MD) write_MD_core(s, tseq, qseq, r, tmp, write_tag);
else write_cs_core(s, tseq, qseq, r, tmp, no_iden, write_tag);
kfree(km, qseq); kfree(km, tseq); kfree(km, tmp);
}
int mm_gen_cs_or_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int is_MD, int no_iden)
{
mm_bseq1_t t;
kstring_t str;
str.s = *buf, str.l = 0, str.m = *max_len;
t.l_seq = strlen(seq);
t.seq = (char*)seq;
write_cs_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0);
*max_len = str.m;
*buf = str.s;
return str.l;
}
int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden)
{
return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 0, no_iden);
}
int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq)
{
return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0);
}
static inline void write_tags(kstring_t *s, const mm_reg1_t *r)
{
int type;
@@ -277,7 +303,7 @@ void mm_write_paf(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, const m
mm_sprintf_lite(s, "%d%c", r->p->cigar[k]>>4, "MIDNSHP=XB"[r->p->cigar[k]&0xf]);
}
if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_MD)))
write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD);
write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1);
if ((opt_flag & MM_F_COPY_COMMENT) && t->comment)
mm_sprintf_lite(s, "\t%s", t->comment);
}
@@ -394,9 +420,11 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
cigar_in_tag = 1;
}
if (cigar_in_tag) {
if (flag & 0x100) mm_sprintf_lite(s, "0S"); // secondary alignment
else if (flag & 0x800) mm_sprintf_lite(s, "%dS", r->re - r->rs); // supplementary alignment
else mm_sprintf_lite(s, "%dS", t->l_seq);
int slen;
if ((flag & 0x900) == 0 || (opt_flag & MM_F_SOFTCLIP)) slen = t->l_seq;
else if (flag & 0x100) slen = 0;
else slen = r->qe - r->qs;
mm_sprintf_lite(s, "%dS%dN", slen, r->re - r->rs);
} else write_sam_cigar(s, flag, 0, t->l_seq, r, opt_flag);
}
@@ -476,7 +504,7 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se
}
}
if (r->p && (opt_flag & (MM_F_OUT_CS|MM_F_OUT_MD)))
write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD);
write_cs_or_MD(km, s, mi, t, r, !(opt_flag&MM_F_OUT_CS_LONG), opt_flag&MM_F_OUT_MD, 1);
if (cigar_in_tag)
write_sam_cigar(s, flag, 1, t->l_seq, r, opt_flag);
}
-216
View File
@@ -1,216 +0,0 @@
#include <stddef.h>
#include <stdio.h>
#include <string.h>
#include "getopt.h"
char *optarg;
int optind=1, opterr=1, optopt, __optpos, optreset=0;
#define optpos __optpos
static void __getopt_msg(const char *a, const char *b, const char *c, size_t l)
{
FILE *f = stderr;
#if !defined(WIN32) && !defined(_WIN32)
flockfile(f);
#endif
fputs(a, f);
fwrite(b, strlen(b), 1, f);
fwrite(c, 1, l, f);
fputc('\n', f);
#if !defined(WIN32) && !defined(_WIN32)
funlockfile(f);
#endif
}
int getopt(int argc, char * const argv[], const char *optstring)
{
int i, c, d;
int k, l;
char *optchar;
if (!optind || optreset) {
optreset = 0;
__optpos = 0;
optind = 1;
}
if (optind >= argc || !argv[optind])
return -1;
if (argv[optind][0] != '-') {
if (optstring[0] == '-') {
optarg = argv[optind++];
return 1;
}
return -1;
}
if (!argv[optind][1])
return -1;
if (argv[optind][1] == '-' && !argv[optind][2])
return optind++, -1;
if (!optpos) optpos++;
c = argv[optind][optpos], k = 1;
optchar = argv[optind]+optpos;
optopt = c;
optpos += k;
if (!argv[optind][optpos]) {
optind++;
optpos = 0;
}
if (optstring[0] == '-' || optstring[0] == '+')
optstring++;
i = 0;
d = 0;
do {
d = optstring[i], l = 1;
if (l>0) i+=l; else i++;
} while (l && d != c);
if (d != c) {
if (optstring[0] != ':' && opterr)
__getopt_msg(argv[0], ": unrecognized option: ", optchar, k);
return '?';
}
if (optstring[i] == ':') {
if (optstring[i+1] == ':') optarg = 0;
else if (optind >= argc) {
if (optstring[0] == ':') return ':';
if (opterr) __getopt_msg(argv[0],
": option requires an argument: ",
optchar, k);
return '?';
}
if (optstring[i+1] != ':' || optpos) {
optarg = argv[optind++] + optpos;
optpos = 0;
}
}
return c;
}
static void permute(char *const *argv, int dest, int src)
{
char **av = (char **)argv;
char *tmp = av[src];
int i;
for (i=src; i>dest; i--)
av[i] = av[i-1];
av[dest] = tmp;
}
static int __getopt_long_core(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx, int longonly)
{
optarg = 0;
if (longopts && argv[optind][0] == '-' &&
((longonly && argv[optind][1] && argv[optind][1] != '-') ||
(argv[optind][1] == '-' && argv[optind][2])))
{
int colon = optstring[optstring[0]=='+'||optstring[0]=='-']==':';
int i, cnt, match = -1;
char *opt;
for (cnt=i=0; longopts[i].name; i++) {
const char *name = longopts[i].name;
opt = argv[optind]+1;
if (*opt == '-') opt++;
for (; *name && *name == *opt; name++, opt++);
if (*opt && *opt != '=') continue;
match = i;
if (!*name) {
cnt = 1;
break;
}
cnt++;
}
if (cnt==1) {
i = match;
optind++;
optopt = longopts[i].val;
if (*opt == '=') {
if (!longopts[i].has_arg) {
if (colon || !opterr)
return '?';
__getopt_msg(argv[0],
": option does not take an argument: ",
longopts[i].name,
strlen(longopts[i].name));
return '?';
}
optarg = opt+1;
} else if (longopts[i].has_arg == required_argument) {
if (!(optarg = argv[optind])) {
if (colon) return ':';
if (!opterr) return '?';
__getopt_msg(argv[0],
": option requires an argument: ",
longopts[i].name,
strlen(longopts[i].name));
return '?';
}
optind++;
}
if (idx) *idx = i;
if (longopts[i].flag) {
*longopts[i].flag = longopts[i].val;
return 0;
}
return longopts[i].val;
}
if (argv[optind][1] == '-') {
if (!colon && opterr)
__getopt_msg(argv[0], cnt ?
": option is ambiguous: " :
": unrecognized option: ",
argv[optind]+2,
strlen(argv[optind]+2));
optind++;
return '?';
}
}
return getopt(argc, argv, optstring);
}
static int __getopt_long(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx, int longonly)
{
int ret, skipped, resumed;
if (!optind || optreset) {
optreset = 0;
__optpos = 0;
optind = 1;
}
if (optind >= argc || !argv[optind]) return -1;
skipped = optind;
if (optstring[0] != '+' && optstring[0] != '-') {
int i;
for (i=optind; ; i++) {
if (i >= argc || !argv[i]) return -1;
if (argv[i][0] == '-' && argv[i][1]) break;
}
optind = i;
}
resumed = optind;
ret = __getopt_long_core(argc, argv, optstring, longopts, idx, longonly);
if (resumed > skipped) {
int i, cnt = optind-resumed;
for (i=0; i<cnt; i++)
permute(argv, skipped, optind-1);
optind = skipped + cnt;
}
return ret;
}
int getopt_long(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx)
{
return __getopt_long(argc, argv, optstring, longopts, idx, 0);
}
int getopt_long_only(int argc, char *const *argv, const char *optstring, const struct option *longopts, int *idx)
{
return __getopt_long(argc, argv, optstring, longopts, idx, 1);
}
-53
View File
@@ -1,53 +0,0 @@
/*
Copyright 2005-2014 Rich Felker, et al.
Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be
included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
#ifndef _GETOPT_H
#define _GETOPT_H
#ifdef __cplusplus
extern "C" {
#endif
int getopt(int, char * const [], const char *);
extern char *optarg;
extern int optind, opterr, optopt, optreset;
struct option {
const char *name;
int has_arg;
int *flag;
int val;
};
int getopt_long(int, char *const *, const char *, const struct option *, int *);
int getopt_long_only(int, char *const *, const char *, const struct option *, int *);
#define no_argument 0
#define required_argument 1
#define optional_argument 2
#ifdef __cplusplus
}
#endif
#endif
+3 -1
View File
@@ -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;
+12 -1
View File
@@ -372,7 +372,9 @@ mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const cha
uint64_t sum_len = 0;
mm128_v a = {0,0,0};
mm_idx_t *mi;
khash_t(str) *h;
int i, flag = 0;
if (n <= 0) return 0;
for (i = 0; i < n; ++i) // get the total length
sum_len += strlen(seq[i]);
@@ -383,13 +385,17 @@ mm_idx_t *mm_idx_str(int w, int k, int is_hpc, int bucket_bits, int n, const cha
mi->n_seq = n;
mi->seq = (mm_idx_seq_t*)kcalloc(mi->km, n, sizeof(mm_idx_seq_t)); // ->seq is allocated from km
mi->S = (uint32_t*)calloc((sum_len + 7) / 8, 4);
mi->h = h = kh_init(str);
for (i = 0, sum_len = 0; i < n; ++i) {
const char *s = seq[i];
mm_idx_seq_t *p = &mi->seq[i];
uint32_t j;
if (name && name[i]) {
int absent;
p->name = (char*)kmalloc(mi->km, strlen(name[i]) + 1);
strcpy(p->name, name[i]);
kh_put(str, h, p->name, &absent);
assert(absent);
}
p->offset = sum_len;
p->len = strlen(s);
@@ -512,14 +518,19 @@ mm_idx_t *mm_idx_load(FILE *fp)
int64_t mm_idx_is_idx(const char *fn)
{
int fd, is_idx = 0;
off_t ret, off_end;
int64_t ret, off_end;
char magic[4];
if (strcmp(fn, "-") == 0) return 0; // read from pipe; not an index
fd = open(fn, O_RDONLY);
if (fd < 0) return -1; // error
#ifdef WIN32
if ((off_end = _lseeki64(fd, 0, SEEK_END)) >= 4) {
_lseeki64(fd, 0, SEEK_SET);
#else
if ((off_end = lseek(fd, 0, SEEK_END)) >= 4) {
lseek(fd, 0, SEEK_SET);
#endif // WIN32
ret = read(fd, magic, 4);
if (ret == 4 && strncmp(magic, MM_IDX_MAGIC, 4) == 0)
is_idx = 1;
+116
View File
@@ -0,0 +1,116 @@
#ifndef KETOPT_H
#define KETOPT_H
#include <string.h> /* for strchr() and strncmp() */
#define ko_no_argument 0
#define ko_required_argument 1
#define ko_optional_argument 2
typedef struct {
int ind; /* equivalent to optind */
int opt; /* equivalent to optopt */
char *arg; /* equivalent to optarg */
int longidx; /* index of a long option; or -1 if short */
/* private variables not intended for external uses */
int i, pos, n_args;
} ketopt_t;
typedef struct {
char *name;
int has_arg;
int val;
} ko_longopt_t;
static ketopt_t KETOPT_INIT = { 1, 0, 0, -1, 1, 0, 0 };
static void ketopt_permute(char *argv[], int j, int n) /* move argv[j] over n elements to the left */
{
int k;
char *p = argv[j];
for (k = 0; k < n; ++k)
argv[j - k] = argv[j - k - 1];
argv[j - k] = p;
}
/**
* Parse command-line options and arguments
*
* This fuction has a similar interface to GNU's getopt_long(). Each call
* parses one option and returns the option name. s->arg points to the option
* argument if present. The function returns -1 when all command-line arguments
* are parsed. In this case, s->ind is the index of the first non-option
* argument.
*
* @param s status; shall be initialized to KETOPT_INIT on the first call
* @param argc length of argv[]
* @param argv list of command-line arguments; argv[0] is ignored
* @param permute non-zero to move options ahead of non-option arguments
* @param ostr option string
* @param longopts long options
*
* @return ASCII for a short option; ko_longopt_t::val for a long option; -1 if
* argv[] is fully processed; '?' for an unknown option or an ambiguous
* long option; ':' if an option argument is missing
*/
static int ketopt(ketopt_t *s, int argc, char *argv[], int permute, const char *ostr, const ko_longopt_t *longopts)
{
int opt = -1, i0, j;
if (permute) {
while (s->i < argc && (argv[s->i][0] != '-' || argv[s->i][1] == '\0'))
++s->i, ++s->n_args;
}
s->arg = 0, s->longidx = -1, i0 = s->i;
if (s->i >= argc || argv[s->i][0] != '-' || argv[s->i][1] == '\0') {
s->ind = s->i - s->n_args;
return -1;
}
if (argv[s->i][0] == '-' && argv[s->i][1] == '-') { /* "--" or a long option */
if (argv[s->i][2] == '\0') { /* a bare "--" */
ketopt_permute(argv, s->i, s->n_args);
++s->i, s->ind = s->i - s->n_args;
return -1;
}
s->opt = 0, opt = '?', s->pos = -1;
if (longopts) { /* parse long options */
int k, n_matches = 0;
const ko_longopt_t *o = 0;
for (j = 2; argv[s->i][j] != '\0' && argv[s->i][j] != '='; ++j) {} /* find the end of the option name */
for (k = 0; longopts[k].name != 0; ++k)
if (strncmp(&argv[s->i][2], longopts[k].name, j - 2) == 0)
++n_matches, o = &longopts[k];
if (n_matches == 1) {
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') {
if (s->i < argc - 1) s->arg = argv[++s->i];
else opt = ':'; /* missing option argument */
}
}
}
} else { /* a short option */
char *p;
if (s->pos == 0) s->pos = 1;
opt = s->opt = argv[s->i][s->pos++];
p = strchr(ostr, opt);
if (p == 0) {
opt = '?'; /* unknown option */
} else if (p[1] == ':') {
if (argv[s->i][s->pos] == 0) {
if (s->i < argc - 1) s->arg = argv[++s->i];
else opt = ':'; /* missing option argument */
} else s->arg = &argv[s->i][s->pos];
s->pos = -1;
}
}
if (s->pos < 0 || argv[s->i][s->pos] == 0) {
++s->i, s->pos = 0;
if (s->n_args > 0) /* permute */
for (j = i0; j < s->i; ++j)
ketopt_permute(argv, j, s->n_args);
}
s->ind = s->i - s->n_args;
return opt;
}
#endif
+14 -15
View File
@@ -17,18 +17,20 @@
void __cpuidex(int cpuid[4], int func_id, int subfunc_id)
{
#if defined(__x86_64__)
asm volatile ("cpuid"
__asm__ volatile ("cpuid"
: "=a" (cpuid[0]), "=b" (cpuid[1]), "=c" (cpuid[2]), "=d" (cpuid[3])
: "0" (func_id), "2" (subfunc_id));
#else // on 32bit, ebx can NOT be used as PIC code
asm volatile ("xchgl %%ebx, %1; cpuid; xchgl %%ebx, %1"
__asm__ volatile ("xchgl %%ebx, %1; cpuid; xchgl %%ebx, %1"
: "=a" (cpuid[0]), "=r" (cpuid[1]), "=c" (cpuid[2]), "=d" (cpuid[3])
: "0" (func_id), "2" (subfunc_id));
#endif
}
#endif
int x86_simd(void)
static int ksw_simd = -1;
static int x86_simd(void)
{
int flag = 0, cpuid[4], max_id;
__cpuidex(cpuid, 0, 0);
@@ -54,11 +56,10 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
{
extern void ksw_extz2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
extern void ksw_extz2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t q, int8_t e, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
unsigned simd;
simd = x86_simd();
if (simd & SIMD_SSE4_1)
if (ksw_simd < 0) ksw_simd = x86_simd();
if (ksw_simd & SIMD_SSE4_1)
ksw_extz2_sse41(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, end_bonus, flag, ez);
else if (simd & SIMD_SSE2)
else if (ksw_simd & SIMD_SSE2)
ksw_extz2_sse2(km, qlen, query, tlen, target, m, mat, q, e, w, zdrop, end_bonus, flag, ez);
else abort();
}
@@ -70,11 +71,10 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
extern void ksw_extd2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez);
unsigned simd;
simd = x86_simd();
if (simd & SIMD_SSE4_1)
if (ksw_simd < 0) ksw_simd = x86_simd();
if (ksw_simd & SIMD_SSE4_1)
ksw_extd2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, ez);
else if (simd & SIMD_SSE2)
else if (ksw_simd & SIMD_SSE2)
ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, ez);
else abort();
}
@@ -86,11 +86,10 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
extern void ksw_exts2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat,
int8_t q, int8_t e, int8_t q2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez);
unsigned simd;
simd = x86_simd();
if (simd & SIMD_SSE4_1)
if (ksw_simd < 0) ksw_simd = x86_simd();
if (ksw_simd & SIMD_SSE4_1)
ksw_exts2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
else if (simd & SIMD_SSE2)
else if (ksw_simd & SIMD_SSE2)
ksw_exts2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, noncan, zdrop, flag, ez);
else abort();
}
+1 -1
View File
@@ -100,7 +100,7 @@ void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin
for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF;
}
if (with_cigar) {
mem2 = (uint8_t*)kmalloc(km, ((qlen + tlen - 1) * n_col_ + 1) * 16);
mem2 = (uint8_t*)kmalloc(km, ((size_t)(qlen + tlen - 1) * n_col_ + 1) * 16);
p = (__m128i*)(((size_t)mem2 + 15) >> 4 << 4);
off = (int*)kmalloc(km, (qlen + tlen - 1) * sizeof(int) * 2);
off_end = off + qlen + tlen - 1;
+152 -139
View File
@@ -4,13 +4,9 @@
#include "bseq.h"
#include "minimap.h"
#include "mmpriv.h"
#ifdef HAVE_GETOPT
#include <getopt.h>
#else
#include "getopt.h"
#endif
#include "ketopt.h"
#define MM_VERSION "2.11-r815-dirty"
#define MM_VERSION "2.14-r883"
#ifdef __linux__
#include <sys/resource.h>
@@ -26,51 +22,54 @@ void liftrlimit()
void liftrlimit() {}
#endif
static struct option long_options[] = {
{ "bucket-bits", required_argument, 0, 0 },
{ "mb-size", required_argument, 0, 'K' },
{ "seed", required_argument, 0, 0 },
{ "no-kalloc", no_argument, 0, 0 },
{ "print-qname", no_argument, 0, 0 },
{ "no-self", no_argument, 0, 'D' },
{ "print-seeds", no_argument, 0, 0 },
{ "max-chain-skip", required_argument, 0, 0 },
{ "min-dp-len", required_argument, 0, 0 },
{ "print-aln-seq", no_argument, 0, 0 },
{ "splice", no_argument, 0, 0 },
{ "cost-non-gt-ag", required_argument, 0, 'C' },
{ "no-long-join", no_argument, 0, 0 },
{ "sr", no_argument, 0, 0 },
{ "frag", required_argument, 0, 0 },
{ "secondary", required_argument, 0, 0 },
{ "cs", optional_argument, 0, 0 },
{ "end-bonus", required_argument, 0, 0 },
{ "no-pairing", no_argument, 0, 0 },
{ "splice-flank", required_argument, 0, 0 },
{ "idx-no-seq", no_argument, 0, 0 },
{ "end-seed-pen", required_argument, 0, 0 }, // 21
{ "for-only", no_argument, 0, 0 }, // 22
{ "rev-only", no_argument, 0, 0 }, // 23
{ "heap-sort", required_argument, 0, 0 }, // 24
{ "all-chain", no_argument, 0, 'P' },
{ "dual", required_argument, 0, 0 }, // 26
{ "max-clip-ratio", required_argument, 0, 0 }, // 27
{ "min-occ-floor", required_argument, 0, 0 }, // 28
{ "MD", no_argument, 0, 0 }, // 29
{ "lj-min-ratio", required_argument, 0, 0 }, // 30
{ "score-N", required_argument, 0, 0 }, // 31
{ "eqx", no_argument, 0, 0 }, // 32
{ "paf-no-hit", no_argument, 0, 0 }, // 33
{ "split-prefix", required_argument, 0, 0 }, // 34
{ "help", no_argument, 0, 'h' },
{ "max-intron-len", required_argument, 0, 'G' },
{ "version", no_argument, 0, 'V' },
{ "min-count", required_argument, 0, 'n' },
{ "min-chain-score",required_argument, 0, 'm' },
{ "mask-level", required_argument, 0, 'M' },
{ "min-dp-score", required_argument, 0, 's' },
{ "sam", no_argument, 0, 'a' },
{ 0, 0, 0, 0}
static ko_longopt_t long_options[] = {
{ "bucket-bits", ko_required_argument, 300 },
{ "mb-size", ko_required_argument, 'K' },
{ "seed", ko_required_argument, 302 },
{ "no-kalloc", ko_no_argument, 303 },
{ "print-qname", ko_no_argument, 304 },
{ "no-self", ko_no_argument, 'D' },
{ "print-seeds", ko_no_argument, 306 },
{ "max-chain-skip", ko_required_argument, 307 },
{ "min-dp-len", ko_required_argument, 308 },
{ "print-aln-seq", ko_no_argument, 309 },
{ "splice", ko_no_argument, 310 },
{ "cost-non-gt-ag", ko_required_argument, 'C' },
{ "no-long-join", ko_no_argument, 312 },
{ "sr", ko_no_argument, 313 },
{ "frag", ko_required_argument, 314 },
{ "secondary", ko_required_argument, 315 },
{ "cs", ko_optional_argument, 316 },
{ "end-bonus", ko_required_argument, 317 },
{ "no-pairing", ko_no_argument, 318 },
{ "splice-flank", ko_required_argument, 319 },
{ "idx-no-seq", ko_no_argument, 320 },
{ "end-seed-pen", ko_required_argument, 321 },
{ "for-only", ko_no_argument, 322 },
{ "rev-only", ko_no_argument, 323 },
{ "heap-sort", ko_required_argument, 324 },
{ "all-chain", ko_no_argument, 'P' },
{ "dual", ko_required_argument, 326 },
{ "max-clip-ratio", ko_required_argument, 327 },
{ "min-occ-floor", ko_required_argument, 328 },
{ "MD", ko_no_argument, 329 },
{ "lj-min-ratio", ko_required_argument, 330 },
{ "score-N", ko_required_argument, 331 },
{ "eqx", ko_no_argument, 332 },
{ "paf-no-hit", ko_no_argument, 333 },
{ "split-prefix", ko_required_argument, 334 },
{ "no-end-flt", ko_no_argument, 335 },
{ "hard-mask-level",ko_no_argument, 336 },
{ "cap-sw-mem", ko_required_argument, 337 },
{ "help", ko_no_argument, 'h' },
{ "max-intron-len", ko_required_argument, 'G' },
{ "version", ko_no_argument, 'V' },
{ "min-count", ko_required_argument, 'n' },
{ "min-chain-score",ko_required_argument, 'm' },
{ "mask-level", ko_required_argument, 'M' },
{ "min-dp-score", ko_required_argument, 's' },
{ "sam", ko_no_argument, 'a' },
{ 0, 0, 0 }
};
static inline int64_t mm_parse_num(const char *str)
@@ -99,10 +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, long_idx;
int i, c, n_threads = 3, n_parts, old_best_n = -1;
char *fnw = 0, *rg = 0, *s;
FILE *fp_help = stderr;
mm_idx_reader_t *idx_rdr;
@@ -113,30 +113,36 @@ int main(int argc, char *argv[])
mm_realtime0 = realtime();
mm_set_opt(0, &ipt, &opt);
while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) // apply option -x/preset first
while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) { // test command line options and apply option -x/preset first
if (c == 'x') {
if (mm_set_opt(optarg, &ipt, &opt) < 0) {
fprintf(stderr, "[ERROR] unknown preset '%s'\n", optarg);
if (mm_set_opt(o.arg, &ipt, &opt) < 0) {
fprintf(stderr, "[ERROR] unknown preset '%s'\n", o.arg);
return 1;
}
break;
} else if (c == ':') {
fprintf(stderr, "[ERROR] missing option argument\n");
return 1;
} else if (c == '?') {
fprintf(stderr, "[ERROR] unknown option in \"%s\"\n", argv[o.i - 1]);
return 1;
}
optind = 0; // for musl getopt, optind=0 has the same effect as optreset=1; older libc doesn't have optreset
}
o = KETOPT_INIT;
while ((c = getopt_long(argc, argv, opt_str, long_options, &long_idx)) >= 0) {
if (c == 'w') ipt.w = atoi(optarg);
else if (c == 'k') ipt.k = atoi(optarg);
while ((c = ketopt(&o, argc, argv, 1, opt_str, long_options)) >= 0) {
if (c == 'w') ipt.w = atoi(o.arg);
else if (c == 'k') ipt.k = atoi(o.arg);
else if (c == 'H') ipt.flag |= MM_I_HPC;
else if (c == 'd') fnw = optarg; // the above are indexing related options, except -I
else if (c == 'r') opt.bw = (int)mm_parse_num(optarg);
else if (c == 't') n_threads = atoi(optarg);
else if (c == 'v') mm_verbose = atoi(optarg);
else if (c == 'g') opt.max_gap = (int)mm_parse_num(optarg);
else if (c == 'G') mm_mapopt_max_intron_len(&opt, (int)mm_parse_num(optarg));
else if (c == 'F') opt.max_frag_len = (int)mm_parse_num(optarg);
else if (c == 'N') opt.best_n = atoi(optarg);
else if (c == 'p') opt.pri_ratio = atof(optarg);
else if (c == 'M') opt.mask_level = atof(optarg);
else if (c == 'd') fnw = o.arg; // the above are indexing related options, except -I
else if (c == 'r') opt.bw = (int)mm_parse_num(o.arg);
else if (c == 't') n_threads = atoi(o.arg);
else if (c == 'v') mm_verbose = atoi(o.arg);
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') 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;
else if (c == 'D') opt.flag |= MM_F_NO_DIAG;
else if (c == 'P') opt.flag |= MM_F_ALL_CHAINS;
@@ -146,64 +152,67 @@ int main(int argc, char *argv[])
else if (c == 'Y') opt.flag |= MM_F_SOFTCLIP;
else if (c == 'L') opt.flag |= MM_F_LONG_CIGAR;
else if (c == 'y') opt.flag |= MM_F_COPY_COMMENT;
else if (c == 'T') opt.sdust_thres = atoi(optarg);
else if (c == 'n') opt.min_cnt = atoi(optarg);
else if (c == 'm') opt.min_chain_score = atoi(optarg);
else if (c == 'A') opt.a = atoi(optarg);
else if (c == 'B') opt.b = atoi(optarg);
else if (c == 's') opt.min_dp_max = atoi(optarg);
else if (c == 'C') opt.noncan = atoi(optarg);
else if (c == 'I') ipt.batch_size = mm_parse_num(optarg);
else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(optarg);
else if (c == 'R') rg = optarg;
else if (c == 'T') opt.sdust_thres = atoi(o.arg);
else if (c == 'n') opt.min_cnt = atoi(o.arg);
else if (c == 'm') opt.min_chain_score = atoi(o.arg);
else if (c == 'A') opt.a = atoi(o.arg);
else if (c == 'B') opt.b = atoi(o.arg);
else if (c == 's') opt.min_dp_max = atoi(o.arg);
else if (c == 'C') opt.noncan = atoi(o.arg);
else if (c == 'I') ipt.batch_size = mm_parse_num(o.arg);
else if (c == 'K') opt.mini_batch_size = (int)mm_parse_num(o.arg);
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 == 0 && long_idx == 0) ipt.bucket_bits = atoi(optarg); // --bucket-bits
else if (c == 0 && long_idx == 2) opt.seed = atoi(optarg); // --seed
else if (c == 0 && long_idx == 3) mm_dbg_flag |= MM_DBG_NO_KALLOC; // --no-kalloc
else if (c == 0 && long_idx == 4) mm_dbg_flag |= MM_DBG_PRINT_QNAME; // --print-qname
else if (c == 0 && long_idx == 6) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_SEED, n_threads = 1; // --print-seed
else if (c == 0 && long_idx == 7) opt.max_chain_skip = atoi(optarg); // --max-chain-skip
else if (c == 0 && long_idx == 8) opt.min_ksw_len = atoi(optarg); // --min-dp-len
else if (c == 0 && long_idx == 9) mm_dbg_flag |= MM_DBG_PRINT_QNAME | MM_DBG_PRINT_ALN_SEQ, n_threads = 1; // --print-aln-seq
else if (c == 0 && long_idx ==10) opt.flag |= MM_F_SPLICE; // --splice
else if (c == 0 && long_idx ==12) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
else if (c == 0 && long_idx ==13) opt.flag |= MM_F_SR; // --sr
else if (c == 0 && long_idx ==17) opt.end_bonus = atoi(optarg); // --end-bonus
else if (c == 0 && long_idx ==18) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
else if (c == 0 && long_idx ==20) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
else if (c == 0 && long_idx ==21) opt.anchor_ext_shift = atoi(optarg); // --end-seed-pen
else if (c == 0 && long_idx ==22) opt.flag |= MM_F_FOR_ONLY; // --for-only
else if (c == 0 && long_idx ==23) opt.flag |= MM_F_REV_ONLY; // --rev-only
else if (c == 0 && long_idx ==27) opt.max_clip_ratio = atof(optarg); // --max-clip-ratio
else if (c == 0 && long_idx ==28) opt.min_mid_occ = atoi(optarg); // --min-occ-floor
else if (c == 0 && long_idx ==29) opt.flag |= MM_F_OUT_MD; // --MD
else if (c == 0 && long_idx ==30) opt.min_join_flank_ratio = atof(optarg); // --lj-min-ratio
else if (c == 0 && long_idx ==31) opt.sc_ambi = atoi(optarg); // --score-N
else if (c == 0 && long_idx ==32) opt.flag |= MM_F_EQX; // --eqx
else if (c == 0 && long_idx ==33) opt.flag |= MM_F_PAF_NO_HIT; // --paf-no-hit
else if (c == 0 && long_idx ==34) opt.split_prefix = optarg; // --split-prefix
else if (c == 0 && long_idx == 14) { // --frag
yes_or_no(&opt, MM_F_FRAG_MODE, long_idx, optarg, 1);
} else if (c == 0 && long_idx == 15) { // --secondary
yes_or_no(&opt, MM_F_NO_PRINT_2ND, long_idx, optarg, 0);
} else if (c == 0 && long_idx == 16) { // --cs
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 == 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
else if (c == 312) opt.flag |= MM_F_NO_LJOIN; // --no-long-join
else if (c == 313) opt.flag |= MM_F_SR; // --sr
else if (c == 317) opt.end_bonus = atoi(o.arg); // --end-bonus
else if (c == 318) opt.flag |= MM_F_INDEPEND_SEG; // --no-pairing
else if (c == 320) ipt.flag |= MM_I_NO_SEQ; // --idx-no-seq
else if (c == 321) opt.anchor_ext_shift = atoi(o.arg); // --end-seed-pen
else if (c == 322) opt.flag |= MM_F_FOR_ONLY; // --for-only
else if (c == 323) opt.flag |= MM_F_REV_ONLY; // --rev-only
else if (c == 327) opt.max_clip_ratio = atof(o.arg); // --max-clip-ratio
else if (c == 328) opt.min_mid_occ = atoi(o.arg); // --min-occ-floor
else if (c == 329) opt.flag |= MM_F_OUT_MD; // --MD
else if (c == 330) opt.min_join_flank_ratio = atof(o.arg); // --lj-min-ratio
else if (c == 331) opt.sc_ambi = atoi(o.arg); // --score-N
else if (c == 332) opt.flag |= MM_F_EQX; // --eqx
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) opt.flag |= MM_F_HARD_MLEVEL; // --hard-mask-level
else if (c == 337) opt.max_sw_mat = mm_parse_num(o.arg); // --cap-sw-mat
else if (c == 314) { // --frag
yes_or_no(&opt, MM_F_FRAG_MODE, o.longidx, o.arg, 1);
} else if (c == 315) { // --secondary
yes_or_no(&opt, MM_F_NO_PRINT_2ND, o.longidx, o.arg, 0);
} else if (c == 316) { // --cs
opt.flag |= MM_F_OUT_CS | MM_F_CIGAR;
if (optarg == 0 || strcmp(optarg, "short") == 0) {
if (o.arg == 0 || strcmp(o.arg, "short") == 0) {
opt.flag &= ~MM_F_OUT_CS_LONG;
} else if (strcmp(optarg, "long") == 0) {
} else if (strcmp(o.arg, "long") == 0) {
opt.flag |= MM_F_OUT_CS_LONG;
} else if (strcmp(optarg, "none") == 0) {
} else if (strcmp(o.arg, "none") == 0) {
opt.flag &= ~MM_F_OUT_CS;
} else if (mm_verbose >= 2) {
fprintf(stderr, "[WARNING]\033[1;31m --cs only takes 'short' or 'long'. Invalid values are assumed to be 'short'.\033[0m\n");
}
} else if (c == 0 && long_idx == 19) { // --splice-flank
yes_or_no(&opt, MM_F_SPLICE_FLANK, long_idx, optarg, 1);
} else if (c == 0 && long_idx == 24) { // --heap-sort
yes_or_no(&opt, MM_F_HEAP_SORT, long_idx, optarg, 1);
} else if (c == 0 && long_idx == 26) { // --dual
yes_or_no(&opt, MM_F_NO_DUAL, long_idx, optarg, 0);
} else if (c == 319) { // --splice-flank
yes_or_no(&opt, MM_F_SPLICE_FLANK, o.longidx, o.arg, 1);
} else if (c == 324) { // --heap-sort
yes_or_no(&opt, MM_F_HEAP_SORT, o.longidx, o.arg, 1);
} else if (c == 326) { // --dual
yes_or_no(&opt, MM_F_NO_DUAL, o.longidx, o.arg, 0);
} else if (c == 'S') {
opt.flag |= MM_F_OUT_CS | MM_F_CIGAR | MM_F_OUT_CS_LONG;
if (mm_verbose >= 2)
@@ -214,27 +223,27 @@ int main(int argc, char *argv[])
} else if (c == 'f') {
double x;
char *p;
x = strtod(optarg, &p);
x = strtod(o.arg, &p);
if (x < 1.0) opt.mid_occ_frac = x, opt.mid_occ = 0;
else opt.mid_occ = (int)(x + .499);
if (*p == ',') opt.max_occ = (int)(strtod(p+1, &p) + .499);
} else if (c == 'u') {
if (*optarg == 'b') opt.flag |= MM_F_SPLICE_FOR|MM_F_SPLICE_REV; // both strands
else if (*optarg == 'f') opt.flag |= MM_F_SPLICE_FOR, opt.flag &= ~MM_F_SPLICE_REV; // match GT-AG
else if (*optarg == 'r') opt.flag |= MM_F_SPLICE_REV, opt.flag &= ~MM_F_SPLICE_FOR; // match CT-AC (reverse complement of GT-AG)
else if (*optarg == 'n') opt.flag &= ~(MM_F_SPLICE_FOR|MM_F_SPLICE_REV); // don't try to match the GT-AG signal
if (*o.arg == 'b') opt.flag |= MM_F_SPLICE_FOR|MM_F_SPLICE_REV; // both strands
else if (*o.arg == 'f') opt.flag |= MM_F_SPLICE_FOR, opt.flag &= ~MM_F_SPLICE_REV; // match GT-AG
else if (*o.arg == 'r') opt.flag |= MM_F_SPLICE_REV, opt.flag &= ~MM_F_SPLICE_FOR; // match CT-AC (reverse complement of GT-AG)
else if (*o.arg == 'n') opt.flag &= ~(MM_F_SPLICE_FOR|MM_F_SPLICE_REV); // don't try to match the GT-AG signal
else {
fprintf(stderr, "[ERROR]\033[1;31m unrecognized cDNA direction\033[0m\n");
return 1;
}
} else if (c == 'z') {
opt.zdrop = opt.zdrop_inv = strtol(optarg, &s, 10);
opt.zdrop = opt.zdrop_inv = strtol(o.arg, &s, 10);
if (*s == ',') opt.zdrop_inv = strtol(s + 1, &s, 10);
} else if (c == 'O') {
opt.q = opt.q2 = strtol(optarg, &s, 10);
opt.q = opt.q2 = strtol(o.arg, &s, 10);
if (*s == ',') opt.q2 = strtol(s + 1, &s, 10);
} else if (c == 'E') {
opt.e = opt.e2 = strtol(optarg, &s, 10);
opt.e = opt.e2 = strtol(o.arg, &s, 10);
if (*s == ',') opt.e2 = strtol(s + 1, &s, 10);
}
}
@@ -246,8 +255,12 @@ 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 == optind || fp_help == stdout) {
if (argc == o.ind || fp_help == stdout) {
fprintf(fp_help, "Usage: minimap2 [options] <target.fa>|<target.idx> [query.fa] [...]\n");
fprintf(fp_help, "Options:\n");
fprintf(fp_help, " Indexing:\n");
@@ -297,20 +310,20 @@ int main(int argc, char *argv[])
fprintf(fp_help, " - asm5/asm10/asm20: asm-to-ref mapping, for ~0.1/1/5%% sequence divergence\n");
fprintf(fp_help, " - splice: long-read spliced alignment\n");
fprintf(fp_help, " - sr: genomic short-read mapping\n");
fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of command-line options.\n");
fprintf(fp_help, "\nSee `man ./minimap2.1' for detailed description of these and other advanced command-line options.\n");
return fp_help == stdout? 0 : 1;
}
if ((opt.flag & MM_F_SR) && argc - optind > 3) {
if ((opt.flag & MM_F_SR) && argc - o.ind > 3) {
fprintf(stderr, "[ERROR] incorrect input: in the sr mode, please specify no more than two query files.\n");
return 1;
}
idx_rdr = mm_idx_reader_open(argv[optind], &ipt, fnw);
idx_rdr = mm_idx_reader_open(argv[o.ind], &ipt, fnw);
if (idx_rdr == 0) {
fprintf(stderr, "[ERROR] failed to open file '%s'\n", argv[optind]);
fprintf(stderr, "[ERROR] failed to open file '%s'\n", argv[o.ind]);
return 1;
}
if (!idx_rdr->is_idx && fnw == 0 && argc - optind < 2) {
if (!idx_rdr->is_idx && fnw == 0 && argc - o.ind < 2) {
fprintf(stderr, "[ERROR] missing input: please specify a query file to map or option -d to keep the index\n");
mm_idx_reader_close(idx_rdr);
return 1;
@@ -336,13 +349,13 @@ int main(int argc, char *argv[])
if (mm_verbose >= 3)
fprintf(stderr, "[M::%s::%.3f*%.2f] loaded/built the index for %d target sequence(s)\n",
__func__, realtime() - mm_realtime0, cputime() / (realtime() - mm_realtime0), mi->n_seq);
if (argc != optind + 1) mm_mapopt_update(&opt, mi);
if (argc != o.ind + 1) mm_mapopt_update(&opt, mi);
if (mm_verbose >= 3) mm_idx_stat(mi);
if (!(opt.flag & MM_F_FRAG_MODE)) {
for (i = optind + 1; i < argc; ++i)
for (i = o.ind + 1; i < argc; ++i)
mm_map_file(mi, argv[i], &opt, n_threads);
} else {
mm_map_file_frag(mi, argc - (optind + 1), (const char**)&argv[optind + 1], &opt, n_threads);
mm_map_file_frag(mi, argc - (o.ind + 1), (const char**)&argv[o.ind + 1], &opt, n_threads);
}
mm_idx_destroy(mi);
}
@@ -350,7 +363,7 @@ int main(int argc, char *argv[])
mm_idx_reader_close(idx_rdr);
if (opt.split_prefix)
mm_split_merge(argc - (optind + 1), (const char**)&argv[optind + 1], &opt, n_parts);
mm_split_merge(argc - (o.ind + 1), (const char**)&argv[o.ind + 1], &opt, n_parts);
if (fflush(stdout) == EOF) {
fprintf(stderr, "[ERROR] failed to write the results\n");
@@ -362,7 +375,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;
}
+11 -5
View File
@@ -29,6 +29,11 @@ void mm_tbuf_destroy(mm_tbuf_t *b)
free(b);
}
void *mm_tbuf_get_km(mm_tbuf_t *b)
{
return b->km;
}
static int mm_dust_minier(void *km, int n, mm128_t *a, int l_seq, const char *seq, int sdust_thres)
{
int n_dreg, j, k, u = 0;
@@ -243,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
@@ -256,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);
}
@@ -354,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);
}
@@ -497,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]);
@@ -596,6 +601,7 @@ static void *worker_pipeline(void *shared, int step, void *in)
free(s->reg[i]);
free(s->seq[i].seq); free(s->seq[i].name);
if (s->seq[i].qual) free(s->seq[i].qual);
if (s->seq[i].comment) free(s->seq[i].comment);
}
}
free(s->reg); free(s->n_reg); free(s->seq); // seg_off, n_seg, rep_len and frag_gap were allocated with reg; no memory leak here
@@ -682,7 +688,7 @@ int mm_split_merge(int n_segs, const char **fn, const mm_mapopt_t *opt, int n_sp
for (pl.rid_shift[0] = 0, i = 1; i < n_split_idx; ++i)
pl.rid_shift[i] += pl.rid_shift[i - 1];
if (opt->flag & MM_F_OUT_SAM)
for (i = 0; i < pl.mi->n_seq; ++i)
for (i = 0; i < (int32_t)pl.mi->n_seq; ++i)
printf("@SQ\tSN:%s\tLN:%d\n", pl.mi->seq[i].name, pl.mi->seq[i].len);
kt_pipeline(2, worker_pipeline, &pl, 3);
+21
View File
@@ -33,6 +33,8 @@
#define MM_F_COPY_COMMENT 0x2000000
#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
@@ -137,6 +139,7 @@ typedef struct {
int32_t mid_occ; // ignore seeds with occurrences above this threshold
int32_t max_occ;
int mini_batch_size; // size of a batch of query bases to process in parallel
int64_t max_sw_mat;
const char *split_prefix;
} mm_mapopt_t;
@@ -301,6 +304,8 @@ mm_tbuf_t *mm_tbuf_init(void);
*/
void mm_tbuf_destroy(mm_tbuf_t *b);
void *mm_tbuf_get_km(mm_tbuf_t *b);
/**
* Align a query sequence against an index
*
@@ -337,6 +342,22 @@ int mm_map_file(const mm_idx_t *idx, const char *fn, const mm_mapopt_t *opt, int
int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_mapopt_t *opt, int n_threads);
/**
* Generate the cs tag (new in 2.12)
*
* @param km memory blocks; set to NULL if unsure
* @param buf buffer to write the cs/MD tag; typicall NULL on the first call
* @param max_len max length of the buffer; typically set to 0 on the first call
* @param mi index
* @param r alignment
* @param seq query sequence
* @param no_iden true to use : instead of =
*
* @return the length of cs
*/
int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden);
int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq);
// query sequence name and sequence in the minimap2 index
int mm_idx_index_name(mm_idx_t *mi);
int mm_idx_name2id(const mm_idx_t *mi, const char *name);
+21 -3
View File
@@ -1,4 +1,4 @@
.TH minimap2 1 "20 June 2018" "minimap2-2.11 (r797)" "Bioinformatics tools"
.TH minimap2 1 "5 November 2018" "minimap2-2.14 (r883)" "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
@@ -269,6 +274,10 @@ Only map to the reverse complement strand of the reference sequences.
.BR --heap-sort = no | yes
If yes, sort anchors with heap merge, instead of radix sort. Heap merge is
faster for short reads, but slower for long reads. [no]
.TP
.B --no-pairing
Treat two reads in a pair as independent reads. The mate related fields in SAM
are still properly populated.
.SS Alignment options
.TP 10
.BI -A \ INT
@@ -360,6 +369,15 @@ the length of the terminal gap in the chain. This option is only effective
with
.BR --splice .
It helps to avoid tiny terminal exons. [6]
.TP
.B --no-end-flt
Don't filter seeds towards the ends of chains before performing base-level
alignment.
.TP
.BI --cap-sw-mem \ NUM
Skip alignment if the DP matrix size is above
.IR NUM .
Set 0 to disable [0].
.SS Input/output options
.TP 10
.B -a
@@ -484,7 +502,7 @@ Up to 10% sequence divergence.
.B asm20
Long assembly to reference mapping
.RB ( -k19
.B -w10 -A1 -B6 -O6,26 -E2,1 -s200 -z200
.B -w10 -A1 -B4 -O6,26 -E2,1 -s200 -z200
.BR --min-occ-floor=100 ).
Up to 20% sequence divergence.
.TP
+1 -2
View File
@@ -116,8 +116,7 @@ long peakrss(void)
double realtime(void)
{
struct timeval tp;
struct timezone tzp;
gettimeofday(&tp, &tzp);
gettimeofday(&tp, NULL);
return tp.tv_sec + tp.tv_usec * 1e-6;
}
+173 -15
View File
@@ -1,6 +1,6 @@
#!/usr/bin/env k8
var paftools_version = 'r767';
var paftools_version = '2.14-r883';
/*****************************
***** 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=<ID=QSTART,Number=1,Type=Integer,Description="Query start">');
print('##INFO=<ID=QSTRAND,Number=1,Type=String,Description="Query strand">');
print('##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">');
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 = [];
@@ -562,10 +564,12 @@ function paf_call(args)
function paf_asmstat(args)
{
var c, min_seg_len = 10000, max_diff = 0.01;
while ((c = getopt(args, "l:d:")) != null) {
var c, min_seg_len = 10000, max_diff = 0.01, bp_flank_len = 0, bp_gap_len = 0;
while ((c = getopt(args, "l:d:b:g:")) != null) {
if (c == 'l') min_seg_len = parseInt(getopt.arg);
else if (c == 'd') max_diff = parseFloat(getopt.arg);
else if (c == 'b') bp_flank_len = parseInt(getopt.arg);
else if (c == 'g') bp_gap_len = parseInt(getopt.arg);
}
if (getopt.ind == args.length) {
print("Usage: paftools.js asmstat [options] <ref.fa.fai> <asm1.paf> [...]");
@@ -585,7 +589,7 @@ function paf_asmstat(args)
}
file.close();
function process_query(qblocks, qblock_len, bp) {
function process_query(qblocks, qblock_len, bp, qi) {
qblocks.sort(function(a,b) { return a[0]-b[0]; });
var last_k = null, last_blen = null, st = -1, en = -1, qcov = 0;
for (var k = 0; k < qblocks.length; ++k) {
@@ -610,6 +614,7 @@ function paf_asmstat(args)
var min = blen < last_blen? blen : last_blen;
var flank = k == 0? min : blen;
bp.push([flank, gap]);
qi.bp.push([flank, gap]);
}
last_k = k, last_blen = blen;
}
@@ -662,16 +667,22 @@ function paf_asmstat(args)
for (var i = 0; i < n_asm; ++i) {
var n_breaks = 0, qcov = 0;
var fn = args[getopt.ind + 1 + i];
header.push(fn.replace(/.paf(.gz)?$/, ""));
var label = fn.replace(/.paf(.gz)?$/, "");
header.push(label);
var ref_blocks = [], qblock_len = [], qblocks = [], bp = [];
var query = {};
var query = {}, qinfo = {};
var last_qname = null;
file = new File(fn);
while (file.readline(buf) >= 0) {
var m, line = buf.toString();
var t = line.split("\t");
t[1] = parseInt(t[1]);
if (t.length >= 2) query[t[0]] = t[1];
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 (!/\ttp:A:[PI]/.test(line)) continue;
if ((m = /\tcg:Z:(\S+)/.exec(line)) == null) continue;
@@ -688,7 +699,7 @@ function paf_asmstat(args)
if (t[3] - t[2] < min_seg_len) continue;
if (t[0] != last_qname) {
if (last_qname != null)
qcov += process_query(qblocks, qblock_len, bp);
qcov += process_query(qblocks, qblock_len, bp, qinfo[last_qname]);
qblocks = [];
last_qname = t[0];
}
@@ -696,7 +707,7 @@ function paf_asmstat(args)
qblocks.push([t[2], t[3], t[4], t[5], t[7], t[8]]);
}
if (last_qname != null)
qcov += process_query(qblocks, qblock_len, bp);
qcov += process_query(qblocks, qblock_len, bp, qinfo[last_qname]);
file.close();
// compute NG50
@@ -731,12 +742,157 @@ function paf_asmstat(args)
rst[5][i] = n_breaks;
rst[6][i] = count_bp(bp, 500, 0);
rst[7][i] = count_bp(bp, 500, 10000);
// nb-plot
var qa = [];
for (var qn in qinfo)
qa.push([qinfo[qn].len, qinfo[qn].bp]);
qa = qa.sort(function(a, b) { return b[0] - a[0] });
var sum = 0, n_bp = 0, next_quantile = 0.1;
for (var j = 0; j < qa.length; ++j) {
sum += qa[j][0];
for (var k = 0; k < qa[j][1].length; ++k)
if (qa[j][1][k][0] >= bp_flank_len && qa[j][1][k][1] >= bp_gap_len)
++n_bp;
if (sum >= ref_len * next_quantile) {
print(label, Math.floor(next_quantile * 100 + .5), qa[j][0], (sum / n_bp).toFixed(0), n_bp);
next_quantile += 0.1;
if (next_quantile >= 1.0) break;
}
}
}
buf.destroy();
if (bp_flank_len <= 0) {
print(header.join("\t"));
for (var i = 0; i < labels.length; ++i)
print(labels[i], rst[i].join("\t"));
}
}
function paf_asmgene(args)
{
var c, opt = { min_cov:0.99, min_iden:0.99 }, print_err = false;
while ((c = getopt(args, "i:c:e")) != null)
if (c == 'i') opt.min_iden = parseFloat(getopt.arg);
else if (c == 'c') opt.min_cov = parseFloat(getopt.arg);
else if (c == 'e') print_err = true;
var n_fn = args.length - getopt.ind;
if (n_fn < 2) {
print("Usage: paftools.js asmgene [options] <ref-splice.paf> <asm-splice.paf> [...]");
print("Options:");
print(" -i FLOAT min identity [" + opt.min_iden + "]");
print(" -c FLOAT min coverage [" + opt.min_cov + "]");
print(" -e print fragmented/missing genes");
exit(1);
}
print(header.join("\t"));
for (var i = 0; i < labels.length; ++i)
print(labels[i], rst[i].join("\t"));
function process_query(opt, a) {
var b = [], cnt = [0, 0, 0];
for (var j = 0; j < a.length; ++j) {
if (a[j][4] < a[j][5] * opt.min_iden)
continue;
b.push(a[j].slice(0));
}
if (b.length == 0) return cnt;
// count full
var n_full = 0;
for (var j = 0; j < b.length; ++j)
if (b[j][3] - b[j][2] >= b[j][1] * opt.min_cov)
++n_full;
cnt[0] = n_full;
// compute coverage
b = b.sort(function(x, y) { return x[2] - y[2] });
var l_cov = 0, st = b[0][2], en = b[0][3];
for (var j = 1; j < b.length; ++j) {
if (b[j][2] <= en)
en = b[j][3] > en? b[j][3] : en;
else l_cov += en - st;
}
l_cov += en - st;
cnt[1] = l_cov / b[0][1];
cnt[2] = b.length;
return cnt;
}
var buf = new Bytes();
var gene = {}, header = [], refpos = {};
for (var i = getopt.ind; i < args.length; ++i) {
var fn = args[i];
var label = fn.replace(/.paf(.gz)?$/, "");
header.push(label);
var file = new File(fn), a = [];
while (file.readline(buf) >= 0) {
var t = buf.toString().split("\t");
var ql = parseInt(t[1]), qs = parseInt(t[2]), qe = parseInt(t[3]), mlen = parseInt(t[9]), blen = parseInt(t[10]), mapq = parseInt(t[11]);
if (i == getopt.ind) refpos[t[0]] = [t[0], t[1], t[5], t[7], t[8]];
if (gene[t[0]] == null) gene[t[0]] = [];
if (a.length && t[0] != a[0][0]) {
gene[t[0]][i - getopt.ind] = process_query(opt, a);
a = [];
}
a.push([t[0], ql, qs, qe, mlen, blen]);
}
if (a.length)
gene[t[0]][i - getopt.ind] = process_query(opt, a);
file.close();
}
// select the longest genes (not optimal, but should be good enough)
var gene_list = [], gene_nr = {};
for (var g in refpos)
gene_list.push(refpos[g]);
gene_list = gene_list.sort(function(a, b) { return a[2] < b[2]? -1 : a[2] > b[2]? 1 : a[3] - b[3] });
var last = 0;
for (var j = 1; j < gene_list.length; ++j) {
if (gene_list[j][2] != gene_list[last][2] || gene_list[j][3] >= gene_list[last][4]) {
gene_nr[gene_list[last][0]] = 1;
last = j;
} else if (gene_list[j][1] > gene_list[last][1]) {
last = j;
}
}
gene_nr[gene_list[last][0]] = 1;
// count and print
var col1 = ["full_sgl", "full_dup", "frag", "part50+", "part10+", "part10-"];
var rst = [];
for (var k = 0; k < col1.length; ++k) {
rst[k] = [];
for (var i = 0; i < n_fn; ++i)
rst[k][i] = 0;
}
for (var g in gene) {
if (gene[g][0] == null || gene[g][0][0] != 1) continue;
if (gene_nr[g] == null) continue;
for (var i = 0; i < n_fn; ++i) {
if (gene[g][i] == null) {
rst[4][i]++;
if (print_err) print('M', header[i], refpos[g].join("\t"));
} else if (gene[g][i][0] == 1) rst[0][i]++;
else if (gene[g][i][0] > 1) {
rst[1][i]++;
if (print_err) print('D', header[i], refpos[g].join("\t"));
} else if (gene[g][i][1] >= opt.min_cov) {
rst[2][i]++;
if (print_err) print('F', header[i], refpos[g].join("\t"));
} else if (gene[g][i][1] >= 0.5) {
rst[3][i]++;
if (print_err) print('5', header[i], refpos[g].join("\t"));
} else if (gene[g][i][1] >= 0.1) {
rst[4][i]++;
if (print_err) print('1', header[i], refpos[g].join("\t"));
} else {
rst[5][i]++;
if (print_err) print('0', header[i], refpos[g].join("\t")); // TODO: reduce code duplicates...
}
}
}
print('H', 'Metric', header.join("\t"));
for (var k = 0; k < rst.length; ++k) {
print('X', col1[k], rst[k].join("\t"));
}
buf.destroy();
}
@@ -2187,6 +2343,7 @@ function main(args)
print("");
print(" stat collect basic mapping information in PAF/SAM");
print(" asmstat collect basic assembly information");
print(" asmgene evaluate gene completeness (EXPERIMENTAL)");
print(" liftover simplistic liftOver");
print(" call call variants from asm-to-ref alignment with the cs tag");
print(" bedcov compute the number of bases covered");
@@ -2208,6 +2365,7 @@ function main(args)
else if (cmd == 'gff2bed') paf_gff2bed(args);
else if (cmd == 'stat') paf_stat(args);
else if (cmd == 'asmstat') paf_asmstat(args);
else if (cmd == 'asmgene') paf_asmgene(args);
else if (cmd == 'liftover' || cmd == 'liftOver') paf_liftover(args);
else if (cmd == 'call') paf_call(args);
else if (cmd == 'mapeval') paf_mapeval(args);
+1 -1
View File
@@ -75,7 +75,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);
+15
View File
@@ -132,6 +132,11 @@ int mm_set_opt(const char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
{
if (mo->split_prefix && (mo->flag & (MM_F_OUT_CS|MM_F_OUT_MD))) {
if (mm_verbose >= 1)
fprintf(stderr, "[ERROR]\033[1;31m --cs or --MD doesn't work with --split-prefix\033[0m\n");
return -6;
}
if (io->k <= 0 || io->w <= 0) {
if (mm_verbose >= 1)
fprintf(stderr, "[ERROR]\033[1;31m -k and -w must be positive\033[0m\n");
@@ -154,6 +159,11 @@ int mm_check_opt(const mm_idxopt_t *io, const mm_mapopt_t *mo)
fprintf(stderr, "[ERROR]\033[1;31m --for-only and --rev-only can't be applied at the same time\033[0m\n");
return -3;
}
if (mo->e <= 0 || mo->q <= 0) {
if (mm_verbose >= 1)
fprintf(stderr, "[ERROR]\033[1;31m -O and -E must be positive\033[0m\n");
return -1;
}
if ((mo->q != mo->q2 || mo->e != mo->e2) && !(mo->e > mo->e2 && mo->q + mo->e < mo->q2 + mo->e2)) {
if (mm_verbose >= 1)
fprintf(stderr, "[ERROR]\033[1;31m dual gap penalties violating E1>E2 and O1+E1<O2+E2\033[0m\n");
@@ -169,5 +179,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;
}
+1 -1
View File
@@ -54,7 +54,7 @@ void mm_set_pe_thru(const int *qlens, int *n_regs, mm_reg1_t **regs)
if (n_pri[0] == 1 && n_pri[1] == 1) {
mm_reg1_t *p = &regs[0][pri[0]];
mm_reg1_t *q = &regs[1][pri[1]];
if (p->rid == q->rid && p->rev == q->rev && abs(p->rs - q->rs) < 3 && abs(p->re - p->re) < 3
if (p->rid == q->rid && p->rev == q->rev && abs(p->rs - q->rs) < 3 && abs(p->re - q->re) < 3
&& ((p->qs == 0 && qlens[1] - q->qe == 0) || (q->qs == 0 && qlens[0] - p->qe == 0)))
{
p->pe_thru = q->pe_thru = 1;
+26 -6
View File
@@ -43,20 +43,24 @@ The following Python script demonstrates the key functionality of mappy:
APIs
----
Mappy implements two classes and one global function.
Mappy implements two classes and two global function.
Class mappy.Aligner
~~~~~~~~~~~~~~~~~~~
.. code:: python
mappy.Aligner(fn_idx_in, preset=None, ...)
mappy.Aligner(fn_idx_in=None, preset=None, ...)
This constructor accepts the following arguments:
* **fn_idx_in**: index or sequence file name. Minimap2 automatically tests the
file type. If a sequence file is provided, minimap2 builds an index. The
sequence file can be optionally gzip'd.
sequence file can be optionally gzip'd. This option has no effect if **seq**
is set.
* **seq**: a single sequence to index. The sequence name will be set to
:code:`N/A`.
* **preset**: minimap2 preset. Currently, minimap2 supports the following
presets: **sr** for single-end short reads; **map-pb** for PacBio
@@ -79,17 +83,28 @@ This constructor accepts the following arguments:
* **n_threads**: number of indexing threads; 3 by default
* **fn_idx_out**: name of file to which the index is written
* **extra_flags**: additional flags defined in minimap.h
* **fn_idx_out**: name of file to which the index is written. This parameter
has no effect if **seq** is set.
* **scoring**: scoring system. It is a tuple/list consisting of 4, 6 or 7
positive integers. The first 4 elements specify match scoring, mismatch
penalty, gap open and gap extension penalty. The 5th and 6th elements, if
present, set long-gap open and long-gap extension penalty. The 7th sets a
mismatch penalty involving ambiguous bases.
.. code:: python
mappy.Aligner.map(seq, seq2=None)
mappy.Aligner.map(seq, seq2=None, cs=False, MD=False)
This method aligns :code:`seq` against the index. It is a generator, *yielding*
a series of :code:`mappy.Alignment` objects. If :code:`seq2` is present, mappy
performs paired-end alignment, assuming the two ends are in the FR orientation.
Alignments of the two ends can be distinguished by the :code:`read_num` field
(see Class mappy.Alignment below).
(see Class mappy.Alignment below). Argument :code:`cs` asks mappy to generate
the :code:`cs` tag; :code:`MD` is similar. These two arguments might slightly
degrade performance and are not enabled by default.
.. code:: python
@@ -139,6 +154,11 @@ properties:
* **cigar**: CIGAR returned as an array of shape :code:`(n_cigar,2)`. The two
numbers give the length and the operator of each CIGAR operation.
* **MD**: the :code:`MD` tag as in the SAM format. It is an empty string unless
the :code:`MD` argument is applied when calling :code:`mappy.Aligner.map()`.
* **cs**: the :code:`cs` tag.
An :code:`Alignment` object can be converted to a string with :code:`str()` in
the following format:
+14 -2
View File
@@ -126,8 +126,8 @@ static char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int e
*len = 0;
rid = mm_idx_name2id(mi, name);
if (rid < 0) return 0;
if (st >= mi->seq[rid].len || st >= en) return 0;
if (en < 0 || en > mi->seq[rid].len)
if ((uint32_t)st >= mi->seq[rid].len || st >= en) return 0;
if (en < 0 || (uint32_t)en > mi->seq[rid].len)
en = mi->seq[rid].len;
s = (char*)malloc(en - st + 1);
*len = mm_idx_getseq(mi, rid, st, en, (uint8_t*)s);
@@ -137,4 +137,16 @@ static char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int e
return s;
}
static mm_idx_t *mappy_idx_seq(int w, int k, int is_hpc, int bucket_bits, const char *seq, int len)
{
const char *fake_name = "N/A";
char *s;
mm_idx_t *mi;
s = (char*)calloc(len + 1, 1);
memcpy(s, seq, len);
mi = mm_idx_str(w, k, is_hpc, bucket_bits, 1, (const char**)&s, (const char**)&fake_name);
free(s);
return mi;
}
#endif
+6
View File
@@ -40,6 +40,8 @@ cdef extern from "minimap.h":
int32_t mid_occ
int32_t max_occ
int mini_batch_size
int64_t max_sw_mat
const char *split_prefix
int mm_set_opt(char *preset, mm_idxopt_t *io, mm_mapopt_t *mo)
int mm_verbose
@@ -86,6 +88,9 @@ cdef extern from "minimap.h":
mm_tbuf_t *mm_tbuf_init()
void mm_tbuf_destroy(mm_tbuf_t *b)
void *mm_tbuf_get_km(mm_tbuf_t *b)
int mm_gen_cs(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq, int no_iden)
int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq)
#
# Helper header (because it is hard to expose mm_reg1_t with Cython)
@@ -106,6 +111,7 @@ cdef extern from "cmappy.h":
void mm_free_reg1(mm_reg1_t *r)
mm_reg1_t *mm_map_aux(const mm_idx_t *mi, const char *seq1, const char *seq2, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt)
char *mappy_fetch_seq(const mm_idx_t *mi, const char *name, int st, int en, int *l)
mm_idx_t *mappy_idx_seq(int w, int k, int is_hpc, int bucket_bits, const char *seq, int l)
ctypedef struct kstring_t:
unsigned l, m
+56 -17
View File
@@ -3,7 +3,7 @@ from libc.stdlib cimport free
cimport cmappy
import sys
__version__ = '2.11'
__version__ = '2.14'
cmappy.mm_reset_timer()
@@ -14,9 +14,9 @@ cdef class Alignment:
cdef int8_t _strand, _trans_strand
cdef uint8_t _mapq, _is_primary
cdef int _seg_id
cdef _ctg, _cigar # these are python objects
cdef _ctg, _cigar, _cs, _MD # these are python objects
def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand, seg_id):
def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand, seg_id, cs_str, MD_str):
self._ctg = ctg if isinstance(ctg, str) else ctg.decode()
self._ctg_len, self._r_st, self._r_en = cl, cs, ce
self._strand, self._q_st, self._q_en = strand, qs, qe
@@ -26,6 +26,8 @@ cdef class Alignment:
self._is_primary = is_primary
self._trans_strand = trans_strand
self._seg_id = seg_id
self._cs = cs_str
self._MD = MD_str
@property
def ctg(self): return self._ctg
@@ -72,6 +74,12 @@ cdef class Alignment:
@property
def read_num(self): return self._seg_id + 1
@property
def cs(self): return self._cs
@property
def MD(self): return self._MD
@property
def cigar_str(self):
return "".join(map(lambda x: str(x[0]) + 'MIDNSH'[x[1]], self._cigar))
@@ -85,8 +93,10 @@ cdef class Alignment:
if self._trans_strand > 0: ts = 'ts:A:+'
elif self._trans_strand < 0: ts = 'ts:A:-'
else: ts = 'ts:A:.'
return "\t".join([str(self._q_st), str(self._q_en), strand, self._ctg, str(self._ctg_len), str(self._r_st), str(self._r_en),
str(self._mlen), str(self._blen), str(self._mapq), tp, ts, "cg:Z:" + self.cigar_str])
a = [str(self._q_st), str(self._q_en), strand, self._ctg, str(self._ctg_len), str(self._r_st), str(self._r_en),
str(self._mlen), str(self._blen), str(self._mapq), tp, ts, "cg:Z:" + self.cigar_str]
if self._cs != "": a.append("cs:Z:" + self._cs)
return "\t".join(a)
cdef class ThreadBuffer:
cdef cmappy.mm_tbuf_t *_b
@@ -102,7 +112,7 @@ cdef class Aligner:
cdef cmappy.mm_idxopt_t idx_opt
cdef cmappy.mm_mapopt_t map_opt
def __cinit__(self, fn_idx_in, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None):
def __cinit__(self, fn_idx_in=None, preset=None, k=None, w=None, min_cnt=None, min_chain_score=None, min_dp_score=None, bw=None, best_n=None, n_threads=3, fn_idx_out=None, max_frag_len=None, extra_flags=None, seq=None, scoring=None):
cmappy.mm_set_opt(NULL, &self.idx_opt, &self.map_opt) # set the default options
if preset is not None:
cmappy.mm_set_opt(str.encode(preset), &self.idx_opt, &self.map_opt) # apply preset
@@ -116,17 +126,32 @@ cdef class Aligner:
if bw is not None: self.map_opt.bw = bw
if best_n is not None: self.map_opt.best_n = best_n
if max_frag_len is not None: self.map_opt.max_frag_len = max_frag_len
if extra_flags is not None: self.map_opt.flag |= extra_flags
if scoring is not None and len(scoring) >= 4:
self.map_opt.a, self.map_opt.b = scoring[0], scoring[1]
self.map_opt.q, self.map_opt.e = scoring[2], scoring[3]
self.map_opt.q2, self.map_opt.e2 = self.map_opt.q, self.map_opt.e
if len(scoring) >= 6:
self.map_opt.q2, self.map_opt.e2 = scoring[4], scoring[5]
if len(scoring) >= 7:
self.map_opt.sc_ambi = scoring[6]
cdef cmappy.mm_idx_reader_t *r;
if fn_idx_out is None:
r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, NULL)
if seq is None:
if fn_idx_out is None:
r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, NULL)
else:
r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, fn_idx_out)
if r is not NULL:
self._idx = cmappy.mm_idx_reader_read(r, n_threads) # NB: ONLY read the first part
cmappy.mm_idx_reader_close(r)
cmappy.mm_mapopt_update(&self.map_opt, self._idx)
cmappy.mm_idx_index_name(self._idx)
else:
r = cmappy.mm_idx_reader_open(str.encode(fn_idx_in), &self.idx_opt, fn_idx_out)
if r is not NULL:
self._idx = cmappy.mm_idx_reader_read(r, n_threads) # NB: ONLY read the first part
cmappy.mm_idx_reader_close(r)
self._idx = cmappy.mappy_idx_seq(self.idx_opt.w, self.idx_opt.k, self.idx_opt.flag&1, self.idx_opt.bucket_bits, str.encode(seq), len(seq))
cmappy.mm_mapopt_update(&self.map_opt, self._idx)
cmappy.mm_idx_index_name(self._idx)
self.map_opt.mid_occ = 1000 # don't filter high-occ seeds
def __dealloc__(self):
if self._idx is not NULL:
@@ -135,18 +160,24 @@ cdef class Aligner:
def __bool__(self):
return (self._idx != NULL)
def map(self, seq, seq2=None, buf=None, max_frag_len=None):
def map(self, seq, seq2=None, buf=None, cs=False, MD=False, max_frag_len=None, extra_flags=None):
cdef cmappy.mm_reg1_t *regs
cdef cmappy.mm_hitpy_t h
cdef ThreadBuffer b
cdef int n_regs
cdef char *cs_str = NULL
cdef int l_cs_str, m_cs_str = 0
cdef void *km
cdef cmappy.mm_mapopt_t map_opt
map_opt = self.map_opt
if max_frag_len is not None: map_opt.max_frag_len = max_frag_len
if extra_flags is not None: map_opt.flag |= extra_flags
if self._idx is NULL: return None
if buf is None: b = ThreadBuffer()
else: b = buf
km = cmappy.mm_tbuf_get_km(b._b)
_seq = seq if isinstance(seq, bytes) else seq.encode()
if seq2 is None:
@@ -157,13 +188,21 @@ cdef class Aligner:
for i in range(n_regs):
cmappy.mm_reg2hitpy(self._idx, &regs[i], &h)
cigar = []
for k in range(h.n_cigar32):
cigar, _cs, _MD = [], '', ''
for k in range(h.n_cigar32): # convert the 32-bit CIGAR encoding to Python array
c = h.cigar32[k]
cigar.append([c>>4, c&0xf])
yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand, h.seg_id)
if cs or MD: # generate the cs and/or the MD tag, if requested
if cs:
l_cs_str = cmappy.mm_gen_cs(km, &cs_str, &m_cs_str, self._idx, &regs[i], _seq, 1)
_cs = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
if MD:
l_cs_str = cmappy.mm_gen_MD(km, &cs_str, &m_cs_str, self._idx, &regs[i], _seq)
_MD = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
yield Alignment(h.ctg, h.ctg_len, h.ctg_start, h.ctg_end, h.strand, h.qry_start, h.qry_end, h.mapq, cigar, h.is_primary, h.mlen, h.blen, h.NM, h.trans_strand, h.seg_id, _cs, _MD)
cmappy.mm_free_reg1(&regs[i])
free(regs)
free(cs_str)
def seq(self, str name, int start=0, int end=0x7fffffff):
cdef int l
+8 -4
View File
@@ -1,10 +1,11 @@
#!/usr/bin/env python
import sys, getopt
import sys
import getopt
import mappy as mp
def main(argv):
opts, args = getopt.getopt(argv[1:], "x:n:m:k:w:r:")
opts, args = getopt.getopt(argv[1:], "x:n:m:k:w:r:c")
if len(args) < 2:
print("Usage: minimap2.py [options] <ref.fa>|<ref.mmi> <query.fq>")
print("Options:")
@@ -14,9 +15,11 @@ def main(argv):
print(" -k INT k-mer length")
print(" -w INT minimizer window length")
print(" -r INT band width")
print(" -c output the cs tag")
sys.exit(1)
preset, min_cnt, min_sc, k, w, bw = None, None, None, None, None, None
preset = min_cnt = min_sc = k = w = bw = None
out_cs = False
for opt, arg in opts:
if opt == '-x': preset = arg
elif opt == '-n': min_cnt = int(arg)
@@ -24,11 +27,12 @@ def main(argv):
elif opt == '-r': bw = int(arg)
elif opt == '-k': k = int(arg)
elif opt == '-w': w = int(arg)
elif opt == '-c': out_cs = True
a = mp.Aligner(args[0], preset=preset, min_cnt=min_cnt, min_chain_score=min_sc, k=k, w=w, bw=bw)
if not a: raise Exception("ERROR: failed to load/build index file '{}'".format(args[0]))
for name, seq, qual in mp.fastx_read(args[1]): # read one sequence
for h in a.map(seq): # traverse hits
for h in a.map(seq, cs=out_cs): # traverse hits
print('{}\t{}\t{}'.format(name, len(seq), h))
if __name__ == "__main__":
+7 -6
View File
@@ -177,7 +177,7 @@ uint64_t *sdust(void *km, const uint8_t *seq, int l_seq, int T, int W, int *n)
#ifdef _SDUST_MAIN
#include <zlib.h>
#include <stdio.h>
#include "getopt.h"
#include "ketopt.h"
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)
@@ -186,16 +186,17 @@ int main(int argc, char *argv[])
gzFile fp;
kseq_t *ks;
int W = 64, T = 20, c;
ketopt_t o = KETOPT_INIT;
while ((c = getopt(argc, argv, "w:t:")) >= 0) {
if (c == 'w') W = atoi(optarg);
else if (c == 't') T = atoi(optarg);
while ((c = ketopt(&o, argc, argv, 1, "w:t:", 0)) >= 0) {
if (c == 'w') W = atoi(o.arg);
else if (c == 't') T = atoi(o.arg);
}
if (optind == argc) {
if (o.ind == argc) {
fprintf(stderr, "Usage: sdust [-w %d] [-t %d] <in.fa>\n", W, T);
return 1;
}
fp = strcmp(argv[optind], "-")? gzopen(argv[optind], "r") : gzdopen(fileno(stdin), "r");
fp = strcmp(argv[o.ind], "-")? gzopen(argv[o.ind], "r") : gzdopen(fileno(stdin), "r");
ks = kseq_init(fp);
while (kseq_read(ks) >= 0) {
uint64_t *r;
+17 -7
View File
@@ -14,16 +14,26 @@ else: # with Cython
module_src = 'python/mappy.pyx'
cmdclass['build_ext'] = build_ext
import sys
import sys, platform
sys.path.append('python')
extra_compile_args = ['-DHAVE_KALLOC']
include_dirs = ["."]
if platform.machine() in ["aarch64", "arm64"]:
include_dirs.append("sse2neon/")
extra_compile_args.extend(['-ftree-vectorize', '-DKSW_SSE2_ONLY', '-D__SSE2__'])
else:
extra_compile_args.append('-msse4.1') # WARNING: ancient x86_64 CPUs don't have SSE4
def readme():
with open('python/README.rst') as f:
return f.read()
with open('python/README.rst') as f:
return f.read()
setup(
name = 'mappy',
version = '2.11',
version = '2.14',
url = 'https://github.com/lh3/minimap2',
description = 'Minimap2 python binding',
long_description = readme(),
@@ -35,12 +45,12 @@ setup(
ext_modules = [Extension('mappy',
sources = [module_src, 'align.c', 'bseq.c', 'chain.c', 'format.c', 'hit.c', 'index.c', 'pe.c', 'options.c',
'ksw2_extd2_sse.c', 'ksw2_exts2_sse.c', 'ksw2_extz2_sse.c', 'ksw2_ll_sse.c',
'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c', 'esterr.c'],
'kalloc.c', 'kthread.c', 'map.c', 'misc.c', 'sdust.c', 'sketch.c', 'esterr.c', 'splitidx.c'],
depends = ['minimap.h', 'bseq.h', 'kalloc.h', 'kdq.h', 'khash.h', 'kseq.h', 'ksort.h',
'ksw2.h', 'kthread.h', 'kvec.h', 'mmpriv.h', 'sdust.h',
'python/cmappy.h', 'python/cmappy.pxd'],
extra_compile_args = ['-DHAVE_KALLOC', '-msse4.1'], # WARNING: ancient x86_64 CPUs don't have SSE4
include_dirs = ['.'],
extra_compile_args = extra_compile_args,
include_dirs = include_dirs,
libraries = ['z', 'm', 'pthread'])],
classifiers = [
'Development Status :: 5 - Production/Stable',