* simutaneously, to achieve random access and to use the BED interface.
* To compile this program separately, you may:
*
* simutaneously, to achieve random access and to use the BED interface.
* To compile this program separately, you may:
*
- int i, n, tid, beg, end, pos, *n_plp;
+ int i, n, tid, beg, end, pos, *n_plp, baseQ = 0;
const bam_pileup1_t **plp;
char *reg = 0; // specified region
void *bed = 0; // BED data structure
const bam_pileup1_t **plp;
char *reg = 0; // specified region
void *bed = 0; // BED data structure
switch (n) {
case 'r': reg = strdup(optarg); break; // parsing a region requires a BAM header
case 'b': bed = bed_read(optarg); break; // BED or position list file can be parsed now
switch (n) {
case 'r': reg = strdup(optarg); break; // parsing a region requires a BAM header
case 'b': bed = bed_read(optarg); break; // BED or position list file can be parsed now
// the core multi-pileup loop
mplp = bam_mplp_init(n, read_bam, (void**)data); // initialization
n_plp = calloc(n, sizeof(int)); // n_plp[i] is the number of covering reads from the i-th BAM
// the core multi-pileup loop
mplp = bam_mplp_init(n, read_bam, (void**)data); // initialization
n_plp = calloc(n, sizeof(int)); // n_plp[i] is the number of covering reads from the i-th BAM
while (bam_mplp_auto(mplp, &tid, &pos, n_plp, plp) > 0) { // come to the next covered position
if (pos < beg || pos >= end) continue; // out of range; skip
if (bed && bed_overlap(bed, h->target_name[tid], pos, pos + 1) == 0) continue; // not in BED; skip
while (bam_mplp_auto(mplp, &tid, &pos, n_plp, plp) > 0) { // come to the next covered position
if (pos < beg || pos >= end) continue; // out of range; skip
if (bed && bed_overlap(bed, h->target_name[tid], pos, pos + 1) == 0) continue; // not in BED; skip
- const bam_pileup1_t *p = plp[i];
- for (j = 0; j < n_plp[i]; ++j) // this loop counts #reads having deletions or refskip at tid:pos
- if (p->is_del || p->is_refskip) ++m;
+ for (j = 0; j < n_plp[i]; ++j) {
+ const bam_pileup1_t *p = plp[i] + j; // DON'T modfity plp[][] unless you really know
+ if (p->is_del || p->is_refskip) ++m; // having dels or refskips at tid:pos
+ else if (bam1_qual(p->b)[p->qpos] < baseQ) ++m; // low base quality
+ }