From d01369ec47ff61e5fe3d3477fdf378b6bb42931d Mon Sep 17 00:00:00 2001 From: Heng Li Date: Mon, 17 Sep 2018 12:31:42 -0400 Subject: [PATCH] pass BED regions to extd2; not tested --- align.c | 26 ++++++++++++++------------ index.c | 17 ++++++++++++++++- ksw2.h | 2 +- ksw2_dispatch.c | 10 +++++----- ksw2_extd2_sse.c | 6 +++--- mmpriv.h | 1 + 6 files changed, 40 insertions(+), 22 deletions(-) diff --git a/align.c b/align.c index cd28d6a..9bd320e 100644 --- a/align.c +++ b/align.c @@ -290,7 +290,7 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t r->p = p; } -static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez) +static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint8_t *qseq, int tlen, const uint8_t *tseq, const int8_t *ob, const int8_t *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez) { if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) { int i; @@ -305,7 +305,7 @@ static void mm_align_pair(void *km, const mm_mapopt_t *opt, int qlen, const uint else if (opt->q == opt->q2 && opt->e == opt->e2) ksw_extz2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, w, zdrop, end_bonus, flag, ez); else - ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ez); + ksw_extd2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->e2, w, zdrop, end_bonus, flag, ob, ez); if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) { int i; fprintf(stderr, "score=%d, cigar=", ez->score); @@ -547,7 +547,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int int32_t i, l, bw, dropped = 0, extra_flag = 0, rs0, re0, qs0, qe0; int32_t rs, re, qs, qe; int32_t rs1, qs1, re1, qe1; - int8_t mat[25]; + int8_t mat[25], *ob; if (is_sr) assert(!(mi->flag & MM_I_HPC)); // HPC won't work with SR because with HPC we can't easily tell if there is a gap @@ -660,13 +660,14 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int assert(re0 > rs0); tseq = (uint8_t*)kmalloc(km, re0 - rs0); + ob = (int8_t*)kmalloc(km, re0 - rs0); if (qs > 0 && rs > 0) { // left extension qseq = &qseq0[rev][qs0]; - mm_idx_getseq(mi, rid, rs0, rs, tseq); + mm_idx_getseq2(mi, rid, rs0, rs, tseq, ob); mm_seq_rev(qs - qs0, qseq); mm_seq_rev(rs - rs0, tseq); - mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez); + mm_align_pair(km, opt, qs - qs0, qseq, rs - rs0, tseq, ob, mat, bw, opt->end_bonus, r->split_inv? opt->zdrop_inv : opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY|KSW_EZ_RIGHT|KSW_EZ_REV_CIGAR, ez); if (ez->n_cigar > 0) { mm_append_cigar(r, ez->n_cigar, ez->cigar); r->p->dp_score += ez->max; @@ -691,7 +692,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int bw1 = qe - qs > re - rs? qe - qs : re - rs; // perform alignment qseq = &qseq0[rev][qs]; - mm_idx_getseq(mi, rid, rs, re, tseq); + mm_idx_getseq2(mi, rid, rs, re, tseq, ob); if (is_sr) { // perform ungapped alignment assert(qe - qs == re - rs); ksw_reset_extz(ez); @@ -701,11 +702,11 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int } ez->cigar = ksw_push_cigar(km, &ez->n_cigar, &ez->m_cigar, ez->cigar, 0, qe - qs); } else { // perform normal gapped alignment - mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop + mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, ob, mat, bw1, -1, opt->zdrop, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop } // test Z-drop and inversion Z-drop if ((zdrop_code = mm_test_zdrop(km, opt, qseq, tseq, ez->n_cigar, ez->cigar, mat)) != 0) - mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate + mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, ob, mat, bw1, -1, zdrop_code == 2? opt->zdrop_inv : opt->zdrop, extra_flag, ez); // second pass: lift approximate // update CIGAR if (ez->n_cigar > 0) mm_append_cigar(r, ez->n_cigar, ez->cigar); @@ -730,8 +731,8 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int if (!dropped && qe < qe0 && re < re0) { // right extension qseq = &qseq0[rev][qe]; - mm_idx_getseq(mi, rid, re, re0, tseq); - mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez); + mm_idx_getseq2(mi, rid, re, re0, tseq, ob); + mm_align_pair(km, opt, qe0 - qe, qseq, re0 - re, tseq, ob, mat, bw, opt->end_bonus, opt->zdrop, extra_flag|KSW_EZ_EXTZ_ONLY, ez); if (ez->n_cigar > 0) { mm_append_cigar(r, ez->n_cigar, ez->cigar); r->p->dp_score += ez->max; @@ -747,7 +748,7 @@ 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) { - mm_idx_getseq(mi, rid, rs1, re1, tseq); + mm_idx_getseq2(mi, rid, rs1, re1, tseq, ob); mm_update_extra(r, &qseq0[r->rev][qs1], tseq, mat, opt->q, opt->e); if (opt->flag & MM_F_EQX) mm_update_cigar_eqx(r, &qseq0[r->rev][qs1], tseq); if (rev && r->p->trans_strand) @@ -755,6 +756,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int } kfree(km, tseq); + kfree(km, ob); } static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez) @@ -788,7 +790,7 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i mm_seq_rev(tl, tseq); if (score < opt->min_dp_max) goto end_align1_inv; q_off = ql - (q_off + 1), t_off = tl - (t_off + 1); - mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, mat, (int)(opt->bw * 1.5), -1, opt->zdrop, KSW_EZ_EXTZ_ONLY, ez); + mm_align_pair(km, opt, ql - q_off, qseq + q_off, tl - t_off, tseq + t_off, 0, mat, (int)(opt->bw * 1.5), -1, opt->zdrop, KSW_EZ_EXTZ_ONLY, ez); if (ez->n_cigar == 0) goto end_align1_inv; // should never be here mm_append_cigar(r_inv, ez->n_cigar, ez->cigar); r_inv->p->dp_score = ez->max; diff --git a/index.c b/index.c index 5be102c..b47901d 100644 --- a/index.c +++ b/index.c @@ -134,7 +134,7 @@ int mm_idx_name2id(const mm_idx_t *mi, const char *name) return k == kh_end(h)? -1 : kh_val(h, k); } -int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq) +int mm_idx_getseq2(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq, int8_t *b) { uint64_t i, st1, en1; if (rid >= mi->n_seq || st >= mi->seq[rid].len) return -1; @@ -143,9 +143,24 @@ int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, ui en1 = mi->seq[rid].offset + en; for (i = st1; i < en1; ++i) seq[i - st1] = mm_seq4_get(mi->S, i); + if (b && mi->R) { + uint32_t i, z; + memset(b, 0, en - st); + z = mm_idx_bed_query(mi, (uint64_t)rid << 32 | st); + for (i = z + 1; i < mi->n_R; ++i) { + uint32_t j, rr, rs, re; + rr = mi->R[i].x >> 32, rs = (uint32_t)mi->R[i].x, re = mi->R[i].end; + if (rr != rid || rs >= en) break; + assert(rs >= st); + for (j = rs; j < re; ++j) + b[j - st] = 1; + } + } return en - st; } +int mm_idx_getseq(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq) { return mm_idx_getseq2(mi, rid, st, en, seq, 0); } + int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f) { int i; diff --git a/ksw2.h b/ksw2.h index 213c27f..afe5325 100644 --- a/ksw2.h +++ b/ksw2.h @@ -58,7 +58,7 @@ void ksw_extd(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int flag, ksw_extz_t *ez); void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, - int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez); + int8_t gapo, int8_t gape, int8_t gapo2, int8_t gape2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez); void ksw_exts2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, int8_t gapo, int8_t gape, int8_t gapo2, int8_t noncan, int zdrop, int flag, ksw_extz_t *ez); diff --git a/ksw2_dispatch.c b/ksw2_dispatch.c index c4e3e07..0e18ffb 100644 --- a/ksw2_dispatch.c +++ b/ksw2_dispatch.c @@ -64,18 +64,18 @@ void ksw_extz2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uin } void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, - int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez) + int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez) { extern void ksw_extd2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, - int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez); + int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez); extern void ksw_extd2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, - int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez); + int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez); unsigned simd; simd = x86_simd(); if (simd & SIMD_SSE4_1) - ksw_extd2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, ez); + ksw_extd2_sse41(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, qd, ez); else if (simd & SIMD_SSE2) - ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, ez); + ksw_extd2_sse2(km, qlen, query, tlen, target, m, mat, q, e, q2, e2, w, zdrop, end_bonus, flag, qd, ez); else abort(); } diff --git a/ksw2_extd2_sse.c b/ksw2_extd2_sse.c index b578274..7c2e0d6 100644 --- a/ksw2_extd2_sse.c +++ b/ksw2_extd2_sse.c @@ -17,14 +17,14 @@ #ifdef KSW_CPU_DISPATCH #ifdef __SSE4_1__ void ksw_extd2_sse41(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, - int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez) + int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez) #else void ksw_extd2_sse2(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, - int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez) + int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez) #endif #else void ksw_extd2_sse(void *km, int qlen, const uint8_t *query, int tlen, const uint8_t *target, int8_t m, const int8_t *mat, - int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, ksw_extz_t *ez) + int8_t q, int8_t e, int8_t q2, int8_t e2, int w, int zdrop, int end_bonus, int flag, const int8_t *qd, ksw_extz_t *ez) #endif // ~KSW_CPU_DISPATCH { #define __dp_code_block1 \ diff --git a/mmpriv.h b/mmpriv.h index 9e4a95b..2e56bb0 100644 --- a/mmpriv.h +++ b/mmpriv.h @@ -74,6 +74,7 @@ void mm_write_sam2(kstring_t *s, const mm_idx_t *mi, const mm_bseq1_t *t, int se void mm_idxopt_init(mm_idxopt_t *opt); const uint64_t *mm_idx_get(const mm_idx_t *mi, uint64_t minier, int *n); int32_t mm_idx_cal_max_occ(const mm_idx_t *mi, float f); +int mm_idx_getseq2(const mm_idx_t *mi, uint32_t rid, uint32_t st, uint32_t en, uint8_t *seq, int8_t *b); mm128_t *mm_chain_dp(int max_dist_x, int max_dist_y, int bw, int max_skip, int min_cnt, int min_sc, int is_cdna, int n_segs, int64_t n, mm128_t *a, int *n_u_, uint64_t **_u, void *km); mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, int *n_regs_, mm_reg1_t *regs, mm128_t *a);