* samtools-0.1.9-18 (r825)

 * change the behaviour of calmd such that by default it does not change the base quality

diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c
index 518768d4af3095d4ad178539acae0cbb0ed736bb..819c1e7b1371296b7da5536939fab691f53c1722 100644 (file)
--- a/bam2bcf_indel.c
+++ b/bam2bcf_indel.c
@@ -286,7 +286,7 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla
                                        const uint8_t *qual = bam1_qual(p->b), *bq;
                                        uint8_t *qq;
                                        qq = calloc(qend - qbeg, 1);
-                                       bq = (uint8_t*)bam_aux_get(p->b, "BQ");
+                                       bq = (uint8_t*)bam_aux_get(p->b, "ZQ");
                                        if (bq) ++bq; // skip type
                                        for (l = qbeg; l < qend; ++l) {
                                                qq[l - qbeg] = bq? qual[l] + (bq[l] - 64) : qual[l];

diff --git a/bam_md.c b/bam_md.c
index 29b315f4e9235f4e83516d6fb544821185d19a82..e1e756d1d77907a9d67e9bb65337cc87705338c0 100644 (file)
--- a/bam_md.c
+++ b/bam_md.c
@@ -162,14 +162,31 @@ int bam_cap_mapQ(bam1_t *b, char *ref, int thres)
        return (int)(t + .499);
 }
 
-int bam_prob_realn_core(bam1_t *b, const char *ref, int write_bq)
+int bam_prob_realn_core(bam1_t *b, const char *ref, int apply_baq)
 {
        int k, i, bw, x, y, yb, ye, xb, xe;
        uint32_t *cigar = bam1_cigar(b);
        bam1_core_t *c = &b->core;
        kpa_par_t conf = kpa_par_def;
-       if (c->flag & BAM_FUNMAP) return -1;
-       if (bam_aux_get(b, "BQ")) return -2;
+       uint8_t *bq = 0, *zq = 0, *qual = bam1_qual(b);
+       if (c->flag & BAM_FUNMAP) return -1; // do nothing
+       // test if BQ or ZQ is present
+       if ((bq = bam_aux_get(b, "BQ")) != 0) ++bq;
+       if ((zq = bam_aux_get(b, "ZQ")) != 0 && *zq == 'Z') ++zq;
+       if (bq && zq) return -2; // do not know what to do...
+       if (bq || zq) {
+               if ((apply_baq && zq) || (!apply_baq && bq)) return -3; // in both cases, do nothing
+               if (bq && apply_baq) { // then convert BQ to ZQ
+                       for (i = 0; i < c->l_qseq; ++i)
+                               qual[i] = qual[i] + 64 < bq[i]? 0 : qual[i] - ((int)bq[i] - 64);
+                       *(bq - 3) = 'Z';
+               } else if (zq && !apply_baq) { // then convert ZQ to BQ
+                       for (i = 0; i < c->l_qseq; ++i)
+                               qual[i] += (int)zq[i] - 64;
+                       *(zq - 3) = 'B';
+               }
+               return 0;
+       }
        // find the start and end of the alignment      
        x = c->pos, y = 0, yb = ye = xb = xe = -1;
        for (k = 0; k < c->n_cigar; ++k) {
@@ -182,7 +199,7 @@ int bam_prob_realn_core(bam1_t *b, const char *ref, int write_bq)
                        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;
+               else if (op == BAM_CREF_SKIP) return -1; // do nothing if there is a reference skip
        }
        // set bandwidth and the start and the end
        bw = 7;
@@ -194,12 +211,10 @@ int bam_prob_realn_core(bam1_t *b, const char *ref, int write_bq)
        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), *bq = 0;
+               uint8_t *s, *r, *q, *seq = bam1_seq(b), *bq;
                int *state;
-               if (write_bq) {
-                       bq = calloc(c->l_qseq + 1, 1);
-                       memcpy(bq, qual, c->l_qseq);
-               }
+               bq = calloc(c->l_qseq + 1, 1);
+               memcpy(bq, qual, c->l_qseq);
                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);
@@ -212,40 +227,40 @@ int bam_prob_realn_core(bam1_t *b, const char *ref, int write_bq)
                        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];
+                                       if ((state[i]&3) != 0 || state[i]>>2 != x - xb + (i - y)) bq[i] = 0;
+                                       else bq[i] = bq[i] < q[i]? bq[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;
                }
-               if (write_bq) {
-                       for (i = 0; i < c->l_qseq; ++i) bq[i] = bq[i] - qual[i] + 64;
-                       bam_aux_append(b, "BQ", 'Z', c->l_qseq + 1, bq);
-                       free(bq);
-               }
-               free(s); free(r); free(q); free(state);
+               for (i = 0; i < c->l_qseq; ++i) bq[i] = qual[i] - bq[i] + 64; // finalize BQ
+               if (apply_baq) {
+                       for (i = 0; i < c->l_qseq; ++i) qual[i] -= bq[i] - 64; // modify qual
+                       bam_aux_append(b, "ZQ", 'Z', c->l_qseq + 1, bq);
+               } else bam_aux_append(b, "BQ", 'Z', c->l_qseq + 1, bq);
+               free(bq); free(s); free(r); free(q); free(state);
        }
        return 0;
 }
 
 int bam_prob_realn(bam1_t *b, const char *ref)
 {
-       return bam_prob_realn_core(b, ref, 0);
+       return bam_prob_realn_core(b, ref, 1);
 }
 
 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, capQ = 0;
+       int c, is_equal, tid = -2, ret, len, is_bam_out, is_sam_in, is_uncompressed, max_nm, is_realn, capQ, apply_baq;
        samfile_t *fp, *fpout = 0;
        faidx_t *fai;
        char *ref = 0, mode_w[8], mode_r[8];
        bam1_t *b;
 
-       is_bam_out = is_sam_in = is_uncompressed = is_realn = 0;
+       is_equal = is_bam_out = is_sam_in = is_uncompressed = is_realn = max_nm = apply_baq = capQ = 0;
        mode_w[0] = mode_r[0] = 0;
        strcpy(mode_r, "r"); strcpy(mode_w, "w");
-       while ((c = getopt(argc, argv, "reubSC:n:")) >= 0) {
+       while ((c = getopt(argc, argv, "reubSC:n:A")) >= 0) {
                switch (c) {
                case 'r': is_realn = 1; break;
                case 'e': is_equal = 1; break;
@@ -254,6 +269,7 @@ int bam_fillmd(int argc, char *argv[])
                case 'S': is_sam_in = 1; break;
                case 'n': max_nm = atoi(optarg); break;
                case 'C': capQ = atoi(optarg); break;
+               case 'A': apply_baq = 1; break;
                default: fprintf(stderr, "[bam_fillmd] unrecognized option '-%c'\n", c); return 1;
                }
        }
@@ -268,6 +284,7 @@ int bam_fillmd(int argc, char *argv[])
                fprintf(stderr, "         -u       uncompressed BAM output (for piping)\n");
                fprintf(stderr, "         -b       compressed BAM output\n");
                fprintf(stderr, "         -S       the input is SAM with header\n");
+               fprintf(stderr, "         -A       modify the quality string\n");
                fprintf(stderr, "         -r       read-independent local realignment\n\n");
                return 1;
        }
@@ -291,7 +308,7 @@ 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_prob_realn_core(b, ref, 1);
+                       if (is_realn) bam_prob_realn_core(b, ref, apply_baq);
                        if (capQ > 10) {
                                int q = bam_cap_mapQ(b, ref, capQ);
                                if (b->core.qual > q) b->core.qual = q;

diff --git a/bamtk.c b/bamtk.c
index 7923b04007fbce7359c504bb6d105b4e518bd129..4db24532bd14f70c06817ecece5788b44429df06 100644 (file)
--- a/bamtk.c
+++ b/bamtk.c
@@ -9,7 +9,7 @@
 #endif
 
 #ifndef PACKAGE_VERSION
-#define PACKAGE_VERSION "0.1.9-17 (r824)"
+#define PACKAGE_VERSION "0.1.9-18 (r825)"
 #endif
 
 int bam_taf2baf(int argc, char *argv[]);

diff --git a/samtools.1 b/samtools.1
index 17db53758b2e97c95b9ee450f14784a0ed0a7d74..e87bef7e788f09d1c0c38219f231062868a010d4 100644 (file)
--- a/samtools.1
+++ b/samtools.1
@@ -364,6 +364,11 @@ tag. Output SAM by default.
 .B OPTIONS:
 .RS
 .TP 8
+.B -A
+When used jointly with
+.B -r
+this option overwrites the original base quality.
+.TP 8
 .B -e
 Convert a the read base to = if it is identical to the aligned reference
 base. Indel caller does not support the = bases at the moment.
@@ -383,8 +388,7 @@ Coefficient to cap mapping quality of poorly mapped reads. See the
 command for details. [0]
 .TP
 .B -r
-Perform probabilistic realignment to compute BAQ, which will be used to
-cap base quality.
+Compute the BQ tag without changing the base quality.
 .RE
 
 .TP
@@ -666,7 +670,7 @@ commands estimate AFS by EM.
 .IP o 2
 Dump BAQ applied alignment for other SNP callers:
 
-  samtools calmd -br aln.bam > aln.baq.bam
+  samtools calmd -bAr aln.bam > aln.baq.bam
 
 It adds and corrects the
 .B NM