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
- name: Compile with ${{ matrix.compiler }}
uses: uraimo/run-on-arch-action@v2
uses: uraimo/run-on-arch-action@v3
with:
arch: aarch64
distro: ubuntu22.04
@@ -65,4 +65,4 @@ jobs:
run: |
make CC=${{ matrix.compiler }} arm_neon=1 aarch64=1 -j
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
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:
(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)
---------------------------------
+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)
[![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
**ALERT:** `minimap2.com` is a [phishing site](https://github.com/lh3/minimap2/issues/1316). Please don't use anything from that website.
```sh
git clone https://github.com/lh3/minimap2
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
the [release page][release] with:
```sh
curl -L https://github.com/lh3/minimap2/releases/download/v2.30/minimap2-2.30_x64-linux.tar.bz2 | tar -jxvf -
./minimap2-2.30_x64-linux/minimap2
curl -L https://github.com/lh3/minimap2/releases/download/v2.31/minimap2-2.31_x64-linux.tar.bz2 | tar -jxvf -
./minimap2-2.31_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
+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:
```sh
# install minimap2 executables
curl -L https://github.com/lh3/minimap2/releases/download/v2.30/minimap2-2.30_x64-linux.tar.bz2 | tar jxf -
cp minimap2-2.30_x64-linux/{minimap2,k8,paftools.js} . # copy executables
curl -L https://github.com/lh3/minimap2/releases/download/v2.31/minimap2-2.31_x64-linux.tar.bz2 | tar jxf -
cp minimap2-2.31_x64-linux/{minimap2,k8,paftools.js} . # copy executables
export PATH="$PATH:"`pwd` # put the current directory on PATH
# download example datasets
curl -L https://github.com/lh3/minimap2/releases/download/v2.10/cookbook-data.tgz | tar zxf -
+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);
}
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;
kstring_t str;
str.s = *buf, str.l = 0, str.m = *max_len;
t.l_seq = strlen(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;
*buf = str.s;
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)
{
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)
{
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) {
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;
keep[i] = 0;
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) {
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;
else if (r[i].p) free(r[i].p);
keep[i] = 1, ++n_2nd;
} 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[k++] = r[i], ++n_2nd;
} else if (r[i].p) free(r[i].p);
keep[i] = 1, ++n_2nd;
}
}
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()
*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 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;
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];
else ++k;
}
}
free(keep);
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;
for (i = 0; i < slen; ++i)
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 {
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;
for (i = 0; i < slen; ++i)
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;
Submodule lib/simde deleted from b30129b3b4
+2 -1
View File
@@ -5,7 +5,7 @@
#include <stdio.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_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
*/
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);
// 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
.PP
minimap2 - mapping and alignment between collections of DNA sequences
+100 -4
View File
@@ -1,6 +1,6 @@
#!/usr/bin/env k8
var paftools_version = '2.30-r1287';
var paftools_version = '2.31-r1302';
/*****************************
***** Library functions *****
@@ -1740,10 +1740,11 @@ function paf_gff2bed(args)
function paf_sam2paf(args)
{
var c, pri_only = false, long_cs = false, pri_pri_only = false;
while ((c = getopt(args, "pPL")) != null) {
var c, pri_only = false, long_cs = false, pri_pri_only = false, allow_unmapped = false;
while ((c = getopt(args, "pPUL")) != null) {
if (c == 'p') 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;
}
if (args.length == getopt.ind) {
@@ -1751,6 +1752,7 @@ function paf_sam2paf(args)
print("Options:");
print(" -p convert primary or supplementary 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");
exit(1);
}
@@ -1775,7 +1777,15 @@ function paf_sam2paf(args)
var flag = parseInt(t[1]);
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);
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_pri_only && (flag&0x900)) continue;
var tlen = ctg_len[t[2]];
@@ -2338,6 +2348,41 @@ function paf_mason2fq(args)
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
function paf_pbsim2fq(args)
{
@@ -2395,6 +2440,53 @@ function paf_pbsim2fq(args)
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)
{
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(" pafcmp compare two PAF files");
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(" badread2fa convert Baderead FASTQ to FASTA");
print(" junceval evaluate splice junction 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");
@@ -3718,7 +3812,9 @@ function main(args)
else if (cmd == 'pafcmp') paf_pafcmp(args);
else if (cmd == 'bedcov') paf_bedcov(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 == 'badread2fa') paf_badread2fa(args);
else if (cmd == 'junceval') paf_junceval(args);
else if (cmd == 'exoneval') paf_exoneval(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) {
int i, k, n = *n_, n_2nd = 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;
if (r[i].parent == i) { // primary
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 (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);
}
kfree(km, keep);
if (k != n) mm_sync_regs(km, k, r); // removing hits requires sync()
*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_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_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)
#
+16 -8
View File
@@ -3,7 +3,7 @@ from libc.stdlib cimport free
cimport cmappy
import sys
__version__ = '2.30'
__version__ = '2.31'
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, _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_len, self._r_st, self._r_en = cl, cs, ce
self._strand, self._q_st, self._q_en = strand, qs, qe
@@ -27,6 +27,7 @@ cdef class Alignment:
self._trans_strand = trans_strand
self._seg_id = seg_id
self._cs = cs_str
self._ds = ds_str
self._MD = MD_str
@property
@@ -77,6 +78,9 @@ cdef class Alignment:
@property
def cs(self): return self._cs
@property
def ds(self): return self._ds
@property
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),
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._ds != "": a.append("ds:Z:" + self._ds)
if self._MD != "": a.append("MD:Z:" + self._MD)
return "\t".join(a)
@@ -165,7 +170,7 @@ cdef class Aligner:
def __bool__(self):
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_hitpy_t h
cdef ThreadBuffer b
@@ -184,7 +189,6 @@ cdef class Aligner:
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 name is not None:
@@ -206,19 +210,23 @@ cdef class Aligner:
i = 0
while i < n_regs:
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
c = h.cigar32[k]
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
if cs:
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()
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:
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()
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])
i += 1
finally:
+5 -3
View File
@@ -5,7 +5,7 @@ import getopt
import mappy as mp
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:
print("Usage: minimap2.py [options] <ref.fa>|<ref.mmi> <query.fq>")
print("Options:")
@@ -16,11 +16,12 @@ def main(argv):
print(" -w INT minimizer window length")
print(" -r INT band width")
print(" -c output the cs tag")
print(" -d output the ds tag")
print(" -M output the MD tag")
sys.exit(1)
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:
if opt == '-x': preset = arg
elif opt == '-n': min_cnt = int(arg)
@@ -29,12 +30,13 @@ def main(argv):
elif opt == '-k': k = int(arg)
elif opt == '-w': w = int(arg)
elif opt == '-c': out_cs = True
elif opt == '-d': out_ds = 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)
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, 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))
if __name__ == "__main__":
+1 -1
View File
@@ -23,7 +23,7 @@ def readme():
setup(
name = 'mappy',
version = '2.30',
version = '2.31',
url = 'https://github.com/lh3/minimap2',
description = 'Minimap2 python binding',
long_description = readme(),