X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bam_md.c;h=f65b845b5d594980d4b4dba36122672c71b19cc9;hb=3bedbd7fb801e7d0ffdedf2eb53dcce9ff08d31a;hp=fb4673263044cd80831288a1b7a286ec62486d14;hpb=221f82f662dd770a17d1bd2181de46b8c3b3fdfd;p=samtools.git diff --git a/bam_md.c b/bam_md.c index fb46732..f65b845 100644 --- a/bam_md.c +++ b/bam_md.c @@ -2,6 +2,7 @@ #include #include #include +#include #include "faidx.h" #include "sam.h" #include "kstring.h" @@ -109,60 +110,117 @@ void bam_fillmd1(bam1_t *b, char *ref, int is_equal) bam_fillmd1_core(b, ref, is_equal, 0); } -// local realignment - -#define MIN_REF_LEN 10 - -int bam_realn(bam1_t *b, const char *ref) +int bam_cap_mapQ(bam1_t *b, char *ref, int thres) { - int k, l_ref, score, n_cigar; + uint8_t *seq = bam1_seq(b), *qual = bam1_qual(b); uint32_t *cigar = bam1_cigar(b); - uint8_t *s_ref = 0, *s_read = 0, *seq; - ka_param_t par; bam1_core_t *c = &b->core; - // set S/W parameters - par = ka_param_blast; par.gap_open = 4; par.gap_ext = 1; par.gap_end_open = par.gap_end_ext = 0; - // calculate the length of the reference in the alignment - for (k = l_ref = 0; k < c->n_cigar; ++k) { - if ((cigar[k]&0xf) == BAM_CREF_SKIP) break; // do not do realignment if there is an `N' operation - if ((cigar[k]&0xf) == BAM_CMATCH || (cigar[k]&0xf) == BAM_CDEL) - l_ref += cigar[k]>>4; - } - if (k != c->n_cigar || l_ref < MIN_REF_LEN) return -1; - for (k = 0; k < l_ref; ++k) - if (ref[c->pos + k] == 0) return -1; // the read stands out of the reference - // allocate - s_ref = calloc(l_ref, 1); - s_read = calloc(c->l_qseq, 1); - for (k = 0, seq = bam1_seq(b); k < c->l_qseq; ++k) - s_read[k] = bam_nt16_nt4_table[bam1_seqi(seq, k)]; - for (k = 0; k < l_ref; ++k) - s_ref[k] = bam_nt16_nt4_table[(int)bam_nt16_table[(int)ref[c->pos + k]]]; - // do alignment - cigar = ka_global_core(s_ref, l_ref, s_read, c->l_qseq, &par, &score, &n_cigar); - if (score <= 0) { // realignment failed - free(cigar); free(s_ref); free(s_read); - return -1; + int i, x, y, mm, q, len, clip_l, clip_q; + double t; + if (thres < 0) thres = 40; // set the default + mm = q = len = clip_l = clip_q = 0; + for (i = y = 0, x = c->pos; i < c->n_cigar; ++i) { + int j, l = cigar[i]>>4, op = cigar[i]&0xf; + if (op == BAM_CMATCH) { + for (j = 0; j < l; ++j) { + int z = y + j; + int c1 = bam1_seqi(seq, z), c2 = bam_nt16_table[(int)ref[x+j]]; + if (ref[x+j] == 0) break; // out of boundary + if (c2 != 15 && c1 != 15 && qual[z] >= 13) { // not ambiguous + ++len; + if (c1 && c1 != c2 && qual[z] >= 13) { // mismatch + ++mm; + q += qual[z] > 33? 33 : qual[z]; + } + } + } + if (j < l) break; + x += l; y += l; len += l; + } else if (op == BAM_CDEL) { + for (j = 0; j < l; ++j) + if (ref[x+j] == 0) break; + if (j < l) break; + x += l; + } else if (op == BAM_CSOFT_CLIP) { + for (j = 0; j < l; ++j) clip_q += qual[y+j]; + clip_l += l; + y += l; + } else if (op == BAM_CHARD_CLIP) { + clip_q += 13 * l; + clip_l += l; + } else if (op == BAM_CINS) y += l; + else if (op == BAM_CREF_SKIP) x += l; } - // copy over the alignment - if (4 * (n_cigar - (int)c->n_cigar) + b->data_len > b->m_data) { // enlarge b->data - b->m_data = 4 * (n_cigar - (int)c->n_cigar) + b->data_len; - kroundup32(b->m_data); - b->data = realloc(b->data, b->m_data); + for (i = 0, t = 1; i < mm; ++i) + t *= (double)len / (i+1); + t = q - 4.343 * log(t) + clip_q / 5.; + if (t > thres) return -1; + if (t < 0) t = 0; + t = sqrt((thres - t) / thres) * thres; +// fprintf(stderr, "%s %lf %d\n", bam1_qname(b), t, q); + return (int)(t + .499); +} + +int bam_prob_realn(bam1_t *b, const char *ref) +{ + int k, i, bw, x, y, yb, ye, xb, xe; + uint32_t *cigar = bam1_cigar(b); + bam1_core_t *c = &b->core; + ka_probpar_t conf = ka_probpar_def; + // find the start and end of the alignment + if (c->flag & BAM_FUNMAP) return -1; + x = c->pos, y = 0, yb = ye = xb = xe = -1; + for (k = 0; k < c->n_cigar; ++k) { + int op, l; + op = cigar[k]&0xf; l = cigar[k]>>4; + if (op == BAM_CMATCH) { + if (yb < 0) yb = y; + if (xb < 0) xb = x; + ye = y + l; xe = x + l; + x += l; y += l; + } else if (op == BAM_CSOFT_CLIP || op == BAM_CINS) y += l; + else if (op == BAM_CDEL) x += l; + else if (op == BAM_CREF_SKIP) return -1; } - if (n_cigar != (int)c->n_cigar) { // move data - memmove(b->data + c->l_qname + 4 * n_cigar, bam1_seq(b), b->data_len - c->l_qname - 4 * c->n_cigar); - b->data_len += 4 * (n_cigar - (int)c->n_cigar); + // set bandwidth and the start and the end + bw = 7; + if (abs((xe - xb) - (ye - yb)) > bw) + bw = abs((xe - xb) - (ye - yb)) + 3; + conf.bw = bw; + xb -= yb + bw/2; if (xb < 0) xb = 0; + xe += c->l_qseq - ye + bw/2; + if (xe - xb - c->l_qseq > bw) + xb += (xe - xb - c->l_qseq - bw) / 2, xe -= (xe - xb - c->l_qseq - bw) / 2; + { // glocal + uint8_t *s, *r, *q, *seq = bam1_seq(b), *qual = bam1_qual(b); + int *state; + s = calloc(c->l_qseq, 1); + for (i = 0; i < c->l_qseq; ++i) s[i] = bam_nt16_nt4_table[bam1_seqi(seq, i)]; + r = calloc(xe - xb, 1); + for (i = xb; i < xe; ++i) + r[i-xb] = bam_nt16_nt4_table[bam_nt16_table[(int)ref[i]]]; + state = calloc(c->l_qseq, sizeof(int)); + q = calloc(c->l_qseq, 1); + ka_prob_glocal(r, xe-xb, s, c->l_qseq, qual, &conf, state, q); + for (k = 0, x = c->pos, y = 0; k < c->n_cigar; ++k) { + int op = cigar[k]&0xf, l = cigar[k]>>4; + if (op == BAM_CMATCH) { + for (i = y; i < y + l; ++i) { + if ((state[i]&3) != 0 || state[i]>>2 != x - xb + (i - y)) qual[i] = 0; + else qual[i] = qual[i] < q[i]? qual[i] : q[i]; + } + x += l; y += l; + } else if (op == BAM_CSOFT_CLIP || op == BAM_CINS) y += l; + else if (op == BAM_CDEL) x += l; + } + free(s); free(r); free(q); free(state); } - memcpy(bam1_cigar(b), cigar, n_cigar * 4); - c->n_cigar = n_cigar; - free(s_ref); free(s_read); free(cigar); return 0; } int bam_fillmd(int argc, char *argv[]) { - int c, is_equal = 0, tid = -2, ret, len, is_bam_out, is_sam_in, is_uncompressed, max_nm = 0, is_realn; + int c, is_equal = 0, tid = -2, ret, len, is_bam_out, is_sam_in, is_uncompressed, max_nm = 0, is_realn, capQ = 0; samfile_t *fp, *fpout = 0; faidx_t *fai; char *ref = 0, mode_w[8], mode_r[8]; @@ -171,7 +229,7 @@ int bam_fillmd(int argc, char *argv[]) is_bam_out = is_sam_in = is_uncompressed = is_realn = 0; mode_w[0] = mode_r[0] = 0; strcpy(mode_r, "r"); strcpy(mode_w, "w"); - while ((c = getopt(argc, argv, "reubSn:")) >= 0) { + while ((c = getopt(argc, argv, "reubSC:n:")) >= 0) { switch (c) { case 'r': is_realn = 1; break; case 'e': is_equal = 1; break; @@ -179,6 +237,7 @@ int bam_fillmd(int argc, char *argv[]) case 'u': is_uncompressed = is_bam_out = 1; break; case 'S': is_sam_in = 1; break; case 'n': max_nm = atoi(optarg); break; + case 'C': capQ = atoi(optarg); break; default: fprintf(stderr, "[bam_fillmd] unrecognized option '-%c'\n", c); return 1; } } @@ -216,7 +275,12 @@ int bam_fillmd(int argc, char *argv[]) fprintf(stderr, "[bam_fillmd] fail to find sequence '%s' in the reference.\n", fp->header->target_name[tid]); } - if (is_realn) bam_realn(b, ref); +// if (is_realn) bam_realn(b, ref); + if (is_realn) bam_prob_realn(b, ref); + if (capQ > 10) { + int q = bam_cap_mapQ(b, ref, capQ); + if (b->core.qual > q) b->core.qual = q; + } if (ref) bam_fillmd1_core(b, ref, is_equal, max_nm); } samwrite(fpout, b);