mirror of
https://github.com/lh3/minimap2.git
synced 2026-09-23 23:58:12 +08:00
Compare commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
3c28777e7e | ||
|
|
d8371a403b | ||
|
|
e5066c7976 | ||
|
|
f8381755f1 | ||
|
|
80d92c686f | ||
|
|
37a650f58c | ||
|
|
4ca5a951ca | ||
|
|
dd5b2c04b8 | ||
|
|
ca6f4dc5e2 | ||
|
|
de3c6ec646 | ||
|
|
e2542e6425 | ||
|
|
9bb4d2bed4 | ||
|
|
6d49eb690f | ||
|
|
bd0cba5012 | ||
|
|
370f3f8236 |
@@ -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
|
||||
|
||||
|
||||
|
||||
@@ -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)
|
||||
|
||||
@@ -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,6 +3,7 @@
|
||||
[](https://pypi.python.org/pypi/mappy)
|
||||
[](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
@@ -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 -
|
||||
|
||||
@@ -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);
|
||||
|
||||
@@ -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
@@ -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;
|
||||
|
||||
-1
Submodule lib/simde deleted from b30129b3b4
@@ -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
@@ -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
@@ -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);
|
||||
|
||||
@@ -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;
|
||||
}
|
||||
|
||||
@@ -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
@@ -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, ®s[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, ®s[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, ®s[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, ®s[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(®s[i])
|
||||
i += 1
|
||||
finally:
|
||||
|
||||
+5
-3
@@ -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__":
|
||||
|
||||
Reference in New Issue
Block a user