fork/minimap2
1
0
Fork 0
mirror of https://github.com/lh3/minimap2.git synced 2026-10-03 01:18:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

r708: implemented double Z-drop thresholds (#112)

Browse Source 
When aligning long reads, we would prefer to align through low-quality
regions. This requires a large Z-drop threshold. However, to find small
inversions, we need to use a small Z-drop. This commit address this
conflict with two Z-drop thresholds. When Z-drop exceeds the smaller
threshold, we perform a local alignment to check if there is a potential
inversion. If there is one, we break the alignment; otherwise we break
the alignment only if Z-drop excess the larger threshold.

This commit also fixes a bug that reported wrong coordinates when the
inversion is on the forward strand (#112).
This commit is contained in:
Heng Li
2018-02-15 10:50:49 -05:00
parent c0e0d5d84b
commit 1372977a37
6 changed files with 77 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files
align.c
+58 -28
View File
@@ -28,39 +28,63 @@ static inline void mm_seq_rev(uint32_t len, uint8_t *seq)
t = seq[i], seq[i] = seq[len - 1 - i], seq[len - 1 - i] = t;
}
static inline int test_zdrop_aux(int32_t score, int i, int j, int32_t *max, int *max_i, int *max_j, int e, int zdrop)
static inline void update_max_zdrop(int32_t score, int i, int j, int32_t *max, int *max_i, int *max_j, int e, int *max_zdrop, int pos[2][2])
{
if (score < *max) {
int li = i - *max_i;
int lj = j - *max_j;
int diff = li > lj? li - lj : lj - li;
if (*max - score > zdrop + diff * e)
return 1;
int z = *max - score - diff * e;
if (z > *max_zdrop) {
*max_zdrop = z;
pos[0][0] = *max_i, pos[0][1] = i + 1;
pos[1][0] = *max_j, pos[1][1] = j + 1;
}
} else *max = score, *max_i = i, *max_j = j;
return 0;
}
static int mm_check_zdrop(const uint8_t *qseq, const uint8_t *tseq, uint32_t n_cigar, uint32_t *cigar, const int8_t *mat, int8_t q, int8_t e, int zdrop)
static int mm_test_zdrop(void *km, const mm_mapopt_t *opt, const uint8_t *qseq, const uint8_t *tseq, uint32_t n_cigar, uint32_t *cigar, const int8_t *mat)
{
uint32_t k;
int32_t score = 0, max = 0, max_i = -1, max_j = -1, i = 0, j = 0;
for (k = 0; k < n_cigar; ++k) {
int32_t score = 0, max = 0, max_i = -1, max_j = -1, i = 0, j = 0, max_zdrop = 0;
int pos[2][2] = {{-1, -1}, {-1, -1}}, q_len, t_len;
// find the score and the region where score drops most along diagonal
for (k = 0, score = 0; k < n_cigar; ++k) {
uint32_t l, op = cigar[k]&0xf, len = cigar[k]>>4;
if (op == 0) {
for (l = 0; l < len; ++l) {
score += mat[tseq[i + l] * 5 + qseq[j + l]];
if (test_zdrop_aux(score, i+l, j+l, &max, &max_i, &max_j, e, zdrop)) return 1;
update_max_zdrop(score, i+l, j+l, &max, &max_i, &max_j, opt->e, &max_zdrop, pos);
}
i += len, j += len;
} else if (op == 1) {
score -= q + e * len, j += len;
if (test_zdrop_aux(score, i, j, &max, &max_i, &max_j, e, zdrop)) return 1;
} else if (op == 2 || op == 3) {
score -= q + e * len, i += len;
if (test_zdrop_aux(score, i, j, &max, &max_i, &max_j, e, zdrop)) return 1;
} else if (op == 1 || op == 2 || op == 3) {
score -= opt->q + opt->e * len;
if (op == 1) j += len; // insertion
else i += len; // deletion
update_max_zdrop(score, i, j, &max, &max_i, &max_j, opt->e, &max_zdrop, pos);
}
}
return 0;
// test if there is an inversion in the most dropped region
q_len = pos[1][1] - pos[1][0], t_len = pos[0][1] - pos[0][0];
if (!(opt->flag&(MM_F_SPLICE|MM_F_SR|MM_F_FOR_ONLY|MM_F_REV_ONLY)) && max_zdrop > opt->zdrop_inv && q_len < opt->max_gap && t_len < opt->max_gap) {
uint8_t *qseq2;
void *qp;
int q_off, t_off;
qseq2 = (uint8_t*)kmalloc(km, q_len);
for (i = 0; i < q_len; ++i) {
int c = qseq[pos[1][1] - i - 1];
qseq2[i] = c >= 4? 4 : 3 - c;
}
qp = ksw_ll_qinit(km, 2, q_len, qseq2, 5, mat);
score = ksw_ll_i16(qp, t_len, tseq + pos[0][0], opt->q, opt->e, &q_off, &t_off);
kfree(km, qseq2);
kfree(km, qp);
if (score >= opt->min_chain_score * opt->a && score >= opt->min_dp_max)
return 2; // there is a potential inversion
}
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)
@@ -193,9 +217,8 @@ static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) //
}
}
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 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 *mat, int w, int end_bonus, int zdrop, int flag, ksw_extz_t *ez)
{
int zdrop = opt->zdrop;
if (mm_dbg_flag & MM_DBG_PRINT_ALN_SEQ) {
int i;
fprintf(stderr, "===> q=(%d,%d), e=(%d,%d), bw=%d, flag=%d, zdrop=%d <===\n", opt->q, opt->q2, opt->e, opt->e2, w, flag, opt->zdrop);
@@ -522,7 +545,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
mm_idx_getseq(mi, rid, rs0, rs, tseq);
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, 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, 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;
@@ -542,9 +565,10 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
} else mm_adjust_minier(mi, qseq0, &a[as1 + i], &re, &qe);
re1 = re, qe1 = qe;
if (i == cnt1 - 1 || (a[as1+i].y&MM_SEED_LONG_JOIN) || (qe - qs >= opt->min_ksw_len && re - rs >= opt->min_ksw_len)) {
int j, bw1 = bw;
int j, bw1 = bw, zdrop_code;
if (a[as1+i].y & MM_SEED_LONG_JOIN)
bw1 = qe - qs > re - rs? qe - qs : re - rs;
// perform alignment
qseq = &qseq0[rev][qs];
mm_idx_getseq(mi, rid, rs, re, tseq);
if (is_sr) { // perform ungapped alignment
@@ -556,10 +580,12 @@ 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, extra_flag|KSW_EZ_APPROX_MAX, ez); // first pass: with approximate Z-drop
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
}
if (mm_check_zdrop(qseq, tseq, ez->n_cigar, ez->cigar, mat, opt->q, opt->e, opt->zdrop))
mm_align_pair(km, opt, qe - qs, qseq, re - rs, tseq, mat, bw1, -1, extra_flag, ez); // second pass: lift approximate
// 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
// update CIGAR
if (ez->n_cigar > 0)
mm_append_cigar(r, ez->n_cigar, ez->cigar);
if (ez->zdropped) { // truncated by Z-drop; TODO: sometimes Z-drop kicks in because the next seed placement is wrong. This can be fixed in principle.
@@ -571,8 +597,10 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
r->p->dp_score += ez->max;
re1 = rs + (ez->max_t + 1);
qe1 = qs + (ez->max_q + 1);
if (cnt1 - (j + 1) >= opt->min_cnt)
if (cnt1 - (j + 1) >= opt->min_cnt) {
mm_split_reg(r, r2, as1 + j + 1 - r->as, qlen, a);
if (zdrop_code == 2) r2->split_inv = 1;
}
break;
} else r->p->dp_score += ez->score;
rs = re, qs = qe;
@@ -582,7 +610,7 @@ 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, extra_flag|KSW_EZ_EXTZ_ONLY, ez);
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);
if (ez->n_cigar > 0) {
mm_append_cigar(r, ez->n_cigar, ez->cigar);
r->p->dp_score += ez->max;
@@ -638,7 +666,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, KSW_EZ_EXTZ_ONLY, ez);
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);
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;
@@ -648,8 +676,10 @@ static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, i
r_inv->rev = !r1->rev;
r_inv->rid = r1->rid;
r_inv->div = -1.0f;
r_inv->qs = r1->qe + q_off, r_inv->qe = r_inv->qs + ez->max_q + 1;
r_inv->rs = r1->re + t_off, r_inv->re = r_inv->rs + ez->max_t + 1;
r_inv->qs = r1->rev? r2->qe + q_off : r1->qe + q_off;
r_inv->qe = r_inv->qs + ez->max_q + 1;
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);
ret = 1;
end_align1_inv:
@@ -710,7 +740,7 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
regs[i].p->trans_strand = opt->flag&MM_F_SPLICE_FOR? 1 : 2;
}
if (r2.cnt > 0) regs = mm_insert_reg(&r2, i, &n_regs, regs);
if (!(opt->flag&(MM_F_SPLICE|MM_F_SR)) && !(opt->flag&(MM_F_FOR_ONLY|MM_F_REV_ONLY)) && i > 0) { // don't try inversion alignment for -xsplice or -xsr, or --for-only/rev-only
if (i > 0 && regs[i].split_inv) {
if (mm_align1_inv(km, opt, mi, qlen, qseq0, &regs[i-1], &regs[i], &r2, &ez)) {
regs = mm_insert_reg(&r2, i, &n_regs, regs);
++i; // skip the inserted INV alignment