fork/minimap2
1
0
Fork 0
mirror of https://github.com/lh3/minimap2.git synced 2026-09-15 13:07:55 +08:00
Code Issues Packages Projects Releases Wiki Activity

r1290: support ds in Python

Browse Source
This commit is contained in:
Heng Li
2025-12-12 11:39:01 -05:00
parent bd0cba5012
commit 9bb4d2bed4
7 changed files with 35 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile
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)

14
format.c
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);

1
lib/simde

Submodule lib/simde deleted from b30129b3b4

3
minimap.h
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.30-r1290-dirty"
#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
python/cmappy.pxd
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)
# #

20
python/mappy.pyx
View File

@@ -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
@@ -206,19 +211,22 @@ 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
_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:

8
python/minimap2.py
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__":