fixed a bug on rev strand; added example

This commit is contained in:
Heng Li 2017-09-16 17:51:13 -04:00
parent b22703a354
commit 10bd4079d1
3 changed files with 58 additions and 9 deletions

View File

@ -56,7 +56,7 @@ ksw2_dispatch.o:ksw2_dispatch.c ksw2.h
$(CC) -c $(CFLAGS) $(CPPFLAGS) -DKSW_CPU_DISPATCH $(INCLUDES) $< -o $@
clean:
rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist minimap2.so minimap2.c
rm -fr gmon.out *.o a.out $(PROG) $(PROG_EXTRA) *~ *.a *.dSYM build dist minimap2.so minimap2.c python/minimap2.c minimap2.egg*
depend:
(LC_ALL=C; export LC_ALL; makedepend -Y -- $(CFLAGS) $(CPPFLAGS) -- *.c)

View File

@ -1,4 +1,4 @@
from libc.stdint cimport uint8_t, int32_t, int64_t, uint32_t, uint64_t
from libc.stdint cimport int8_t, uint8_t, int32_t, int64_t, uint32_t, uint64_t
cdef extern from "minimap.h":
#
@ -79,13 +79,14 @@ cdef extern from "minimap.h":
#
cdef extern from "cminimap2.h":
ctypedef struct mm_hitpy_t:
const char *ctg;
int32_t ctg_start, ctg_end;
int32_t qry_start, qry_end;
int32_t blen, NM, ctg_len;
uint8_t strand, mapq, is_primary, trans_strand;
int32_t n_cigar32;
uint32_t *cigar32;
const char *ctg
int32_t ctg_start, ctg_end
int32_t qry_start, qry_end
int32_t blen, NM, ctg_len
uint8_t mapq, is_primary
int8_t strand, trans_strand
int32_t n_cigar32
uint32_t *cigar32
void mm_reg2hitpy(const mm_idx_t *mi, mm_reg1_t *r, mm_hitpy_t *h)
void mm_free_reg1(mm_reg1_t *r)

48
python/mm2-lite.py 100755
View File

@ -0,0 +1,48 @@
#!/usr/bin/env python
import sys, getopt
import minimap2 as mm
def readfq(fp): # multi-line fasta/fastq parser
last = None
while True:
if not last:
for l in fp:
if l[0] in '>@':
last = l[:-1]
break
if not last: break
name, seqs, last = last[1:].split()[0], [], None
for l in fp:
if l[0] in '@+>':
last = l[:-1]
break
seqs.append(l[:-1])
if not last or last[0] != '+':
yield name, ''.join(seqs), None
if not last: break
else:
seq, leng, seqs = ''.join(seqs), 0, []
for l in fp:
seqs.append(l[:-1])
leng += len(l) - 1
if leng >= len(seq):
last = None
yield name, seq, ''.join(seqs);
break
if last:
yield name, seq, None
break
def main(argv):
opts, args = getopt.getopt(argv[1:], "")
if len(args) < 2:
print("Usage: mm2-lite.py <ref.fa>|<ref.mmi> <query.fq>")
sys.exit(1)
a = mm.Aligner(args[0]) # load/build index
for name, seq, qual in readfq(open(args[1])): # read one sequence
for h in a.map(seq): # traverse hits
print('{}\t{}\t{}'.format(name, len(seq), h))
if __name__ == "__main__":
main(sys.argv)