diff --git a/Makefile b/Makefile index d2d30e1..6d48e6f 100644 --- a/Makefile +++ b/Makefile @@ -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) diff --git a/format.c b/format.c index 2e12669..ae8b1f1 100644 --- a/format.c +++ b/format.c @@ -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); diff --git a/lib/simde b/lib/simde deleted file mode 160000 index b30129b..0000000 --- a/lib/simde +++ /dev/null @@ -1 +0,0 @@ -Subproject commit b30129b3b48a6823013da2b309c50a081177b6b8 diff --git a/minimap.h b/minimap.h index 41fe462..92f5673 100644 --- a/minimap.h +++ b/minimap.h @@ -5,7 +5,7 @@ #include #include -#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_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 diff --git a/python/cmappy.pxd b/python/cmappy.pxd index ded06ae..7378f5b 100644 --- a/python/cmappy.pxd +++ b/python/cmappy.pxd @@ -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) # diff --git a/python/mappy.pyx b/python/mappy.pyx index 5dccc6b..d236a1a 100644 --- a/python/mappy.pyx +++ b/python/mappy.pyx @@ -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 @@ -206,19 +211,22 @@ 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 _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: diff --git a/python/minimap2.py b/python/minimap2.py index db6125d..1848d9a 100755 --- a/python/minimap2.py +++ b/python/minimap2.py @@ -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] | ") 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__":