Compare commits

...
15 Commits
Author SHA1 Message Date
Heng Li 3c28777e7e Release minimap2-2.31 (r1302) 2026-05-19 18:39:04 -04:00
Heng Li d8371a403b removed duplicate check in sim2bed 2026-04-25 00:35:59 -04:00
Heng LiandClaude Opus 4.5 e5066c7976 added sim2bed to convert simulated read names to BED
Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>
2026-04-25 00:22:26 -04:00
Heng Liandgemini-cli f8381755f1 r1299: fix wrong subchain selection
This and the last bug were both reported by Jeremy Wang.

Co-authored-by: gemini-cli <gemini-cli@users.noreply.github.com>
2026-04-24 16:17:45 -04:00
Heng Li 80d92c686f r1298: fixed out-of-bound local alignment
I believe this would not lead to memory violation; it may occasionally include
a few base beyond the intended alignment start/end.
2026-04-24 13:25:19 -04:00
Heng Li 37a650f58c output simulated identity with badread 2026-04-14 12:49:35 -04:00
Heng Li 4ca5a951ca command line help 2026-04-13 23:52:07 -04:00
Heng Li dd5b2c04b8 Merge remote-tracking branch 'remotes/origin/master' 2026-04-13 23:43:43 -04:00
Heng Li ca6f4dc5e2 convert badread output to evaluation format 2026-04-13 23:43:11 -04:00
blawrence-ont de3c6ec646 Fix use-after-free in mappy (#1345)
mm_map_aux() takes in |b._b| which can end up reallocating |km| at the
end of mm_map_frag_core(). Since the address of |km| is cached before
those calls it ends up pointing to freed memory.

This can result in a crash as seen in #1183, however it also happens to
Just Work most of the time since the new allocation often lands at the
same address as the old one. Preloading ASAN or a similar replacement
allocator that doesn't have that behaviour results in a reliable crash.
2026-02-13 12:14:56 -05:00
Heng Li e2542e6425 Merge remote-tracking branch 'remotes/origin/master' 2025-12-12 11:39:33 -05:00
Heng Li 9bb4d2bed4 r1290: support ds in Python 2025-12-12 11:39:01 -05:00
Alex Leonard 6d49eb690f convert unmapped reads in sam2paf (#1324) 2025-10-05 19:43:23 -04:00
Heng Li bd0cba5012 upgrade arm64-linux to v3 2025-09-09 13:52:16 -04:00
Heng Li 370f3f8236 add an alert about phishing 2025-09-09 11:23:09 -04:00
17 changed files with 207 additions and 40 deletions
+2 -2
View File
@@ -35,7 +35,7 @@ jobs:
uses: actions/checkout@v4 uses: actions/checkout@v4
- name: Compile with ${{ matrix.compiler }} - name: Compile with ${{ matrix.compiler }}
uses: uraimo/run-on-arch-action@v2 uses: uraimo/run-on-arch-action@v3
with: with:
arch: aarch64 arch: aarch64
distro: ubuntu22.04 distro: ubuntu22.04
@@ -65,4 +65,4 @@ jobs:
run: | run: |
make CC=${{ matrix.compiler }} arm_neon=1 aarch64=1 -j make CC=${{ matrix.compiler }} arm_neon=1 aarch64=1 -j
file minimap2 | grep arm64 file minimap2 | grep arm64
+1 -1
View File
@@ -102,7 +102,7 @@ ksw2_exts2_neon.o:ksw2_exts2_sse.c ksw2.h kalloc.h
# other non-file targets # other non-file targets
clean: clean:
rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist mappy*.so mappy.c python/mappy.c mappy.egg* rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist mappy*.so mappy.c python/mappy.c mappy.egg* .eggs
depend: depend:
(LC_ALL=C; export LC_ALL; makedepend -Y -- $(CFLAGS) $(CPPFLAGS) -- *.c) (LC_ALL=C; export LC_ALL; makedepend -Y -- $(CFLAGS) $(CPPFLAGS) -- *.c)
+33
View File
@@ -1,3 +1,36 @@
Release 2.31-r1302 (19 May 2026)
--------------------------------
Notable changes to minimap2:
* Bugfix: supplementary and secondary alignments were occasionally flagged
incorrectly.
* Bugfix: Smith-Waterman alignment for inversion alignment led to an
out-of-bound access in rare cases.
Changes to paftools.js:
* New feature: new `sim2bed` subcommand to get a BED file from simulated
reads.
* New feature: new `badread2fa` subcommand to format reads simulated by
the Badread simulator.
Change to the python binding:
* New feature: mappy optionally writes the `ds` tag.
* Bugfix: a use-after-free error (#1345)
The two bugs in minimap2 had existed for years. They were caught by Jeremy Wang
at UNC when he ported minimap2 to Rust. Due to the two bug fixes, this version
occasionally produces alignment different from the last version.
(2.31: 19 May 2026, r1302)
Release 2.30-r1287 (15 June 2025) Release 2.30-r1287 (15 June 2025)
--------------------------------- ---------------------------------
+3 -2
View File
@@ -3,6 +3,7 @@
[![PyPI](https://img.shields.io/pypi/v/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy) [![PyPI](https://img.shields.io/pypi/v/mappy.svg?style=flat)](https://pypi.python.org/pypi/mappy)
[![Build Status](https://github.com/lh3/minimap2/actions/workflows/ci.yaml/badge.svg)](https://github.com/lh3/minimap2/actions) [![Build Status](https://github.com/lh3/minimap2/actions/workflows/ci.yaml/badge.svg)](https://github.com/lh3/minimap2/actions)
## <a name="started"></a>Getting Started ## <a name="started"></a>Getting Started
**ALERT:** `minimap2.com` is a [phishing site](https://github.com/lh3/minimap2/issues/1316). Please don't use anything from that website.
```sh ```sh
git clone https://github.com/lh3/minimap2 git clone https://github.com/lh3/minimap2
cd minimap2 && make cd minimap2 && make
@@ -77,8 +78,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 Minimap2 is optimized for x86-64 CPUs. You can acquire precompiled binaries from
the [release page][release] with: the [release page][release] with:
```sh ```sh
curl -L https://github.com/lh3/minimap2/releases/download/v2.30/minimap2-2.30_x64-linux.tar.bz2 | tar -jxvf - curl -L https://github.com/lh3/minimap2/releases/download/v2.31/minimap2-2.31_x64-linux.tar.bz2 | tar -jxvf -
./minimap2-2.30_x64-linux/minimap2 ./minimap2-2.31_x64-linux/minimap2
``` ```
If you want to compile from the source, you need to have a C compiler, GNU make 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 and zlib development files installed. Then type `make` in the source code
+2 -2
View File
@@ -31,8 +31,8 @@ To acquire the data used in this cookbook and to install minimap2 and paftools,
please follow the command lines below: please follow the command lines below:
```sh ```sh
# install minimap2 executables # install minimap2 executables
curl -L https://github.com/lh3/minimap2/releases/download/v2.30/minimap2-2.30_x64-linux.tar.bz2 | tar jxf - curl -L https://github.com/lh3/minimap2/releases/download/v2.31/minimap2-2.31_x64-linux.tar.bz2 | tar jxf -
cp minimap2-2.30_x64-linux/{minimap2,k8,paftools.js} . # copy executables cp minimap2-2.31_x64-linux/{minimap2,k8,paftools.js} . # copy executables
export PATH="$PATH:"`pwd` # put the current directory on PATH export PATH="$PATH:"`pwd` # put the current directory on PATH
# download example datasets # download example datasets
curl -L https://github.com/lh3/minimap2/releases/download/v2.10/cookbook-data.tgz | tar zxf - curl -L https://github.com/lh3/minimap2/releases/download/v2.10/cookbook-data.tgz | tar zxf -
+12 -2
View File
@@ -361,24 +361,34 @@ static void write_cs_ds_or_MD(void *km, kstring_t *s, const mm_idx_t *mi, const
kfree(km, qseq); kfree(km, tseq); kfree(km, tmp); 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, int is_qstrand) int mm_gen_cs_ds_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 is_ds, int no_iden, int is_qstrand)
{ {
mm_bseq1_t t; mm_bseq1_t t;
kstring_t str; kstring_t str;
str.s = *buf, str.l = 0, str.m = *max_len; str.s = *buf, str.l = 0, str.m = *max_len;
t.l_seq = strlen(seq); t.l_seq = strlen(seq);
t.seq = (char*)seq; t.seq = (char*)seq;
write_cs_ds_or_MD(km, &str, mi, &t, r, no_iden, is_MD, 0, 0, is_qstrand); write_cs_ds_or_MD(km, &str, mi, &t, r, no_iden, is_MD, is_ds, 0, is_qstrand);
*max_len = str.m; *max_len = str.m;
*buf = str.s; *buf = str.s;
return str.l; return str.l;
} }
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, int is_qstrand)
{
return mm_gen_cs_ds_or_MD(km, buf, max_len, mi, r, seq, is_MD, 0, no_iden, is_qstrand);
}
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_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, 0); return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 0, no_iden, 0);
} }
int mm_gen_ds(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_ds_or_MD(km, buf, max_len, mi, r, seq, 0, 1, no_iden, 0);
}
int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq) 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, 0); return mm_gen_cs_or_MD(km, buf, max_len, mi, r, seq, 1, 0, 0);
+19 -8
View File
@@ -256,19 +256,25 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int chec
{ {
if (pri_ratio > 0.0f && *n_ > 0) { if (pri_ratio > 0.0f && *n_ > 0) {
int i, k, n = *n_, n_2nd = 0; int i, k, n = *n_, n_2nd = 0;
for (i = k = 0; i < n; ++i) { uint8_t *keep = (uint8_t*)kmalloc(km, n);
for (i = 0; i < n; ++i) {
int p = r[i].parent; int p = r[i].parent;
keep[i] = 0;
if (p == i || r[i].inv) { // primary or inversion if (p == i || r[i].inv) { // primary or inversion
r[k++] = r[i]; keep[i] = 1;
} else if ((r[i].score >= r[p].score * pri_ratio || r[i].score + min_diff >= r[p].score) && n_2nd < best_n) { } else if ((r[i].score >= r[p].score * pri_ratio || r[i].score + min_diff >= r[p].score) && n_2nd < best_n) {
if (!(r[i].qs == r[p].qs && r[i].qe == r[p].qe && r[i].rid == r[p].rid && r[i].rs == r[p].rs && r[i].re == r[p].re)) // not identical hits if (!(r[i].qs == r[p].qs && r[i].qe == r[p].qe && r[i].rid == r[p].rid && r[i].rs == r[p].rs && r[i].re == r[p].re)) // not identical hits
r[k++] = r[i], ++n_2nd; keep[i] = 1, ++n_2nd;
else if (r[i].p) free(r[i].p);
} else if (check_strand && n_2nd < best_n && r[i].score > min_strand_sc && r[i].rev != r[p].rev) { } else if (check_strand && n_2nd < best_n && r[i].score > min_strand_sc && r[i].rev != r[p].rev) {
r[i].strand_retained = 1; r[i].strand_retained = 1;
r[k++] = r[i], ++n_2nd; keep[i] = 1, ++n_2nd;
} else if (r[i].p) free(r[i].p); }
} }
for (i = k = 0; i < n; ++i) {
if (keep[i]) r[k++] = r[i];
else if (r[i].p) free(r[i].p);
}
kfree(km, keep);
if (k != n) mm_sync_regs(km, k, r); // removing hits requires sync() if (k != n) mm_sync_regs(km, k, r); // removing hits requires sync()
*n_ = k; *n_ = k;
} }
@@ -277,13 +283,18 @@ void mm_select_sub(void *km, float pri_ratio, int min_diff, int best_n, int chec
int mm_filter_strand_retained(int n_regs, mm_reg1_t *r) int mm_filter_strand_retained(int n_regs, mm_reg1_t *r)
{ {
int i, k; int i, k;
for (i = k = 0; i < n_regs; ++i) { uint8_t *keep = (uint8_t*)malloc(n_regs);
for (i = 0; i < n_regs; ++i) {
int p = r[i].parent; int p = r[i].parent;
if (!r[i].strand_retained || r[i].div < r[p].div * 5.0f || r[i].div < 0.01f) { keep[i] = (!r[i].strand_retained || r[i].div < r[p].div * 5.0f || r[i].div < 0.01f);
}
for (i = k = 0; i < n_regs; ++i) {
if (keep[i]) {
if (k < i) r[k++] = r[i]; if (k < i) r[k++] = r[i];
else ++k; else ++k;
} }
} }
free(keep);
return k; return k;
} }
+2 -2
View File
@@ -67,7 +67,7 @@ void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, co
const int8_t *ma = mat + a * m; const int8_t *ma = mat + a * m;
for (i = 0; i < slen; ++i) for (i = 0; i < slen; ++i)
for (k = i; k < nlen; k += slen) // p iterations for (k = i; k < nlen; k += slen) // p iterations
*t++ = (k >= qlen? 0 : ma[query[k]]) + q->shift; *t++ = (k >= qlen? -1 : ma[query[k]]) + q->shift;
} }
} else { } else {
int16_t *t = (int16_t*)q->qp; int16_t *t = (int16_t*)q->qp;
@@ -76,7 +76,7 @@ void *ksw_ll_qinit(void *km, int size, int qlen, const uint8_t *query, int m, co
const int8_t *ma = mat + a * m; const int8_t *ma = mat + a * m;
for (i = 0; i < slen; ++i) for (i = 0; i < slen; ++i)
for (k = i; k < nlen; k += slen) // p iterations for (k = i; k < nlen; k += slen) // p iterations
*t++ = (k >= qlen? 0 : ma[query[k]]); *t++ = (k >= qlen? -1 : ma[query[k]]);
} }
} }
return q; return q;
Submodule lib/simde deleted from b30129b3b4
+2 -1
View File
@@ -5,7 +5,7 @@
#include <stdio.h> #include <stdio.h>
#include <sys/types.h> #include <sys/types.h>
#define MM_VERSION "2.30-r1287" #define MM_VERSION "2.31-r1302"
#define MM_F_NO_DIAG (0x001LL) // no exact diagonal hit #define MM_F_NO_DIAG (0x001LL) // no exact diagonal hit
#define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name #define MM_F_NO_DUAL (0x002LL) // skip pairs where query name is lexicographically larger than target name
@@ -408,6 +408,7 @@ int mm_map_file_frag(const mm_idx_t *idx, int n_segs, const char **fn, const mm_
* @return the length of cs * @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_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_ds(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); 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 // query sequence name and sequence in the minimap2 index
+1 -1
View File
@@ -1,4 +1,4 @@
.TH minimap2 1 "15 June 2025" "minimap2-2.30 (r1287)" "Bioinformatics tools" .TH minimap2 1 "19 May 2026" "minimap2-2.31 (r1302)" "Bioinformatics tools"
.SH NAME .SH NAME
.PP .PP
minimap2 - mapping and alignment between collections of DNA sequences minimap2 - mapping and alignment between collections of DNA sequences
+100 -4
View File
@@ -1,6 +1,6 @@
#!/usr/bin/env k8 #!/usr/bin/env k8
var paftools_version = '2.30-r1287'; var paftools_version = '2.31-r1302';
/***************************** /*****************************
***** Library functions ***** ***** Library functions *****
@@ -1740,10 +1740,11 @@ function paf_gff2bed(args)
function paf_sam2paf(args) function paf_sam2paf(args)
{ {
var c, pri_only = false, long_cs = false, pri_pri_only = false; var c, pri_only = false, long_cs = false, pri_pri_only = false, allow_unmapped = false;
while ((c = getopt(args, "pPL")) != null) { while ((c = getopt(args, "pPUL")) != null) {
if (c == 'p') pri_only = true; if (c == 'p') pri_only = true;
else if (c == 'P') pri_pri_only = pri_only = true; else if (c == 'P') pri_pri_only = pri_only = true;
else if (c == 'U') allow_unmapped = true;
else if (c == 'L') long_cs = true; else if (c == 'L') long_cs = true;
} }
if (args.length == getopt.ind) { if (args.length == getopt.ind) {
@@ -1751,6 +1752,7 @@ function paf_sam2paf(args)
print("Options:"); print("Options:");
print(" -p convert primary or supplementary alignments only"); print(" -p convert primary or supplementary alignments only");
print(" -P convert primary alignments only"); print(" -P convert primary alignments only");
print(" -U convert unmapped reads as well");
print(" -L output the cs tag in the long form"); print(" -L output the cs tag in the long form");
exit(1); exit(1);
} }
@@ -1775,7 +1777,15 @@ function paf_sam2paf(args)
var flag = parseInt(t[1]); var flag = parseInt(t[1]);
if (t[9] != '*' && t[10] != '*' && t[9].length != t[10].length) if (t[9] != '*' && t[10] != '*' && t[9].length != t[10].length)
throw Error("at line " + lineno + ": inconsistent SEQ and QUAL lengths - " + t[9].length + " != " + t[10].length); throw Error("at line " + lineno + ": inconsistent SEQ and QUAL lengths - " + t[9].length + " != " + t[10].length);
if (t[2] == '*' || (flag&4) || t[5] == '*') continue; if (t[2] == '*' || (flag&4) || t[5] == '*') {
if (allow_unmapped) {
// emit an unmapped PAF line instead of skipping
// fields: qname, qlen, qstart, qend, strand, tname, tlen, tstart, tend, n_match, aln_len, mapq
var qlen_val = (t[9] == '*' ? 0 : t[9].length);
print([t[0], qlen_val, 0, 0, '*', '*', 0, 0, 0, 0, 0, 0].join("\t"));
}
continue;
}
if (pri_only && (flag&0x100)) continue; if (pri_only && (flag&0x100)) continue;
if (pri_pri_only && (flag&0x900)) continue; if (pri_pri_only && (flag&0x900)) continue;
var tlen = ctg_len[t[2]]; var tlen = ctg_len[t[2]];
@@ -2338,6 +2348,41 @@ function paf_mason2fq(args)
buf2.destroy(); buf2.destroy();
} }
// convert Mason read names to BED
function paf_sim2bed(args)
{
if (args.length == 0) {
print("Usage: paftools.js sim2bed <sim.txt>");
exit(1);
}
var buf = new Bytes();
var file = new File(args[0]);
while (file.readline(buf) >= 0) {
var line = buf.toString();
var t = line.split("!");
if (t.length < 5) continue;
var chr = t[1], st, en, strand;
if (t[2].indexOf("_") >= 0) { // mason paired-end
var pos = t[2].split("_");
var end = t[3].split("_");
var m = /^(.)(.)\/([12])$/.exec(t[4]);
if (m == null) continue;
strand = m[3] == "1" ? m[1] : m[2];
var read_no = parseInt(m[3]) - 1;
st = parseInt(pos[read_no]);
en = parseInt(end[read_no]);
} else { // badread/pbsim long reads
st = parseInt(t[2]);
en = parseInt(t[3]);
strand = t[4];
}
if (st > en) { var tmp = st; st = en; en = tmp; }
print([chr, st, en, line, 0, strand].join("\t"));
}
file.close();
buf.destroy();
}
// convert pbsim MAF to FASTQ // convert pbsim MAF to FASTQ
function paf_pbsim2fq(args) function paf_pbsim2fq(args)
{ {
@@ -2395,6 +2440,53 @@ function paf_pbsim2fq(args)
buf2.destroy(); buf2.destroy();
} }
function paf_badread2fa(args)
{
if (args.length < 2) {
print("Usage: paftools.js badread2fa <ref.fa.fai> <badread.fq>");
exit(1);
}
var len = {}, file, buf = new Bytes();
file = new File(args[0]);
while (file.readline(buf) >= 0) {
var t = buf.toString().split("\t");
len[t[0]] = parseInt(t[1]);
}
file.close();
var id = 0, n_discard = 0;
file = new File(args[1]);
while (file.readline(buf) >= 0) {
var line = buf.toString();
var m, tag = '', a = null, is_fq = line[0] == '@'? true : false;
if (!/\schimera\s/.test(line) && (m = /\s(\S+),([+-])strand,(\d+)-(\d+).*read_identity=([0-9\.]+)%/.exec(line)) != null) {
if (len[m[1]] == null) throw Error("failed to find the contig length of " + m[1]);
m[3] = parseInt(m[3]);
m[4] = parseInt(m[4]);
if (m[2] == '+')
a = [ "S" + (id+1), m[1], m[3], m[4], m[2] ];
else
a = [ "S" + (id+1), m[1], len[m[1]] - m[4], len[m[1]] - m[3], m[2] ];
tag = "ri:f:" + m[5];
}
file.readline(buf);
var seq = buf.toString();
if (is_fq) {
file.readline(buf);
file.readline(buf);
}
if (a != null) {
print(">" + a.join("!"), tag);
print(seq);
} else ++n_discard;
++id;
}
file.close();
buf.destroy();
warn("WARNING: discarded " + n_discard + " reads");
}
function paf_junceval(args) function paf_junceval(args)
{ {
var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false, chr_only = false, aa = false, is_bed = false; var c, l_fuzzy = 0, print_ovlp = false, print_err_only = false, first_only = false, chr_only = false, aa = false, is_bed = false;
@@ -3693,7 +3785,9 @@ function main(args)
print(" mapeval evaluate mapping accuracy using mason2/PBSIM-simulated FASTQ"); print(" mapeval evaluate mapping accuracy using mason2/PBSIM-simulated FASTQ");
print(" pafcmp compare two PAF files"); print(" pafcmp compare two PAF files");
print(" mason2fq convert mason2-simulated SAM to FASTQ"); print(" mason2fq convert mason2-simulated SAM to FASTQ");
print(" sim2bed convert mason2-simulated read names to BED");
print(" pbsim2fq convert PBSIM-simulated MAF to FASTQ"); print(" pbsim2fq convert PBSIM-simulated MAF to FASTQ");
print(" badread2fa convert Baderead FASTQ to FASTA");
print(" junceval evaluate splice junction consistency with known annotations"); print(" junceval evaluate splice junction consistency with known annotations");
print(" exoneval evaluate exon-level consistency with known annotations"); print(" exoneval evaluate exon-level consistency with known annotations");
print(" ov-eval evaluate read overlap sensitivity using read-to-ref mapping"); print(" ov-eval evaluate read overlap sensitivity using read-to-ref mapping");
@@ -3718,7 +3812,9 @@ function main(args)
else if (cmd == 'pafcmp') paf_pafcmp(args); else if (cmd == 'pafcmp') paf_pafcmp(args);
else if (cmd == 'bedcov') paf_bedcov(args); else if (cmd == 'bedcov') paf_bedcov(args);
else if (cmd == 'mason2fq') paf_mason2fq(args); else if (cmd == 'mason2fq') paf_mason2fq(args);
else if (cmd == 'sim2bed') paf_sim2bed(args);
else if (cmd == 'pbsim2fq') paf_pbsim2fq(args); else if (cmd == 'pbsim2fq') paf_pbsim2fq(args);
else if (cmd == 'badread2fa') paf_badread2fa(args);
else if (cmd == 'junceval') paf_junceval(args); else if (cmd == 'junceval') paf_junceval(args);
else if (cmd == 'exoneval') paf_exoneval(args); else if (cmd == 'exoneval') paf_exoneval(args);
else if (cmd == 'ov-eval') paf_ov_eval(args); else if (cmd == 'ov-eval') paf_ov_eval(args);
+7 -2
View File
@@ -8,7 +8,8 @@ void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int
if (pri_ratio > 0.0f && *n_ > 0) { if (pri_ratio > 0.0f && *n_ > 0) {
int i, k, n = *n_, n_2nd = 0; int i, k, n = *n_, n_2nd = 0;
int max_dist = n_segs == 2? qlens[0] + qlens[1] + max_gap_ref : 0; int max_dist = n_segs == 2? qlens[0] + qlens[1] + max_gap_ref : 0;
for (i = k = 0; i < n; ++i) { uint8_t *keep = (uint8_t*)kmalloc(km, n);
for (i = 0; i < n; ++i) {
int to_keep = 0; int to_keep = 0;
if (r[i].parent == i) { // primary if (r[i].parent == i) { // primary
to_keep = 1; to_keep = 1;
@@ -34,9 +35,13 @@ void mm_select_sub_multi(void *km, float pri_ratio, float pri1, float pri2, int
if (to_keep && r[i].parent != i) { if (to_keep && r[i].parent != i) {
if (n_2nd++ >= best_n) to_keep = 0; // don't keep if there are too many secondary hits if (n_2nd++ >= best_n) to_keep = 0; // don't keep if there are too many secondary hits
} }
if (to_keep) r[k++] = r[i]; keep[i] = to_keep;
}
for (i = k = 0; i < n; ++i) {
if (keep[i]) r[k++] = r[i];
else if (r[i].p) free(r[i].p); else if (r[i].p) free(r[i].p);
} }
kfree(km, keep);
if (k != n) mm_sync_regs(km, k, r); // removing hits requires sync() if (k != n) mm_sync_regs(km, k, r); // removing hits requires sync()
*n_ = k; *n_ = k;
} }
+1
View File
@@ -112,6 +112,7 @@ cdef extern from "minimap.h":
void mm_tbuf_destroy(mm_tbuf_t *b) void mm_tbuf_destroy(mm_tbuf_t *b)
void *mm_tbuf_get_km(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_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_ds(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) int mm_gen_MD(void *km, char **buf, int *max_len, const mm_idx_t *mi, const mm_reg1_t *r, const char *seq)
# #
+16 -8
View File
@@ -3,7 +3,7 @@ from libc.stdlib cimport free
cimport cmappy cimport cmappy
import sys import sys
__version__ = '2.30' __version__ = '2.31'
cmappy.mm_reset_timer() cmappy.mm_reset_timer()
@@ -14,9 +14,9 @@ cdef class Alignment:
cdef int8_t _strand, _trans_strand cdef int8_t _strand, _trans_strand
cdef uint8_t _mapq, _is_primary cdef uint8_t _mapq, _is_primary
cdef int _seg_id cdef int _seg_id
cdef _ctg, _cigar, _cs, _MD # these are python objects cdef _ctg, _cigar, _cs, _ds, _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, cs_str, MD_str): def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand, seg_id, cs_str, ds_str, MD_str):
self._ctg = ctg if isinstance(ctg, str) else ctg.decode() self._ctg = ctg if isinstance(ctg, str) else ctg.decode()
self._ctg_len, self._r_st, self._r_en = cl, cs, ce self._ctg_len, self._r_st, self._r_en = cl, cs, ce
self._strand, self._q_st, self._q_en = strand, qs, qe self._strand, self._q_st, self._q_en = strand, qs, qe
@@ -27,6 +27,7 @@ cdef class Alignment:
self._trans_strand = trans_strand self._trans_strand = trans_strand
self._seg_id = seg_id self._seg_id = seg_id
self._cs = cs_str self._cs = cs_str
self._ds = ds_str
self._MD = MD_str self._MD = MD_str
@property @property
@@ -77,6 +78,9 @@ cdef class Alignment:
@property @property
def cs(self): return self._cs def cs(self): return self._cs
@property
def ds(self): return self._ds
@property @property
def MD(self): return self._MD def MD(self): return self._MD
@@ -96,6 +100,7 @@ cdef class Alignment:
a = [str(self._q_st), str(self._q_en), strand, self._ctg, str(self._ctg_len), str(self._r_st), str(self._r_en), 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] 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) if self._cs != "": a.append("cs:Z:" + self._cs)
if self._ds != "": a.append("ds:Z:" + self._ds)
if self._MD != "": a.append("MD:Z:" + self._MD) if self._MD != "": a.append("MD:Z:" + self._MD)
return "\t".join(a) return "\t".join(a)
@@ -165,7 +170,7 @@ cdef class Aligner:
def __bool__(self): def __bool__(self):
return (self._idx != NULL) return (self._idx != NULL)
def map(self, seq, seq2=None, name=None, buf=None, cs=False, MD=False, max_frag_len=None, extra_flags=None): def map(self, seq, seq2=None, name=None, buf=None, cs=False, ds=False, MD=False, max_frag_len=None, extra_flags=None):
cdef cmappy.mm_reg1_t *regs cdef cmappy.mm_reg1_t *regs
cdef cmappy.mm_hitpy_t h cdef cmappy.mm_hitpy_t h
cdef ThreadBuffer b cdef ThreadBuffer b
@@ -184,7 +189,6 @@ cdef class Aligner:
if self._idx is NULL: return None if self._idx is NULL: return None
if buf is None: b = ThreadBuffer() if buf is None: b = ThreadBuffer()
else: b = buf else: b = buf
km = cmappy.mm_tbuf_get_km(b._b)
_seq = seq if isinstance(seq, bytes) else seq.encode() _seq = seq if isinstance(seq, bytes) else seq.encode()
if name is not None: if name is not None:
@@ -206,19 +210,23 @@ cdef class Aligner:
i = 0 i = 0
while i < n_regs: while i < n_regs:
cmappy.mm_reg2hitpy(self._idx, &regs[i], &h) cmappy.mm_reg2hitpy(self._idx, &regs[i], &h)
cigar, _cs, _MD = [], '', '' cigar, _cs, _ds, _MD = [], '', '', ''
for k in range(h.n_cigar32): # convert the 32-bit CIGAR encoding to Python array for k in range(h.n_cigar32): # convert the 32-bit CIGAR encoding to Python array
c = h.cigar32[k] c = h.cigar32[k]
cigar.append([c>>4, c&0xf]) cigar.append([c>>4, c&0xf])
if cs or MD: # generate the cs and/or the MD tag, if requested if cs or ds or MD: # generate the cs/ds and/or the MD tag, if requested
km = cmappy.mm_tbuf_get_km(b._b)
_cur_seq = _seq2 if h.seg_id > 0 and seq2 is not None else _seq _cur_seq = _seq2 if h.seg_id > 0 and seq2 is not None else _seq
if cs: if cs:
l_cs_str = cmappy.mm_gen_cs(km, &cs_str, &m_cs_str, self._idx, &regs[i], _cur_seq, 1) l_cs_str = cmappy.mm_gen_cs(km, &cs_str, &m_cs_str, self._idx, &regs[i], _cur_seq, 1)
_cs = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode() _cs = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
if ds:
l_cs_str = cmappy.mm_gen_ds(km, &cs_str, &m_cs_str, self._idx, &regs[i], _cur_seq, 1)
_ds = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode()
if MD: if MD:
l_cs_str = cmappy.mm_gen_MD(km, &cs_str, &m_cs_str, self._idx, &regs[i], _cur_seq) l_cs_str = cmappy.mm_gen_MD(km, &cs_str, &m_cs_str, self._idx, &regs[i], _cur_seq)
_MD = cs_str[:l_cs_str] if isinstance(cs_str, str) else cs_str[:l_cs_str].decode() _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) 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, _ds, _MD)
cmappy.mm_free_reg1(&regs[i]) cmappy.mm_free_reg1(&regs[i])
i += 1 i += 1
finally: finally:
+5 -3
View File
@@ -5,7 +5,7 @@ import getopt
import mappy as mp import mappy as mp
def main(argv): def main(argv):
opts, args = getopt.getopt(argv[1:], "x:n:m:k:w:r:cM") opts, args = getopt.getopt(argv[1:], "x:n:m:k:w:r:cdM")
if len(args) < 2: if len(args) < 2:
print("Usage: minimap2.py [options] <ref.fa>|<ref.mmi> <query.fq>") print("Usage: minimap2.py [options] <ref.fa>|<ref.mmi> <query.fq>")
print("Options:") print("Options:")
@@ -16,11 +16,12 @@ def main(argv):
print(" -w INT minimizer window length") print(" -w INT minimizer window length")
print(" -r INT band width") print(" -r INT band width")
print(" -c output the cs tag") print(" -c output the cs tag")
print(" -d output the ds tag")
print(" -M output the MD tag") print(" -M output the MD tag")
sys.exit(1) sys.exit(1)
preset = min_cnt = min_sc = k = w = bw = None preset = min_cnt = min_sc = k = w = bw = None
out_cs = out_MD = False out_cs = out_ds = out_MD = False
for opt, arg in opts: for opt, arg in opts:
if opt == '-x': preset = arg if opt == '-x': preset = arg
elif opt == '-n': min_cnt = int(arg) elif opt == '-n': min_cnt = int(arg)
@@ -29,12 +30,13 @@ def main(argv):
elif opt == '-k': k = int(arg) elif opt == '-k': k = int(arg)
elif opt == '-w': w = int(arg) elif opt == '-w': w = int(arg)
elif opt == '-c': out_cs = True elif opt == '-c': out_cs = True
elif opt == '-d': out_ds = True
elif opt == '-M': out_MD = True elif opt == '-M': out_MD = True
a = mp.Aligner(args[0], preset=preset, min_cnt=min_cnt, min_chain_score=min_sc, k=k, w=w, bw=bw) 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])) 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 name, seq, qual in mp.fastx_read(args[1]): # read one sequence
for h in a.map(seq, cs=out_cs, MD=out_MD): # traverse hits for h in a.map(seq, cs=out_cs, ds=out_ds, MD=out_MD): # traverse hits
print('{}\t{}\t{}'.format(name, len(seq), h)) print('{}\t{}\t{}'.format(name, len(seq), h))
if __name__ == "__main__": if __name__ == "__main__":
+1 -1
View File
@@ -23,7 +23,7 @@ def readme():
setup( setup(
name = 'mappy', name = 'mappy',
version = '2.30', version = '2.31',
url = 'https://github.com/lh3/minimap2', url = 'https://github.com/lh3/minimap2',
description = 'Minimap2 python binding', description = 'Minimap2 python binding',
long_description = readme(), long_description = readme(),