r670: added PE support to mappy

and minor code cleanup

python/README.rst

@@ -81,10 +81,13 @@ This constructor accepts the following arguments:
 
 .. code:: python
 
-	mappy.Aligner.map(seq)
+	mappy.Aligner.map(seq, seq2=None)
 
 This method aligns :code:`seq` against the index. It is a generator, *yielding*
-a series of :code:`mappy.Alignment` objects.
+a series of :code:`mappy.Alignment` objects. If :code:`seq2` is present, mappy
+performs paired-end alignment, assuming the two ends are in the FR orientation.
+Alignments of the two ends can be distinguished by the :code:`read_num` field
+(see below).
 
 Class mappy.Alignment
 ~~~~~~~~~~~~~~~~~~~~~
@@ -118,6 +121,9 @@ properties:
 * **is_primary**: if the alignment is primary (typically the best and the first
   to generate)
 
+* **read_num**: read number that the alignment corresponds to; 1 for the first
+  read and 2 for the second read
+
 * **cigar_str**: CIGAR string
 
 * **cigar**: CIGAR returned as an array of shape :code:`(n_cigar,2)`. The two

python/cmappy.h

@@ -15,6 +15,7 @@ typedef struct {
 	int32_t blen, mlen, NM, ctg_len;
 	uint8_t mapq, is_primary;
 	int8_t strand, trans_strand;
+	int32_t seg_id;
 	int32_t n_cigar32;
 	uint32_t *cigar32;
 } mm_hitpy_t;
@@ -32,6 +33,7 @@ static inline void mm_reg2hitpy(const mm_idx_t *mi, mm_reg1_t *r, mm_hitpy_t *h)
 	h->NM = r->blen - r->mlen + r->p->n_ambi;
 	h->trans_strand = r->p->trans_strand == 1? 1 : r->p->trans_strand == 2? -1 : 0;
 	h->is_primary = (r->id == r->parent);
+	h->seg_id = r->seg_id;
 	h->n_cigar32 = r->p->n_cigar;
 	h->cigar32 = r->p->cigar;
 }
@@ -68,4 +70,35 @@ static inline void mm_reset_timer(void)
 	mm_realtime0 = realtime();
 }
 
+extern unsigned char seq_comp_table[256];
+static inline mm_reg1_t *mm_map_aux(const mm_idx_t *mi, const char *seq1, const char *seq2, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt)
+{
+	if (seq2 == 0) {
+		return mm_map(mi, strlen(seq1), seq1, n_regs, b, opt, NULL);
+	} else {
+		int _n_regs[2];
+		mm_reg1_t *regs[2];
+		char *seq[2];
+		int i, len[2];
+		len[0] = strlen(seq1);
+		len[1] = strlen(seq2);
+		seq[0] = (char*)seq1;
+		seq[1] = strdup(seq2);
+		for (i = 0; i < len[1]>>1; ++i) {
+			int t = seq[1][len[1] - i - 1];
+			seq[1][len[1] - i - 1] = seq_comp_table[(uint8_t)seq[1][i]];
+			seq[1][i] = seq_comp_table[t];
+		}
+		if (len[1]&1) seq[1][len[1]>>1] = seq_comp_table[(uint8_t)seq[1][len[1]>>1]];
+		mm_map_frag(mi, 2, len, (const char**)seq, _n_regs, regs, b, opt, NULL);
+		for (i = 0; i < _n_regs[1]; ++i)
+			regs[1][i].rev = !regs[1][i].rev;
+		*n_regs = _n_regs[0] + _n_regs[1];
+		regs[0] = (mm_reg1_t*)realloc(regs[0], sizeof(mm_reg1_t) * (*n_regs));
+		memcpy(&regs[0][_n_regs[0]], regs[1], _n_regs[1] * sizeof(mm_reg1_t));
+		free(regs[1]);
+		return regs[0];
+	}
+}
+
 #endif

python/cmappy.pxd

@@ -79,7 +79,6 @@ cdef extern from "minimap.h":
 
 	mm_tbuf_t *mm_tbuf_init()
 	void mm_tbuf_destroy(mm_tbuf_t *b)
-	mm_reg1_t *mm_map(const mm_idx_t *mi, int l_seq, const char *seq, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt, const char *name)
 
 #
 # Helper header (because it is hard to expose mm_reg1_t with Cython)
@@ -92,11 +91,13 @@ cdef extern from "cmappy.h":
 		int32_t blen, mlen, NM, ctg_len
 		uint8_t mapq, is_primary
 		int8_t strand, trans_strand
+		int32_t seg_id
 		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)
+	mm_reg1_t *mm_map_aux(const mm_idx_t *mi, const char *seq1, const char *seq2, int *n_regs, mm_tbuf_t *b, const mm_mapopt_t *opt)
 
 	ctypedef struct kstring_t:
 		unsigned l, m

python/mappy.pyx

@@ -10,9 +10,10 @@ cdef class Alignment:
 	cdef int _NM, _mlen, _blen
 	cdef int8_t _strand, _trans_strand
 	cdef uint8_t _mapq, _is_primary
+	cdef int _seg_id
 	cdef _ctg, _cigar # these are python objects
 
-	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand):
+	def __cinit__(self, ctg, cl, cs, ce, strand, qs, qe, mapq, cigar, is_primary, mlen, blen, NM, trans_strand, seg_id):
 		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
@@ -21,6 +22,7 @@ cdef class Alignment:
 		self._cigar = cigar
 		self._is_primary = is_primary
 		self._trans_strand = trans_strand
+		self._seg_id = seg_id
 
 	@property
 	def ctg(self): return self._ctg
@@ -64,6 +66,9 @@ cdef class Alignment:
 	@property
 	def cigar(self): return self._cigar
 
+	@property
+	def read_num(self): return self._seg_id + 1
+
 	@property
 	def cigar_str(self):
 		return "".join(map(lambda x: str(x[0]) + 'MIDNSH'[x[1]], self._cigar))
@@ -125,7 +130,7 @@ cdef class Aligner:
 	def __bool__(self):
 		return (self._idx != NULL)
 
-	def map(self, seq, buf=None):
+	def map(self, seq, seq2=None, buf=None):
 		cdef cmappy.mm_reg1_t *regs
 		cdef cmappy.mm_hitpy_t h
 		cdef ThreadBuffer b
@@ -134,7 +139,8 @@ cdef class Aligner:
 		if self._idx is NULL: return None
 		if buf is None: b = ThreadBuffer()
 		else: b = buf
-		regs = cmappy.mm_map(self._idx, len(seq), str.encode(seq), &n_regs, b._b, &self.map_opt, NULL)
+		if seq2 is None: regs = cmappy.mm_map_aux(self._idx, str.encode(seq), NULL, &n_regs, b._b, &self.map_opt)
+		else: regs = cmappy.mm_map_aux(self._idx, str.encode(seq), str.encode(seq2), &n_regs, b._b, &self.map_opt)
 
 		for i in range(n_regs):
 			cmappy.mm_reg2hitpy(self._idx, &regs[i], &h)
@@ -142,7 +148,7 @@ cdef class Aligner:
 			for k in range(h.n_cigar32):
 				c = h.cigar32[k]
 				cigar.append([c>>4, c&0xf])
-			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)
+			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)
 			cmappy.mm_free_reg1(&regs[i])
 		free(regs)