From 9360df156b23158232dca10ffdb8a8ffdcaf0930 Mon Sep 17 00:00:00 2001 From: Heng Li Date: Wed, 10 Jun 2009 14:56:17 +0000 Subject: [PATCH] * samtools-0.1.4-3 (r324) * make maqcns generate the correct call in repetitive regions. * allow filtering on mapQ at the pileup command line --- bam_maqcns.c | 22 ++++++++++++++++++++++ bam_plcmd.c | 8 +++++--- bamtk.c | 2 +- sam.c | 5 +++-- sam.h | 2 +- 5 files changed, 32 insertions(+), 7 deletions(-) diff --git a/bam_maqcns.c b/bam_maqcns.c index 95cd2b5..44eed4e 100644 --- a/bam_maqcns.c +++ b/bam_maqcns.c @@ -213,6 +213,28 @@ glf1_t *bam_maqcns_glfgen(int _n, const bam_pileup1_t *pl, uint8_t ref_base, bam if (p[j<<2|k] < 0.0) p[j<<2|k] = 0.0; } + { // fix p[k<<2|k] + float max1, max2, min; + int max_k, min_k; + max_k = min_k = -1; + max1 = max2 = -1.0; min = 1e30; + for (k = 0; k < 4; ++k) { + if (b->esum[k] > max1) { + max2 = max1; max1 = b->esum[k]; max_k = k; + } else if (b->esum[k] > max2) { + max2 = b->esum[k]; + } + } + for (k = 0; k < 4; ++k) { + if (min > p[k<<2|k]) { + min = p[k<<2|k]; + min_k = k; + } + } + if (max1 > max2 && min_k != max_k) + p[max_k<<2|max_k] = min > 1.0? min - 1.0 : 0.0; + } + // convert necessary information to glf1_t g->ref_base = ref_base; g->max_mapQ = b->rms_mapQ; g->depth = n > 16777215? 16777215 : n; diff --git a/bam_plcmd.c b/bam_plcmd.c index 46f078a..96f5c45 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -26,7 +26,7 @@ typedef struct { khash_t(64) *hash; uint32_t format; int tid, len, last_pos; - int mask; + int mask, min_mapQ; char *ref; glfFile fp_glf; // for glf output only } pu_data_t; @@ -280,7 +280,7 @@ int bam_pileup(int argc, char *argv[]) d->tid = -1; d->mask = BAM_DEF_MASK; d->c = bam_maqcns_init(); d->ido = bam_maqindel_opt_init(); - while ((c = getopt(argc, argv, "st:f:cT:N:r:l:im:gI:G:vM:")) >= 0) { + while ((c = getopt(argc, argv, "st:f:cT:N:r:l:im:gI:G:vM:q:")) >= 0) { switch (c) { case 's': d->format |= BAM_PLF_SIMPLE; break; case 't': fn_list = strdup(optarg); break; @@ -297,6 +297,7 @@ int bam_pileup(int argc, char *argv[]) case 'g': d->format |= BAM_PLF_GLF; break; case 'I': d->ido->q_indel = atoi(optarg); break; case 'G': d->ido->r_indel = atof(optarg); break; + case 'q': d->min_mapQ = atoi(optarg); break; default: fprintf(stderr, "Unrecognizd option '-%c'.\n", c); return 1; } } @@ -306,6 +307,7 @@ int bam_pileup(int argc, char *argv[]) fprintf(stderr, "Option: -s simple (yet incomplete) pileup format\n"); fprintf(stderr, " -i only show lines/consensus with indels\n"); fprintf(stderr, " -m INT filtering reads with bits in INT [%d]\n", d->mask); + fprintf(stderr, " -q INT filtering reads with mapping quality smaller than INT [%d]\n", d->min_mapQ); fprintf(stderr, " -M INT cap mapping quality at INT [%d]\n", d->c->cap_mapQ); fprintf(stderr, " -t FILE list of reference sequences (assume the input is in SAM)\n"); fprintf(stderr, " -l FILE list of sites at which pileup is output\n"); @@ -342,7 +344,7 @@ int bam_pileup(int argc, char *argv[]) } d->h = fp->header; if (fn_pos) d->hash = load_pos(fn_pos, d->h); - sampileup(fp, d->mask, pileup_func, d); + sampileup(fp, d->mask, d->min_mapQ, pileup_func, d); samclose(fp); // d->h will be destroyed here } diff --git a/bamtk.c b/bamtk.c index 0fb7cb4..8441082 100644 --- a/bamtk.c +++ b/bamtk.c @@ -3,7 +3,7 @@ #include "bam.h" #ifndef PACKAGE_VERSION -#define PACKAGE_VERSION "0.1.4-2 (r316)" +#define PACKAGE_VERSION "0.1.4-3 (r324)" #endif int bam_taf2baf(int argc, char *argv[]); diff --git a/sam.c b/sam.c index 48a4b34..9a6e201 100644 --- a/sam.c +++ b/sam.c @@ -116,7 +116,7 @@ int samwrite(samfile_t *fp, const bam1_t *b) } } -int sampileup(samfile_t *fp, int mask, bam_pileup_f func, void *func_data) +int sampileup(samfile_t *fp, int mask, int min_mapQ, bam_pileup_f func, void *func_data) { bam_plbuf_t *buf; int ret; @@ -125,7 +125,8 @@ int sampileup(samfile_t *fp, int mask, bam_pileup_f func, void *func_data) buf = bam_plbuf_init(func, func_data); bam_plbuf_set_mask(buf, mask); while ((ret = samread(fp, b)) >= 0) - bam_plbuf_push(b, buf); + if (b->core.qual >= min_mapQ) + bam_plbuf_push(b, buf); bam_plbuf_push(0, buf); bam_plbuf_destroy(buf); bam_destroy1(b); diff --git a/sam.h b/sam.h index 93daea6..082fe96 100644 --- a/sam.h +++ b/sam.h @@ -22,7 +22,7 @@ extern "C" { void samclose(samfile_t *fp); int samread(samfile_t *fp, bam1_t *b); int samwrite(samfile_t *fp, const bam1_t *b); - int sampileup(samfile_t *fp, int mask, bam_pileup_f func, void *func_data); + int sampileup(samfile_t *fp, int mask, int min_mapQ, bam_pileup_f func, void *func_data); #ifdef __cplusplus } -- 2.39.5