From edf736e110db9c2cb053a8c11c81d3ca696efd3e Mon Sep 17 00:00:00 2001 From: Heng Li Date: Mon, 17 Sep 2018 13:54:57 -0400 Subject: [PATCH] working on toy examples --- align.c | 13 +++++++---- index.c | 1 + ksw2_extd2_sse.c | 59 ++++++++++++++++++++++++++--------------------- python/cmappy.pxd | 5 ++++ 4 files changed, 47 insertions(+), 31 deletions(-) diff --git a/align.c b/align.c index 9bd320e..cea3618 100644 --- a/align.c +++ b/align.c @@ -88,13 +88,14 @@ static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq, return max_zdrop > opt->zdrop? 1 : 0; } -static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift) +static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift, int left_aln) { mm_extra_t *p = r->p; int32_t toff = 0, qoff = 0, to_shrink = 0; uint32_t k; *qshift = *tshift = 0; if (p->n_cigar <= 1) return; + if (!left_aln) goto end_left_aln; for (k = 0; k < p->n_cigar; ++k) { // indel left alignment uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4; if (len == 0) to_shrink = 1; @@ -123,6 +124,7 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, } } assert(qoff == r->qe - r->qs && toff == r->re - r->rs); +end_left_aln: if (to_shrink) { // squeeze out zero-length operations int32_t l = 0; for (k = 0; k < p->n_cigar; ++k) // squeeze out zero-length operations @@ -147,13 +149,13 @@ static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, } } -static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e) +static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e, int left_aln) { uint32_t k, l; int32_t s = 0, max = 0, qshift, tshift, toff = 0, qoff = 0; mm_extra_t *p = r->p; if (p == 0) return; - mm_fix_cigar(r, qseq, tseq, &qshift, &tshift); + mm_fix_cigar(r, qseq, tseq, &qshift, &tshift, left_aln); qseq += qshift, tseq += tshift; // qseq and tseq may be shifted due to the removal of leading I/D r->blen = r->mlen = 0; for (k = 0; k < p->n_cigar; ++k) { @@ -748,8 +750,9 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int assert(re1 - rs1 <= re0 - rs0); if (r->p) { + int left_aln = (opt->flag & MM_F_SPLICE) || mi->n_R == 0? 1 : 0; mm_idx_getseq2(mi, rid, rs1, re1, tseq, ob); - mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e); + mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e, left_aln); if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r, &qseq0[r->rev][qs1], tseq); if (rev && r->p->trans_strand) r->p->trans_strand ^= 3; // flip to the read strand @@ -809,7 +812,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i } r_inv->rs = r1->re + t_off; r_inv->re = r_inv->rs + ez->max_t + 1; - mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e); + mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, 0); if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r_inv, &qseq[q_off], &tseq[t_off]); ret = 1; end_align1_inv: diff --git a/index.c b/index.c index b47901d..6a3160a 100644 --- a/index.c +++ b/index.c @@ -143,6 +143,7 @@ int mm_idx_getseq2(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, u en1 = mi->seq[rid].offset + en; for (i = st1; i < en1; ++i) seq[i - st1] = mm_seq4_get(mi->S, i); + if (b) memset(b, 0, en - st); if (b && mi->R) { uint32_t i, z; memset(b, 0, en - st); diff --git a/ksw2_extd2_sse.c b/ksw2_extd2_sse.c index 7c2e0d6..bbde0b0 100644 --- a/ksw2_extd2_sse.c +++ b/ksw2_extd2_sse.c @@ -28,6 +28,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin #endif // ~KSW_CPU_DISPATCH { #define __dp_code_block1 \ + dt = _mm_load_si128(&dv[t]); \ z = _mm_load_si128(&s[t]); \ xt1 = _mm_load_si128(&x[t]); /* xt1 <- x[r-1][t..t+15] */ \ tmp = _mm_srli_si128(xt1, 15); /* tmp <- x[r-1][t+15] */ \ @@ -50,10 +51,10 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin #define __dp_code_block2 \ _mm_store_si128(&u[t], _mm_sub_epi8(z, vt1)); /* u[r][t..t+15] <- z - v[r-1][t-1..t+14] */ \ _mm_store_si128(&v[t], _mm_sub_epi8(z, ut)); /* v[r][t..t+15] <- z - u[r-1][t..t+15] */ \ - tmp = _mm_sub_epi8(z, q_); \ + tmp = _mm_sub_epi8(z, _mm_add_epi8(dt, q_)); \ a = _mm_sub_epi8(a, tmp); \ b = _mm_sub_epi8(b, tmp); \ - tmp = _mm_sub_epi8(z, q2_); \ + tmp = _mm_sub_epi8(z, _mm_add_epi8(dt, q2_)); \ a2= _mm_sub_epi8(a2, tmp); \ b2= _mm_sub_epi8(b2, tmp); @@ -61,8 +62,8 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin int with_cigar = !(flag&KSW_EZ_SCORE_ONLY), approx_max = !!(flag&KSW_EZ_APPROX_MAX); int32_t *H = 0, H0 = 0, last_H0_t = 0; uint8_t *qr, *sf, *mem, *mem2 = 0; - __m128i q_, q2_, qe_, qe2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_; - __m128i *u, *v, *x, *y, *x2, *y2, *s, *p = 0; + __m128i q_, q2_, qe_, qe2_, e_, e2_, zero_, sc_mch_, sc_mis_, m1_, sc_N_; + __m128i *u, *v, *x, *y, *x2, *y2, *s, *p = 0, *dv; ksw_reset_extz(ez); if (m <= 1 || qlen <= 0 || tlen <= 0) return; @@ -72,6 +73,8 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin zero_ = _mm_set1_epi8(0); q_ = _mm_set1_epi8(q); q2_ = _mm_set1_epi8(q2); + e_ = _mm_set1_epi8(e); + e2_ = _mm_set1_epi8(e2); qe_ = _mm_set1_epi8(q + e); qe2_ = _mm_set1_epi8(q2 + e2); sc_mch_ = _mm_set1_epi8(mat[0]); @@ -96,16 +99,20 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin ++long_thres; long_diff = long_thres * (e - e2) - (q2 - q) - e2; - mem = (uint8_t*)kcalloc(km, tlen_ * 8 + qlen_ + 1, 16); + mem = (uint8_t*)kcalloc(km, tlen_ * 9 + qlen_ + 1, 16); u = (__m128i*)(((size_t)mem + 15) >> 4 << 4); // 16-byte aligned - v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_; - s = y2 + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16; + v = u + tlen_, x = v + tlen_, y = x + tlen_, x2 = y + tlen_, y2 = x2 + tlen_, dv = y2 + tlen_; + s = dv + tlen_, sf = (uint8_t*)(s + tlen_), qr = sf + tlen_ * 16; memset(u, -q - e, tlen_ * 16); memset(v, -q - e, tlen_ * 16); memset(x, -q - e, tlen_ * 16); memset(y, -q - e, tlen_ * 16); memset(x2, -q2 - e2, tlen_ * 16); memset(y2, -q2 - e2, tlen_ * 16); + if (qd) { + int8_t *tmp = (int8_t*)dv; + for (t = 0; t < tlen; ++t) tmp[t] = -qd[t]; + } if (!approx_max) { H = (int32_t*)kmalloc(km, tlen_ * 16 * 4); for (t = 0; t < tlen_ * 16; ++t) H[t] = KSW_NEG_INF; @@ -182,7 +189,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin assert(en_ - st_ + 1 <= n_col_); if (!with_cigar) { // score only for (t = st_; t <= en_; ++t) { - __m128i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __m128i z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt; __dp_code_block1; #ifdef __SSE4_1__ z = _mm_max_epi8(z, a); @@ -191,10 +198,10 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin z = _mm_max_epi8(z, b2); z = _mm_min_epi8(z, sc_mch_); __dp_code_block2; // save u[] and v[]; update a, b, a2 and b2 - _mm_store_si128(&x[t], _mm_sub_epi8(_mm_max_epi8(a, zero_), qe_)); - _mm_store_si128(&y[t], _mm_sub_epi8(_mm_max_epi8(b, zero_), qe_)); - _mm_store_si128(&x2[t], _mm_sub_epi8(_mm_max_epi8(a2, zero_), qe2_)); - _mm_store_si128(&y2[t], _mm_sub_epi8(_mm_max_epi8(b2, zero_), qe2_)); + _mm_store_si128(&x[t], _mm_sub_epi8(_mm_max_epi8(a, zero_), _mm_add_epi8(dt, qe_))); + _mm_store_si128(&y[t], _mm_sub_epi8(_mm_max_epi8(b, zero_), _mm_add_epi8(dt, qe_))); + _mm_store_si128(&x2[t], _mm_sub_epi8(_mm_max_epi8(a2, zero_), _mm_add_epi8(dt, qe2_))); + _mm_store_si128(&y2[t], _mm_sub_epi8(_mm_max_epi8(b2, zero_), _mm_add_epi8(dt, qe2_))); #else tmp = _mm_cmpgt_epi8(a, z); z = _mm_or_si128(_mm_andnot_si128(tmp, z), _mm_and_si128(tmp, a)); @@ -208,20 +215,20 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin z = _mm_or_si128(_mm_and_si128(tmp, sc_mch_), _mm_andnot_si128(tmp, z)); __dp_code_block2; tmp = _mm_cmpgt_epi8(a, zero_); - _mm_store_si128(&x[t], _mm_sub_epi8(_mm_and_si128(tmp, a), qe_)); + _mm_store_si128(&x[t], _mm_sub_epi8(_mm_and_si128(tmp, a), _mm_add_epi8(dt, qe_))); tmp = _mm_cmpgt_epi8(b, zero_); - _mm_store_si128(&y[t], _mm_sub_epi8(_mm_and_si128(tmp, b), qe_)); + _mm_store_si128(&y[t], _mm_sub_epi8(_mm_and_si128(tmp, b), _mm_add_epi8(dt, qe_))); tmp = _mm_cmpgt_epi8(a2, zero_); - _mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), qe2_)); + _mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), _mm_add_epi8(dt, qe2_))); tmp = _mm_cmpgt_epi8(b2, zero_); - _mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), qe2_)); + _mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), _mm_add_epi8(dt, qe2_))); #endif } } else if (!(flag&KSW_EZ_RIGHT)) { // gap left-alignment __m128i *pr = p + (size_t)r * n_col_ - st_; off[r] = st, off_end[r] = en; for (t = st_; t <= en_; ++t) { - __m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt; __dp_code_block1; #ifdef __SSE4_1__ d = _mm_and_si128(_mm_cmpgt_epi8(a, z), _mm_set1_epi8(1)); // d = a > z? 1 : 0 @@ -251,16 +258,16 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin #endif __dp_code_block2; tmp = _mm_cmpgt_epi8(a, zero_); - _mm_store_si128(&x[t], _mm_sub_epi8(_mm_and_si128(tmp, a), qe_)); + _mm_store_si128(&x[t], _mm_sub_epi8(_mm_and_si128(tmp, a), _mm_add_epi8(dt, qe_))); d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x08))); // d = a > 0? 1<<3 : 0 tmp = _mm_cmpgt_epi8(b, zero_); - _mm_store_si128(&y[t], _mm_sub_epi8(_mm_and_si128(tmp, b), qe_)); + _mm_store_si128(&y[t], _mm_sub_epi8(_mm_and_si128(tmp, b), _mm_add_epi8(dt, qe_))); d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x10))); // d = b > 0? 1<<4 : 0 tmp = _mm_cmpgt_epi8(a2, zero_); - _mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), qe2_)); + _mm_store_si128(&x2[t], _mm_sub_epi8(_mm_and_si128(tmp, a2), _mm_add_epi8(dt, qe2_))); d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x20))); // d = a > 0? 1<<5 : 0 tmp = _mm_cmpgt_epi8(b2, zero_); - _mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), qe2_)); + _mm_store_si128(&y2[t], _mm_sub_epi8(_mm_and_si128(tmp, b2), _mm_add_epi8(dt, qe2_))); d = _mm_or_si128(d, _mm_and_si128(tmp, _mm_set1_epi8(0x40))); // d = b > 0? 1<<6 : 0 _mm_store_si128(&pr[t], d); } @@ -268,7 +275,7 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin __m128i *pr = p + (size_t)r * n_col_ - st_; off[r] = st, off_end[r] = en; for (t = st_; t <= en_; ++t) { - __m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp; + __m128i d, z, a, b, a2, b2, xt1, x2t1, vt1, ut, tmp, dt; __dp_code_block1; #ifdef __SSE4_1__ d = _mm_andnot_si128(_mm_cmpgt_epi8(z, a), _mm_set1_epi8(1)); // d = z > a? 0 : 1 @@ -298,16 +305,16 @@ void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin #endif __dp_code_block2; tmp = _mm_cmpgt_epi8(zero_, a); - _mm_store_si128(&x[t], _mm_sub_epi8(_mm_andnot_si128(tmp, a), qe_)); + _mm_store_si128(&x[t], _mm_sub_epi8(_mm_andnot_si128(tmp, a), _mm_add_epi8(dt, qe_))); d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x08))); // d = a > 0? 1<<3 : 0 tmp = _mm_cmpgt_epi8(zero_, b); - _mm_store_si128(&y[t], _mm_sub_epi8(_mm_andnot_si128(tmp, b), qe_)); + _mm_store_si128(&y[t], _mm_sub_epi8(_mm_andnot_si128(tmp, b), _mm_add_epi8(dt, qe_))); d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x10))); // d = b > 0? 1<<4 : 0 tmp = _mm_cmpgt_epi8(zero_, a2); - _mm_store_si128(&x2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, a2), qe2_)); + _mm_store_si128(&x2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, a2), _mm_add_epi8(dt, qe2_))); d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x20))); // d = a > 0? 1<<5 : 0 tmp = _mm_cmpgt_epi8(zero_, b2); - _mm_store_si128(&y2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, b2), qe2_)); + _mm_store_si128(&y2[t], _mm_sub_epi8(_mm_andnot_si128(tmp, b2), _mm_add_epi8(dt, qe2_))); d = _mm_or_si128(d, _mm_andnot_si128(tmp, _mm_set1_epi8(0x40))); // d = b > 0? 1<<6 : 0 _mm_store_si128(&pr[t], d); } diff --git a/python/cmappy.pxd b/python/cmappy.pxd index a6c2a3d..28ce3f7 100644 --- a/python/cmappy.pxd +++ b/python/cmappy.pxd @@ -53,6 +53,10 @@ cdef extern from "minimap.h": uint64_t offset uint32_t len + ctypedef struct mm_idx_bed_t: + uint64_t x + int32_t end, idx + ctypedef struct mm_idx_bucket_t: pass @@ -62,6 +66,7 @@ cdef extern from "minimap.h": mm_idx_seq_t *seq uint32_t *S mm_idx_bucket_t *B + mm_idx_bed_t *R void *km void *h