From 1a1147a5aef19776be078f3be4bd4cc7b9aa8ceb Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 10 Nov 2010 18:23:02 +0000 Subject: [PATCH] * samtools-0.1.9-8 (r809) * fixed a bug in the indel caller --- bam2bcf.c | 5 ++++- bam2bcf_indel.c | 9 +++++---- bam_plcmd.c | 13 +++++++------ bamtk.c | 2 +- 4 files changed, 17 insertions(+), 12 deletions(-) diff --git a/bam2bcf.c b/bam2bcf.c index 1ba288c..8537e23 100644 --- a/bam2bcf.c +++ b/bam2bcf.c @@ -124,7 +124,10 @@ int bcf_call_combine(int n, const bcf_callret1_t *calls, int ref_base /*4-bit*/, if (((ref4 < 4 && j < 4) || (ref4 == 4 && j < 5)) && i >= 0) call->unseen = j, call->a[j++] = qsum[i]&3; call->n_alleles = j; - } else call->n_alleles = j; + } else { + call->n_alleles = j; + if (call->n_alleles == 1) return -1; // no reliable supporting read. stop doing anything + } // set the PL array if (call->n < n) { call->n = n; diff --git a/bam2bcf_indel.c b/bam2bcf_indel.c index 2d0b5c6..c28b685 100644 --- a/bam2bcf_indel.c +++ b/bam2bcf_indel.c @@ -110,7 +110,7 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla const void *rghash) { extern void ks_introsort_uint32_t(int, uint32_t*); - int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score, N, K, l_run, ref_type; + int i, s, j, k, t, n_types, *types, max_rd_len, left, right, max_ins, *score, N, K, l_run, ref_type, n_alt; char *inscns = 0, *ref2, *query; khash_t(rg) *hash = (khash_t(rg)*)rghash; if (ref == 0 || bca == 0) return -1; @@ -342,13 +342,14 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla memcpy(&bca->inscns[t * bca->maxins], &inscns[(sumq[t]&0x3f) * max_ins], bca->maxins); } // update p->aux - for (s = K = 0; s < n; ++s) { - for (i = 0; i < n_plp[s]; ++i, ++K) { + for (s = n_alt = 0; s < n; ++s) { + for (i = 0; i < n_plp[s]; ++i) { bam_pileup1_t *p = plp[s] + i; int x = types[p->aux>>16&0x3f]; for (j = 0; j < 4; ++j) if (x == bca->indel_types[j]) break; p->aux = j<<16 | (j == 4? 0 : (p->aux&0xffff)); + if ((p->aux>>16&0x3f) > 0) ++n_alt; // fprintf(stderr, "pos=%d read=%d:%d name=%s call=%d q=%d\n", pos, s, i, bam1_qname(p->b), p->aux>>16&63, p->aux&0xff); } } @@ -356,5 +357,5 @@ int bcf_call_gap_prep(int n, int *n_plp, bam_pileup1_t **plp, int pos, bcf_calla free(score); // free free(types); free(inscns); - return 0; + return n_alt > 0? 0 : -1; } diff --git a/bam_plcmd.c b/bam_plcmd.c index 0fbc930..1dc819a 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -739,12 +739,13 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn) if (bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref, rghash) >= 0) { for (i = 0; i < gplp.n; ++i) bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], -1, bca, bcr + i); - bcf_call_combine(gplp.n, bcr, -1, &bc); - b = calloc(1, sizeof(bcf1_t)); - bcf_call2bcf(tid, pos, &bc, b, (conf->flag&(MPLP_FMT_DP|MPLP_FMT_SP))? bcr : 0, - (conf->flag&MPLP_FMT_SP), bca, ref); - bcf_write(bp, bh, b); - bcf_destroy(b); + if (bcf_call_combine(gplp.n, bcr, -1, &bc) >= 0) { + b = calloc(1, sizeof(bcf1_t)); + bcf_call2bcf(tid, pos, &bc, b, (conf->flag&(MPLP_FMT_DP|MPLP_FMT_SP))? bcr : 0, + (conf->flag&MPLP_FMT_SP), bca, ref); + bcf_write(bp, bh, b); + bcf_destroy(b); + } } } else { printf("%s\t%d\t%c", h->target_name[tid], pos + 1, (ref && pos < ref_len)? ref[pos] : 'N'); diff --git a/bamtk.c b/bamtk.c index d2817b9..04a7899 100644 --- a/bamtk.c +++ b/bamtk.c @@ -9,7 +9,7 @@ #endif #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.1.9-8 (r807)" +#define PACKAGE_VERSION "0.1.9-8 (r809)" #endif int bam_taf2baf(int argc, char *argv[]); -- 2.39.2