From: Heng Li Date: Tue, 16 Nov 2010 18:51:33 +0000 (+0000) Subject: * samtools-0.1.9-18 (r825) X-Git-Url: https://git.donarmstrong.com/?a=commitdiff_plain;h=a27cb2ac9367fe6e425226721741172c174d7a78;p=samtools.git * 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 518768d..819c1e7 100644 --- 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 29b315f..e1e756d 100644 --- 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 7923b04..4db2453 100644 --- 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 17db537..e87bef7 100644 --- 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