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

dev-r1079: per-read error rate

Browse Source 
more tuning needed
This commit is contained in:
Heng Li
2021-07-18 20:38:53 -04:00
parent 8a6edab847
commit 161ae7ff73
6 changed files with 85 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files
align.c
+74 -9
View File
@@ -237,7 +237,7 @@ static void mm_update_cigar_eqx(mm_reg1_t *r, const uint8_t *qseq, const uint8_t
r->p = p;
}
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 is_eqx, int b2)
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 is_eqx, int log_gap)
{
uint32_t k, l;
int32_t qshift, tshift, toff = 0, qoff = 0;
@@ -255,7 +255,7 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
int cq = qseq[qoff + l], ct = tseq[toff + l];
if (ct > 3 || cq > 3) ++n_ambi;
else if (ct != cq) ++n_diff;
s += b2 <= 0? mat[ct * 5 + cq] : (ct > 3 || cq > 3)? 0 : ct == cq? 1 : -b2;
s += mat[ct * 5 + cq];
if (s < 0) s = 0;
else max = max > s? max : s;
}
@@ -266,8 +266,8 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
for (l = 0; l < len; ++l)
if (qseq[qoff + l] > 3) ++n_ambi;
r->blen += len - n_ambi, p->n_ambi += n_ambi;
if (b2 > 0) s -= b2 + (double)mg_log2(1.0 + len);
else s -= q + e * len;
if (log_gap) s -= q + (double)e * mg_log2(1.0 + len);
else s -= q + e;
if (s < 0) s = 0;
qoff += len;
} else if (op == MM_CIGAR_DEL) {
@@ -275,8 +275,8 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
for (l = 0; l < len; ++l)
if (tseq[toff + l] > 3) ++n_ambi;
r->blen += len - n_ambi, p->n_ambi += n_ambi;
if (b2 > 0) s -= b2 + (double)mg_log2(1.0 + len);
else s -= q + e * len;
if (log_gap) s -= q + (double)e * mg_log2(1.0 + len);
else s -= q + e;
if (s < 0) s = 0;
toff += len;
} else if (op == MM_CIGAR_N_SKIP) {
@@ -284,7 +284,6 @@ static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *ts
}
}
p->dp_max = (int32_t)(max + .499);
if (b2 > 0) p->dp_max *= mat[0]; // for compatibility with mm_set_mapq()
assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
if (is_eqx) mm_update_cigar_eqx(r, qseq, tseq); // NB: it has to be called here as changes to qseq and tseq are not returned
}
@@ -817,7 +816,7 @@ static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int
mm_idx_getseq(mi, rid, rs1, re1, tseq);
qseq = &qseq0[r->rev][qs1];
}
mm_update_extra(r, qseq, tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX, opt->b2);
mm_update_extra(r, qseq, tseq, mat, opt->q, opt->e, opt->flag & MM_F_EQX, !(opt->flag & MM_F_SR));
if (rev && r->p->trans_strand)
r->p->trans_strand ^= 3; // flip to the read strand
}
@@ -876,7 +875,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, opt->flag & MM_F_EQX, opt->b2);
mm_update_extra(r_inv, &qseq[q_off], &tseq[t_off], mat, opt->q, opt->e, opt->flag & MM_F_EQX, !(opt->flag & MM_F_SR));
ret = 1;
end_align1_inv:
kfree(km, tseq);
@@ -893,6 +892,70 @@ static inline mm_reg1_t *mm_insert_reg(const mm_reg1_t *r, int i, int *n_regs, m
return regs;
}
static inline void mm_count_gaps(const mm_reg1_t *r, int32_t *n_gap_, int32_t *n_gapo_)
{
uint32_t i;
int32_t n_gapo = 0, n_gap = 0;
*n_gap_ = *n_gapo_ = -1;
if (r->p == 0) return;
for (i = 0; i < r->p->n_cigar; ++i) {
int32_t op = r->p->cigar[i] & 0xf, len = r->p->cigar[i] >> 4;
if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL)
++n_gapo, n_gap += len;
}
*n_gap_ = n_gap, *n_gapo_ = n_gapo;
}
double mm_event_identity(const mm_reg1_t *r)
{
int32_t n_gap, n_gapo;
if (r->p == 0) return -1.0f;
mm_count_gaps(r, &n_gap, &n_gapo);
return (double)r->mlen / (r->blen + r->p->n_ambi - n_gap + n_gapo);
}
static int32_t mm_recal_max_dp(const mm_reg1_t *r, double b2, int32_t match_sc)
{
uint32_t i;
int32_t n_gap = 0, n_gapo = 0, n_mis;
double gap_cost = 0.0;
if (r->p == 0) return -1;
for (i = 0; i < r->p->n_cigar; ++i) {
int32_t op = r->p->cigar[i] & 0xf, len = r->p->cigar[i] >> 4;
if (op == MM_CIGAR_INS || op == MM_CIGAR_DEL) {
gap_cost += b2 + (double)mg_log2(1.0 + len);
++n_gapo, n_gap += len;
}
}
n_mis = r->blen + r->p->n_ambi - r->mlen - n_gap;
return (int32_t)(match_sc * (r->mlen - b2 * n_mis - gap_cost) + .499);
}
static void mm_update_dp_max(int qlen, int n_regs, mm_reg1_t *regs, float frac, int a, int b)
{
int32_t max = -1, max2 = -1, i, max_i = -1;
double div, b2;
if (n_regs < 2) return;
for (i = 0; i < n_regs; ++i) {
mm_reg1_t *r = &regs[i];
if (r->p == 0) continue;
if (r->p->dp_max > max) max2 = max, max = r->p->dp_max, max_i = i;
else if (r->p->dp_max > max2) max2 = r->p->dp_max;
}
if (max_i < 0 || max < 0 || max2 < 0) return;
if (regs[max_i].qe - regs[max_i].qs < (double)qlen * frac) return;
if (max2 < (double)max * frac) return;
div = 1. - mm_event_identity(&regs[max_i]);
if (div < 0.02) div = 0.02;
b2 = 0.5 / div; // max value: 25
if (b2 * a < b) b2 = (double)a / b;
for (i = 0; i < n_regs; ++i) {
mm_reg1_t *r = &regs[i];
if (r->p == 0) continue;
r->p->dp_max = mm_recal_max_dp(r, b2, a);
}
}
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)
{
extern unsigned char seq_nt4_table[256];
@@ -947,6 +1010,8 @@ mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *m
kfree(km, qseq0[0]);
kfree(km, ez.cigar);
mm_filter_regs(opt, qlen, n_regs_, regs);
if (!(opt->flag&MM_F_SR) && qlen >= opt->rank_min_len)
mm_update_dp_max(qlen, *n_regs_, regs, opt->rank_frac, opt->a, opt->b);
mm_hit_sort(km, n_regs_, regs, opt->alt_drop);
return regs;
}