Permalink
Find file
Fetching contributors…
Cannot retrieve contributors at this time
192 lines (177 sloc) 5.8 KB
#include <string.h>
#include <stdint.h>
#include <stdlib.h>
#include "bwtaln.h"
#include "stdaln.h"
/*
Here is a delicate example. ref_nt=ATTAAC(RBRBG), read_cs=RBBOG. If we
decode as ATTGAC(RBGOG), there are one color change and one nt change;
if we decode as ATTAAC(RBRBG), there are two color changes.
In DP, if color quality is smaller than COLOR_MM, we will use COLOR_MM
as the penalty; otherwise, we will use color quality as the
penalty. This means we always prefer two consistent color changes over
a nt change, but if a color has high quality, we may prefer one nt
change.
In the above example, the penalties of the two types of decoding are
q(B)+25 and q(B)+q(O), respectively. If q(O)>25, we prefer the first;
otherwise the second. Note that no matter what we choose, the fourth
base will get a low nt quality.
*/
#define COLOR_MM 19
#define NUCL_MM 25
static const int nst_ntnt2cs_table[] = { 4, 0, 0, 1, 0, 2, 3, 4, 0, 3, 2, 4, 1, 4, 4, 4 };
/*
{A,C,G,T,N} -> {0,1,2,3,4}
nt_ref[0..size]: nucleotide reference: 0/1/2/3/4
cs_read[0..size-1]: color read+qual sequence: base<<6|qual; qual==63 for N
nt_read[0..size]: nucleotide read sequence: 0/1/2/3 (returned)
btarray[0..4*size]: backtrack array (working space)
*/
void cs2nt_DP(int size, const uint8_t *nt_ref, const uint8_t *cs_read, uint8_t *nt_read, uint8_t *btarray)
{
int h[8], curr, last;
int x, y, xmin, hmin, k;
// h[0..3] and h[4..7] are the current and last best score array, depending on curr and last
// recursion: initial value
if (nt_ref[0] >= 4) memset(h, 0, sizeof(int) << 2);
else {
for (x = 0; x != 4; ++x) h[x] = NUCL_MM;
h[nt_ref[0]] = 0;
}
// recursion: main loop
curr = 1; last = 0;
for (k = 1; k <= size; ++k) {
for (x = 0; x != 4; ++x) {
int min = 0x7fffffff, ymin = 0;
for (y = 0; y != 4; ++y) {
int s = h[last<<2|y];
if ((cs_read[k-1]&0x3f) != 63 && cs_read[k-1]>>6 != nst_ntnt2cs_table[1<<x|1<<y])
s += ((cs_read[k-1]&0x3f) < COLOR_MM)? COLOR_MM : (cs_read[k-1]&0x3f); // color mismatch
if (nt_ref[k] < 4 && nt_ref[k] != x) s += NUCL_MM; // nt mismatch
if (s < min) {
min = s; ymin = y;
}
}
h[curr<<2|x] = min; btarray[k<<2|x] = ymin;
}
last = curr; curr = 1 - curr; // swap
}
// back trace
hmin = 0x7fffffff; xmin = 0;
for (x = 0; x != 4; ++x) {
if (h[last<<2|x] < hmin) {
hmin = h[last<<2|x]; xmin = x;
}
}
nt_read[size] = xmin;
for (k = size - 1; k >= 0; --k)
nt_read[k] = btarray[(k+1)<<2 | nt_read[k+1]];
}
/*
nt_read[0..size]: nucleotide read sequence: 0/1/2/3
cs_read[0..size-1]: color read+qual sequence: base<<6|qual; qual==63 for N
tarray[0..size*2-1]: temporary array
*/
uint8_t *cs2nt_nt_qual(int size, const uint8_t *nt_read, const uint8_t *cs_read, uint8_t *tarray)
{
int k, c1, c2;
uint8_t *t2array = tarray + size;
// get the color sequence of nt_read
c1 = nt_read[0];
for (k = 1; k <= size; ++k) {
c2 = nt_read[k]; // in principle, there is no 'N' in nt_read[]; just in case
tarray[k-1] = (c1 >= 4 || c2 >= 4)? 4 : nst_ntnt2cs_table[1<<c1 | 1<<c2];
c1 = c2;
}
for (k = 1; k != size; ++k) {
int q = 0;
if (tarray[k-1] == cs_read[k-1]>>6 && tarray[k] == cs_read[k]>>6) {
q = (int)(cs_read[k-1]&0x3f) + (int)(cs_read[k]&0x3f) + 10;
} else if (tarray[k-1] == cs_read[k-1]>>6) {
q = (int)(cs_read[k-1]&0x3f) - (int)(cs_read[k]&0x3f);
} else if (tarray[k] == cs_read[k]>>6) {
q = (int)(cs_read[k]&0x3f) - (int)(cs_read[k-1]&0x3f);
} // else, q = 0
if (q < 0) q = 0;
if (q > 60) q = 60;
t2array[k] = nt_read[k]<<6 | q;
if ((cs_read[k-1]&0x3f) == 63 || (cs_read[k]&0x3f) == 63) t2array[k] = 0;
}
return t2array + 1; // of size-2
}
// this function will be called when p->seq has been reversed by refine_gapped()
void bwa_cs2nt_core(bwa_seq_t *p, bwtint_t l_pac, ubyte_t *pac)
{
uint8_t *ta, *nt_read, *btarray, *tarray, *nt_ref, *cs_read, *new_nt_read;
int i, len;
uint8_t *seq;
// set temporary arrays
if (p->type == BWA_TYPE_NO_MATCH) return;
len = p->len + p->n_gapo + p->n_gape + 100; // leave enough space
ta = (uint8_t*)malloc(len * 7);
nt_ref = ta;
cs_read = nt_ref + len;
nt_read = cs_read + len;
btarray = nt_read + len;
tarray = nt_read + len;
#define __gen_csbase(_cs, _i, _seq) do { \
int q = p->qual[p->strand? p->len - 1 - (_i) : (_i)] - 33; \
if (q > 60) q = 60; \
if (_seq[_i] > 3) q = 63; \
(_cs) = _seq[_i]<<6 | q; \
} while (0)
// generate len, nt_ref[] and cs_read
seq = p->strand? p->rseq : p->seq;
nt_ref[0] = p->pos? bns_pac(pac, p->pos-1) : 4;
if (p->cigar == 0) { // no gap or clipping
len = p->len;
for (i = 0; i < p->len; ++i) {
__gen_csbase(cs_read[i], i, seq);
nt_ref[i+1] = bns_pac(pac, p->pos + i);
}
} else {
int k, z;
bwtint_t x, y;
x = p->pos; y = 0;
for (k = z = 0; k < p->n_cigar; ++k) {
int l = __cigar_len(p->cigar[k]);
if (__cigar_op(p->cigar[k]) == FROM_M) {
for (i = 0; i < l; ++i, ++x, ++y) {
__gen_csbase(cs_read[z], y, seq);
nt_ref[z+1] = bns_pac(pac, x);
++z;
}
} else if (__cigar_op(p->cigar[k]) == FROM_I) {
for (i = 0; i < l; ++i, ++y) {
__gen_csbase(cs_read[z], y, seq);
nt_ref[z+1] = 4;
++z;
}
} else if (__cigar_op(p->cigar[k]) == FROM_S) y += l;
else x += l;
}
len = z;
}
cs2nt_DP(len, nt_ref, cs_read, nt_read, btarray);
new_nt_read = cs2nt_nt_qual(len, nt_read, cs_read, tarray);
// update p
p->len = p->full_len = len - 1;
for (i = 0; i < p->len; ++i) {
if ((new_nt_read[i]&0x3f) == 63) {
p->qual[i] = 33; seq[i] = 4;
} else {
p->qual[i] = (new_nt_read[i]&0x3f) + 33;
seq[i] = new_nt_read[i]>>6;
}
}
p->qual[p->len] = seq[p->len] = 0;
if (p->strand) {
memcpy(p->seq, seq, p->len);
seq_reverse(p->len, p->seq, 1);
seq_reverse(p->len, p->qual, 0);
} else {
memcpy(p->rseq, seq, p->len);
seq_reverse(p->len, p->rseq, 1);
}
free(ta);
}