fork/minimap2
1
0
Fork 0
mirror of https://github.com/lh3/minimap2.git synced 2026-09-16 05:27:55 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
fb8a1b55367d554f65b50e73d496b04183ee10bc
minimap2/align.c
Heng Li 79b0caca95 r537: model the next base to GT/AG 
[PMID:18688272] shows that the base following GT tends to be A or G (i.e. R) in
both human and yeast, and that the base preceeding AG tends to be C or T (i.e.
Y). In the new model, we pay no cost to GTr..yAG, but we pay half of the cost
if there is no r or y. This improves the junction accuracy when mapping to
human and mouse and decreases the accuacy when mapping to SIRV. My guess is
that SIRV does not honor this trend. Need to investigate in future.

Also in this commit, --cost-non-gt-ag is aliased to -C. The default is changed
to 9 instead of 5. I also added --splice-flank to enable the above model. This
may become the default once I confirm my hypothesis on SIRV.
2017-10-28 00:25:01 -04:00

672 lines
24 KiB
C
Raw Blame History

#include <assert.h>
#include <string.h>
#include <stdlib.h>
#include "minimap.h"
#include "mmpriv.h"
#include "ksw2.h"
static void ksw_gen_simple_mat(int m, int8_t *mat, int8_t a, int8_t b)
{
int i, j;
a = a < 0? -a : a;
b = b > 0? -b : b;
for (i = 0; i < m - 1; ++i) {
for (j = 0; j < m - 1; ++j)
mat[i * m + j] = i == j? a : b;
mat[i * m + m - 1] = 0;
}
for (j = 0; j < m; ++j)
mat[(m - 1) * m + j] = 0;
}
static inline void mm_seq_rev(uint32_t len, uint8_t *seq)
{
uint32_t i;
uint8_t t;
for (i = 0; i < len>>1; ++i)
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)
{
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;
} 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)
{
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) {
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;
}
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;
}
}
return 0;
}
static void mm_fix_cigar(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *tseq, int *qshift, int *tshift)
{
mm_extra_t *p = r->p;
int32_t k, toff = 0, qoff = 0, to_shrink = 0;
*qshift = *tshift = 0;
if (p->n_cigar <= 1) return;
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;
if (op == 0) {
toff += len, qoff += len;
} else if (op == 1 || op == 2) { // insertion or deletion
if (k > 0 && k < p->n_cigar - 1 && (p->cigar[k-1]&0xf) == 0 && (p->cigar[k+1]&0xf) == 0) {
int l, prev_len = p->cigar[k-1] >> 4;
if (op == 1) {
for (l = 0; l < prev_len; ++l)
if (qseq[qoff - 1 - l] != qseq[qoff + len - 1 - l])
break;
} else {
for (l = 0; l < prev_len; ++l)
if (tseq[toff - 1 - l] != tseq[toff + len - 1 - l])
break;
}
if (l > 0)
p->cigar[k-1] -= l<<4, p->cigar[k+1] += l<<4, qoff -= l, toff -= l;
if (l == prev_len) to_shrink = 1;
}
if (op == 1) qoff += len;
else toff += len;
} else if (op == 3) {
toff += len;
}
}
assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
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
if (p->cigar[k]>>4 != 0)
p->cigar[l++] = p->cigar[k];
p->n_cigar = l;
for (k = l = 0; k < p->n_cigar; ++k) // merge two adjacent operations if they are the same
if (k == p->n_cigar - 1 || (p->cigar[k]&0xf) != (p->cigar[k+1]&0xf))
p->cigar[l++] = p->cigar[k];
else p->cigar[k+1] += p->cigar[k]>>4<<4; // add length to the next CIGAR operator
p->n_cigar = l;
}
if ((p->cigar[0]&0xf) == 1 || (p->cigar[0]&0xf) == 2) { // get rid of leading I or D
int32_t l = p->cigar[0] >> 4;
if ((p->cigar[0]&0xf) == 1) r->qs += l, *qshift = l;
else r->rs += l, *tshift = l;
--p->n_cigar;
memmove(p->cigar, p->cigar + 1, p->n_cigar * 4);
}
}
static void mm_update_extra(mm_reg1_t *r, const uint8_t *qseq, const uint8_t *qual, const uint8_t *tseq, const int8_t *mat, int8_t q, int8_t e)
{
uint32_t k, l, toff = 0, qoff = 0;
int32_t s = 0, max = 0, qshift, tshift;
mm_extra_t *p = r->p;
if (p == 0) return;
mm_fix_cigar(r, qseq, tseq, &qshift, &tshift);
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) {
uint32_t op = p->cigar[k]&0xf, len = p->cigar[k]>>4;
if (op == 0) { // match/mismatch
int n_ambi = 0, n_diff = 0;
float n_diff2 = 0.0f;
for (l = 0; l < len; ++l) {
int cq = qseq[qoff + l], ct = tseq[toff + l];
if (ct > 3 || cq > 3) ++n_ambi;
else if (ct != cq) {
++n_diff;
n_diff2 += qual == 0 || qual[qoff + l] >= 20? 1.0f : .05f * qual[qoff + l];
}
s += mat[ct * 5 + cq];
if (s < 0) s = 0;
else max = max > s? max : s;
}
r->blen += len - n_ambi, r->mlen += len - (n_ambi + n_diff), p->n_ambi += n_ambi;
p->n_diff2 += n_diff2, p->blen2 += len - n_ambi;
toff += len, qoff += len;
} else if (op == 1) { // insertion
int n_ambi = 0;
for (l = 0; l < len; ++l)
if (qseq[qoff + l] > 3) ++n_ambi;
r->blen += len - n_ambi, p->n_ambi += n_ambi;
p->n_diff2 += 1.0f, ++p->blen2;
s -= q + e * len;
if (s < 0) s = 0;
qoff += len;
} else if (op == 2) { // deletion
int n_ambi = 0;
for (l = 0; l < len; ++l)
if (tseq[toff + l] > 3) ++n_ambi;
r->blen += len - n_ambi, p->n_ambi += n_ambi;
p->n_diff2 += 1.0f, ++p->blen2;
s -= q + e * len;
if (s < 0) s = 0;
toff += len;
} else if (op == 3) { // intron
toff += len;
}
}
p->dp_max = max;
assert(qoff == r->qe - r->qs && toff == r->re - r->rs);
}
static void mm_append_cigar(mm_reg1_t *r, uint32_t n_cigar, uint32_t *cigar) // TODO: this calls the libc realloc()
{
mm_extra_t *p;
if (n_cigar == 0) return;
if (r->p == 0) {
uint32_t capacity = n_cigar + sizeof(mm_extra_t);
kroundup32(capacity);
r->p = (mm_extra_t*)calloc(capacity, 4);
r->p->capacity = capacity;
} else if (r->p->n_cigar + n_cigar + sizeof(mm_extra_t) > r->p->capacity) {
r->p->capacity = r->p->n_cigar + n_cigar + sizeof(mm_extra_t);
kroundup32(r->p->capacity);
r->p = (mm_extra_t*)realloc(r->p, r->p->capacity * 4);
}
p = r->p;
if (p->n_cigar > 0 && (p->cigar[p->n_cigar-1]&0xf) == (cigar[0]&0xf)) { // same CIGAR op at the boundary
p->cigar[p->n_cigar-1] += cigar[0]>>4<<4;
if (n_cigar > 1) memcpy(p->cigar + p->n_cigar, cigar + 1, (n_cigar - 1) * 4);
p->n_cigar += n_cigar - 1;
} else {
memcpy(p->cigar + p->n_cigar, cigar, n_cigar * 4);
p->n_cigar += n_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)
{
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);
for (i = 0; i < tlen; ++i) fputc("ACGTN"[tseq[i]], stderr);
fputc('\n', stderr);
for (i = 0; i < qlen; ++i) fputc("ACGTN"[qseq[i]], stderr);
fputc('\n', stderr);
}
if (opt->flag & MM_F_SPLICE)
ksw_exts2_sse(km, qlen, qseq, tlen, tseq, 5, mat, opt->q, opt->e, opt->q2, opt->noncan, zdrop, flag, ez);
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);
}
static inline int mm_get_hplen_back(const mm_idx_t *mi, uint32_t rid, uint32_t x)
{
int64_t i, off0 = mi->seq[rid].offset, off = off0 + x;
int c = mm_seq4_get(mi->S, off);
for (i = off - 1; i >= off0; --i)
if (mm_seq4_get(mi->S, i) != c) break;
return (int)(off - i);
}
static inline void mm_adjust_minier(const mm_idx_t *mi, uint8_t *const qseq0[2], mm128_t *a, int32_t *r, int32_t *q)
{
if (mi->is_hpc) {
const uint8_t *qseq = qseq0[a->x>>63];
int i, c;
*q = (int32_t)a->y;
for (i = *q - 1, c = qseq[*q]; i > 0; --i)
if (qseq[i] != c) break;
*q = i + 1;
c = mm_get_hplen_back(mi, a->x<<1>>33, (int32_t)a->x);
*r = (int32_t)a->x + 1 - c;
} else {
*r = (int32_t)a->x - (mi->k>>1);
*q = (int32_t)a->y - (mi->k>>1);
}
}
static void mm_filter_bad_seeds(void *km, int as1, int cnt1, mm128_t *a, int min_gap, int diff_thres, int max_ext_len, int max_ext_cnt)
{
int max_st, max_en, n, i, k, max, *K;
for (i = 1, n = 0; i < cnt1; ++i) { // count the number of gaps longer than min_gap
int gap = ((int32_t)a[as1 + i].y - a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - a[as1 + i - 1].x);
if (gap < -min_gap || gap > min_gap) ++n;
}
if (n <= 1) return;
K = (int*)kmalloc(km, n * sizeof(int));
for (i = 1, n = 0; i < cnt1; ++i) { // store the positions of long gaps
int gap = ((int32_t)a[as1 + i].y - a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - a[as1 + i - 1].x);
if (gap < -min_gap || gap > min_gap)
K[n++] = i;
}
max = 0, max_st = max_en = -1;
for (k = 0;; ++k) { // traverse long gaps
int gap, l, n_ins = 0, n_del = 0, qs, rs, max_diff = 0, max_diff_l = -1;
if (k == n || k >= max_en) {
if (max_en > 0)
for (i = K[max_st]; i < K[max_en]; ++i)
a[as1 + i].y |= MM_SEED_IGNORE;
max = 0, max_st = max_en = -1;
if (k == n) break;
}
i = K[k];
gap = ((int32_t)a[as1 + i].y - a[as1 + i - 1].y) - ((int32_t)a[as1 + i].x - a[as1 + i - 1].x);
if (gap > 0) n_ins += gap;
else n_del += -gap;
qs = (int32_t)a[as1 + i - 1].y;
rs = (int32_t)a[as1 + i - 1].x;
for (l = k + 1; l < n && l <= k + max_ext_cnt; ++l) {
int j = K[l], diff;
if ((int32_t)a[as1 + j].y - qs > max_ext_len || (int32_t)a[as1 + j].x - rs > max_ext_len) break;
gap = ((int32_t)a[as1 + j].y - (int32_t)a[as1 + j - 1].y) - (a[as1 + j].x - a[as1 + j - 1].x);
if (gap > 0) n_ins += gap;
else n_del += -gap;
diff = n_ins + n_del - abs(n_ins - n_del);
if (max_diff < diff)
max_diff = diff, max_diff_l = l;
}
if (max_diff > diff_thres && max_diff > max)
max = max_diff, max_st = k, max_en = max_diff_l;
}
kfree(km, K);
}
static void mm_fix_bad_ends(const mm_reg1_t *r, const mm128_t *a, int bw, int32_t *as, int32_t *cnt)
{
int32_t i, l;
*as = r->as, *cnt = r->cnt;
if (r->cnt < 3) return;
l = a[r->as].y >> 32 & 0xff;
for (i = r->as + 1; i < r->as + r->cnt - 1; ++i) {
int32_t lq, lr, min, max;
lr = (int32_t)a[i].x - (int32_t)a[i-1].x;
lq = (int32_t)a[i].y - (int32_t)a[i-1].y;
min = lr < lq? lr : lq;
max = lr > lq? lr : lq;
if (max - min > l >> 1) *as = i;
l += min;
if (l >= bw << 1) break;
}
*cnt = r->as + r->cnt - *as;
l = a[r->as + r->cnt - 1].y >> 32 & 0xff;
for (i = r->as + r->cnt - 2; i > *as; --i) {
int32_t lq, lr, min, max;
lr = (int32_t)a[i+1].x - (int32_t)a[i].x;
lq = (int32_t)a[i+1].y - (int32_t)a[i].y;
min = lr < lq? lr : lq;
max = lr > lq? lr : lq;
if (max - min > l >> 1) *cnt = i + 1 - *as;
l += min;
if (l >= bw) break;
}
}
static void mm_max_stretch(const mm_mapopt_t *opt, const mm_reg1_t *r, const mm128_t *a, int32_t *as, int32_t *cnt)
{
int32_t i, score, max_score, len, max_i, max_len;
*as = r->as, *cnt = r->cnt;
if (r->cnt < 2) return;
max_score = -1, max_i = -1, max_len = 0;
score = a[r->as].y >> 32 & 0xff, len = 1;
for (i = r->as + 1; i < r->as + r->cnt; ++i) {
int32_t lq, lr, q_span;
q_span = a[i].y >> 32 & 0xff;
lr = (int32_t)a[i].x - (int32_t)a[i-1].x;
lq = (int32_t)a[i].y - (int32_t)a[i-1].y;
if (lq == lr) {
score += lq < q_span? lq : q_span;
++len;
} else {
if (score > max_score)
max_score = score, max_len = len, max_i = i - len;
score = q_span, len = 1;
}
}
if (score > max_score)
max_score = score, max_len = len, max_i = i - len;
*as = max_i, *cnt = max_len;
}
static void mm_align1(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], uint8_t *qual0[2], mm_reg1_t *r, mm_reg1_t *r2, int n_a, mm128_t *a, ksw_extz_t *ez, int splice_flag)
{
int is_sr = !!(opt->flag & MM_F_SR), is_splice = !!(opt->flag & MM_F_SPLICE);
int32_t rid = a[r->as].x<<1>>33, rev = a[r->as].x>>63, as1, cnt1;
uint8_t *tseq, *qseq;
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];
if (is_sr) assert(!mi->is_hpc); // HPC won't work with SR because with HPC we can't easily tell if there is a gap
r2->cnt = 0;
if (r->cnt == 0) return;
ksw_gen_simple_mat(5, mat, opt->a, opt->b);
bw = (int)(opt->bw * 1.5 + 1.);
if (is_sr && !mi->is_hpc) {
mm_max_stretch(opt, r, a, &as1, &cnt1);
rs = (int32_t)a[as1].x + 1 - (int32_t)(a[as1].y>>32&0xff);
qs = (int32_t)a[as1].y + 1 - (int32_t)(a[as1].y>>32&0xff);
re = (int32_t)a[as1+cnt1-1].x + 1;
qe = (int32_t)a[as1+cnt1-1].y + 1;
} else {
if (!is_splice)
mm_fix_bad_ends(r, a, opt->bw, &as1, &cnt1);
else as1 = r->as, cnt1 = r->cnt;
mm_filter_bad_seeds(km, as1, cnt1, a, 10, 40, opt->max_gap>>1, 10);
mm_adjust_minier(mi, qseq0, &a[as1], &rs, &qs);
mm_adjust_minier(mi, qseq0, &a[as1 + cnt1 - 1], &re, &qe);
}
assert(cnt1 > 0);
if (is_splice) {
if (splice_flag & MM_F_SPLICE_FOR) extra_flag |= rev? KSW_EZ_SPLICE_REV : KSW_EZ_SPLICE_FOR;
if (splice_flag & MM_F_SPLICE_REV) extra_flag |= rev? KSW_EZ_SPLICE_FOR : KSW_EZ_SPLICE_REV;
if (opt->flag & MM_F_SPLICE_FLANK) extra_flag |= KSW_EZ_SPLICE_FLANK;
}
/* Look for the start and end of regions to perform DP. This sounds easy
* but is in fact tricky. Excessively small regions lead to unnecessary
* clippings and lose alignable sequences. Excessively large regions
* occasionally lead to large overlaps between two chains and may cause
* loss of alignments in corner cases. */
if (is_sr) {
qs0 = 0, qe0 = qlen;
l = qs;
l += l * opt->a + opt->end_bonus > opt->q? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0;
rs0 = rs - l > 0? rs - l : 0;
l = qlen - qe;
l += l * opt->a + opt->end_bonus > opt->q? (l * opt->a + opt->end_bonus - opt->q) / opt->e : 0;
re0 = re + l < mi->seq[rid].len? re + l : mi->seq[rid].len;
} else {
// compute rs0 and qs0
rs0 = (int32_t)a[r->as].x + 1 - (int32_t)(a[r->as].y>>32&0xff);
qs0 = (int32_t)a[r->as].y + 1 - (int32_t)(a[r->as].y>>32&0xff);
if (rs0 < 0) rs0 = 0; // this may happen when HPC is in use
assert(qs0 >= 0); // this should never happen, or it is logic error
rs1 = qs1 = 0;
for (i = r->as - 1, l = 0; i >= 0 && a[i].x>>32 == a[r->as].x>>32; --i) { // inspect nearby seeds
int32_t x = (int32_t)a[i].x + 1 - (int32_t)(a[i].y>>32&0xff);
int32_t y = (int32_t)a[i].y + 1 - (int32_t)(a[i].y>>32&0xff);
if (x < rs0 && y < qs0) {
if (++l > opt->min_cnt) {
l = rs0 - x > qs0 - y? rs0 - x : qs0 - y;
rs1 = rs0 - l, qs1 = qs0 - l;
break;
}
}
}
if (qs > 0 && rs > 0) {
l = qs < opt->max_gap? qs : opt->max_gap;
qs1 = qs1 > qs - l? qs1 : qs - l;
qs0 = qs0 < qs1? qs0 : qs1; // at least include qs0
l += l * opt->a > opt->q? (l * opt->a - opt->q) / opt->e : 0;
l = l < opt->max_gap? l : opt->max_gap;
l = l < rs? l : rs;
rs1 = rs1 > rs - l? rs1 : rs - l;
rs0 = rs0 < rs1? rs0 : rs1;
} else rs0 = rs, qs0 = qs;
// compute re0 and qe0
re0 = (int32_t)a[r->as + r->cnt - 1].x + 1;
qe0 = (int32_t)a[r->as + r->cnt - 1].y + 1;
re1 = mi->seq[rid].len, qe1 = qlen;
for (i = r->as + r->cnt, l = 0; i < n_a && a[i].x>>32 == a[r->as].x>>32; ++i) { // inspect nearby seeds
int32_t x = (int32_t)a[i].x + 1;
int32_t y = (int32_t)a[i].y + 1;
if (x > re0 && y > qe0) {
if (++l > opt->min_cnt) {
l = x - re0 > y - qe0? x - re0 : y - qe0;
re1 = re0 + l, qe1 = qe0 + l;
break;
}
}
}
if (qe < qlen && re < mi->seq[rid].len) {
l = qlen - qe < opt->max_gap? qlen - qe : opt->max_gap;
qe1 = qe1 < qe + l? qe1 : qe + l;
qe0 = qe0 > qe1? qe0 : qe1; // at least include qe0
l += l * opt->a > opt->q? (l * opt->a - opt->q) / opt->e : 0;
l = l < opt->max_gap? l : opt->max_gap;
l = l < mi->seq[rid].len - re? l : mi->seq[rid].len - re;
re1 = re1 < re + l? re1 : re + l;
re0 = re0 > re1? re0 : re1;
} else re0 = re, qe0 = qe;
}
assert(re0 > rs0);
tseq = (uint8_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_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);
if (ez->n_cigar > 0) {
mm_append_cigar(r, ez->n_cigar, ez->cigar);
r->p->dp_score += ez->max;
}
rs1 = rs - (ez->reach_end? ez->mqe_t + 1 : ez->max_t + 1);
qs1 = qs - (ez->reach_end? qs - qs0 : ez->max_q + 1);
mm_seq_rev(qs - qs0, qseq);
} else rs1 = rs, qs1 = qs;
re1 = rs, qe1 = qs;
assert(qs1 >= 0 && rs1 >= 0);
for (i = is_sr? cnt1 - 1 : 1; i < cnt1; ++i) { // gap filling
if ((a[as1+i].y & (MM_SEED_IGNORE|MM_SEED_TANDEM)) && i != cnt1 - 1) continue;
if (is_sr && !mi->is_hpc) {
re = (int32_t)a[as1 + i].x + 1;
qe = (int32_t)a[as1 + i].y + 1;
} 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;
if (a[as1+i].y & MM_SEED_LONG_JOIN)
bw1 = qe - qs > re - rs? qe - qs : re - rs;
qseq = &qseq0[rev][qs];
mm_idx_getseq(mi, rid, rs, re, tseq);
if (is_sr) { // perform ungapped alignment
assert(qe - qs == re - rs);
ksw_reset_extz(ez);
for (j = 0, ez->score = 0; j < qe - qs; ++j) {
if (qseq[j] >= 4 || tseq[j] >= 4) ez->score += opt->e2;
else ez->score += qseq[j] == tseq[j]? opt->a : -opt->b;
}
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
}
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
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.
for (j = i - 1; j >= 0; --j)
if ((int32_t)a[as1 + j].x <= rs + ez->max_t)
break;
dropped = 1;
if (j < 0) j = 0;
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)
mm_split_reg(r, r2, as1 + j + 1 - r->as, qlen, a);
break;
} else r->p->dp_score += ez->score;
rs = re, qs = qe;
}
}
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);
if (ez->n_cigar > 0) {
mm_append_cigar(r, ez->n_cigar, ez->cigar);
r->p->dp_score += ez->max;
}
re1 = re + (ez->reach_end? ez->mqe_t + 1 : ez->max_t + 1);
qe1 = qe + (ez->reach_end? qe0 - qe : ez->max_q + 1);
}
assert(qe1 <= qlen);
r->rs = rs1, r->re = re1;
if (rev) r->qs = qlen - qe1, r->qe = qlen - qs1;
else r->qs = qs1, r->qe = qe1;
assert(re1 - rs1 <= re0 - rs0);
if (r->p) {
mm_idx_getseq(mi, rid, rs1, re1, tseq);
mm_update_extra(r, &qseq0[r->rev][qs1], qual0[r->rev]? &qual0[r->rev][qs1] : 0, tseq, mat, opt->q, opt->e);
if (rev && r->p->trans_strand)
r->p->trans_strand ^= 3; // flip to the read strand
}
kfree(km, tseq);
}
static int mm_align1_inv(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, uint8_t *qseq0[2], uint8_t *qual0[2], const mm_reg1_t *r1, const mm_reg1_t *r2, mm_reg1_t *r_inv, ksw_extz_t *ez)
{
int tl, ql, score, ret = 0, q_off, t_off;
uint8_t *tseq, *qseq, *qual;
int8_t mat[25];
void *qp;
memset(r_inv, 0, sizeof(mm_reg1_t));
if (!(r1->split&1) || !(r2->split&2)) return 0;
if (r1->id != r1->parent && r1->parent != MM_PARENT_TMP_PRI) return 0;
if (r2->id != r2->parent && r2->parent != MM_PARENT_TMP_PRI) return 0;
if (r1->rid != r2->rid || r1->rev != r2->rev) return 0;
ql = r2->qs - r1->qe;
tl = r2->rs - r1->re;
if (ql < opt->min_chain_score || ql > opt->max_gap) return 0;
if (tl < opt->min_chain_score || tl > opt->max_gap) return 0;
ksw_gen_simple_mat(5, mat, opt->a, opt->b);
tseq = (uint8_t*)kmalloc(km, tl);
mm_idx_getseq(mi, r1->rid, r1->re, r2->rs, tseq);
qseq = &qseq0[!r1->rev][qlen - r2->qs];
qual = qual0[!r1->rev]? &qseq0[!r1->rev][qlen - r2->qs] : 0;
mm_seq_rev(ql, qseq);
mm_seq_rev(tl, tseq);
qp = ksw_ll_qinit(km, 2, ql, qseq, 5, mat);
score = ksw_ll_i16(qp, tl, tseq, opt->q, opt->e, &q_off, &t_off);
kfree(km, qp);
mm_seq_rev(ql, qseq);
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);
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;
r_inv->id = -1;
r_inv->parent = MM_PARENT_UNSET;
r_inv->inv = 1;
r_inv->rev = !r1->rev;
r_inv->rid = r1->rid;
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;
mm_update_extra(r_inv, &qseq[q_off], qual? &qual[q_off] : 0, &tseq[t_off], mat, opt->q, opt->e);
ret = 1;
end_align1_inv:
kfree(km, tseq);
return ret;
}
static inline mm_reg1_t *mm_insert_reg(const mm_reg1_t *r, int i, int *n_regs, mm_reg1_t *regs)
{
regs = (mm_reg1_t*)realloc(regs, (*n_regs + 1) * sizeof(mm_reg1_t));
if (i + 1 != *n_regs)
memmove(&regs[i + 2], &regs[i + 1], sizeof(mm_reg1_t) * (*n_regs - i - 1));
regs[i + 1] = *r;
++*n_regs;
return regs;
}
mm_reg1_t *mm_align_skeleton(void *km, const mm_mapopt_t *opt, const mm_idx_t *mi, int qlen, const char *qstr, const char *qual, int *n_regs_, mm_reg1_t *regs, mm128_t *a)
{
extern unsigned char seq_nt4_table[256];
int32_t i, n_regs = *n_regs_, n_a;
uint8_t *qseq0[2], *qual0[2];
ksw_extz_t ez;
// encode the query sequence
qseq0[0] = (uint8_t*)kmalloc(km, qlen * 2);
qseq0[1] = qseq0[0] + qlen;
for (i = 0; i < qlen; ++i) {
qseq0[0][i] = seq_nt4_table[(uint8_t)qstr[i]];
qseq0[1][qlen - 1 - i] = qseq0[0][i] < 4? 3 - qseq0[0][i] : 4;
}
if (qual) {
qual0[0] = (uint8_t*)kmalloc(km, qlen * 2);
qual0[1] = qual0[0] + qlen;
for (i = 0; i < qlen; ++i)
qual0[0][i] = qual0[1][qlen - 1 - i] = qual[i] - 33;
} else qual0[0] = qual0[1] = 0;
// align through seed hits
n_a = mm_squeeze_a(km, n_regs, regs, a);
memset(&ez, 0, sizeof(ksw_extz_t));
for (i = 0; i < n_regs; ++i) {
mm_reg1_t r2;
if ((opt->flag&MM_F_SPLICE) && (opt->flag&MM_F_SPLICE_FOR) && (opt->flag&MM_F_SPLICE_REV)) { // then do two rounds of alignments for both strands
mm_reg1_t s[2], s2[2];
int which, trans_strand;
s[0] = s[1] = regs[i];
mm_align1(km, opt, mi, qlen, qseq0, qual0, &s[0], &s2[0], n_a, a, &ez, MM_F_SPLICE_FOR);
mm_align1(km, opt, mi, qlen, qseq0, qual0, &s[1], &s2[1], n_a, a, &ez, MM_F_SPLICE_REV);
if (s[0].p->dp_score > s[1].p->dp_score) which = 0, trans_strand = 1;
else if (s[0].p->dp_score < s[1].p->dp_score) which = 1, trans_strand = 2;
else trans_strand = 3, which = (qlen + s[0].p->dp_score) & 1; // randomly choose a strand, effectively
if (which == 0) {
regs[i] = s[0], r2 = s2[0];
free(s[1].p);
} else {
regs[i] = s[1], r2 = s2[1];
free(s[0].p);
}
regs[i].p->trans_strand = trans_strand;
} else { // one round of alignment
mm_align1(km, opt, mi, qlen, qseq0, qual0, &regs[i], &r2, n_a, a, &ez, opt->flag);
if (opt->flag&MM_F_SPLICE)
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) && !(opt->flag&MM_F_SR) && i > 0) { // don't try inversion alignment for -xsplice or -xsr
if (mm_align1_inv(km, opt, mi, qlen, qseq0, qual0, &regs[i-1], &regs[i], &r2, &ez)) {
regs = mm_insert_reg(&r2, i, &n_regs, regs);
++i; // skip the inserted INV alignment
}
}
}
*n_regs_ = n_regs;
kfree(km, qseq0[0]);
if (qual0[0]) kfree(km, qual0[0]);
kfree(km, ez.cigar);
mm_filter_regs(km, opt, n_regs_, regs);
mm_hit_sort_by_dp(km, n_regs_, regs);
return regs;
}
Reference in New Issue View Git Blame Copy Permalink