* make maqcns generate the correct call in repetitive regions.
* allow filtering on mapQ at the pileup command line
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;
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;
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;
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;
}
}
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");
}
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
}
#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[]);
}
}
-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;
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);
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
}