+ int ret, skip = 0;
+ do {
+ int has_ref;
+ ret = ma->iter? bam_iter_read(ma->fp, ma->iter, b) : bam_read1(ma->fp, b);
+ if (ret < 0) break;
+ if (b->core.tid < 0 || (b->core.flag&BAM_FUNMAP)) { // exclude unmapped reads
+ skip = 1;
+ continue;
+ }
+ if (ma->conf->bed) { // test overlap
+ skip = !bed_overlap(ma->conf->bed, ma->h->target_name[b->core.tid], b->core.pos, bam_calend(&b->core, bam1_cigar(b)));
+ if (skip) continue;
+ }
+ if (ma->conf->rghash) { // exclude read groups
+ uint8_t *rg = bam_aux_get(b, "RG");
+ skip = (rg && bcf_str2id(ma->conf->rghash, (const char*)(rg+1)) >= 0);
+ if (skip) continue;
+ }
+ if (ma->conf->flag & MPLP_ILLUMINA13) {
+ int i;
+ uint8_t *qual = bam1_qual(b);
+ for (i = 0; i < b->core.l_qseq; ++i)
+ qual[i] = qual[i] > 31? qual[i] - 31 : 0;
+ }
+ has_ref = (ma->ref && ma->ref_id == b->core.tid)? 1 : 0;
+ skip = 0;
+ if (has_ref && (ma->conf->flag&MPLP_REALN)) bam_prob_realn_core(b, ma->ref, (ma->conf->flag & MPLP_EXT_BAQ)? 3 : 1);
+ if (has_ref && ma->conf->capQ_thres > 10) {
+ int q = bam_cap_mapQ(b, ma->ref, ma->conf->capQ_thres);
+ if (q < 0) skip = 1;
+ else if (b->core.qual > q) b->core.qual = q;
+ }
+ else if (b->core.qual < ma->conf->min_mq) skip = 1;
+ else if ((ma->conf->flag&MPLP_NO_ORPHAN) && (b->core.flag&1) && !(b->core.flag&2)) skip = 1;
+ } while (skip);
+ return ret;
+}
+
+static void group_smpl(mplp_pileup_t *m, bam_sample_t *sm, kstring_t *buf,
+ int n, char *const*fn, int *n_plp, const bam_pileup1_t **plp)
+{
+ int i, j;
+ memset(m->n_plp, 0, m->n * sizeof(int));
+ for (i = 0; i < n; ++i) {
+ for (j = 0; j < n_plp[i]; ++j) {
+ const bam_pileup1_t *p = plp[i] + j;
+ uint8_t *q;
+ int id = -1;
+ q = bam_aux_get(p->b, "RG");
+ if (q) id = bam_smpl_rg2smid(sm, fn[i], (char*)q+1, buf);
+ if (id < 0) id = bam_smpl_rg2smid(sm, fn[i], 0, buf);
+ if (id < 0 || id >= m->n) {
+ assert(q); // otherwise a bug
+ fprintf(stderr, "[%s] Read group %s used in file %s but not defined in the header.\n", __func__, (char*)q+1, fn[i]);
+ exit(1);
+ }
+ if (m->n_plp[id] == m->m_plp[id]) {
+ m->m_plp[id] = m->m_plp[id]? m->m_plp[id]<<1 : 8;
+ m->plp[id] = realloc(m->plp[id], sizeof(bam_pileup1_t) * m->m_plp[id]);
+ }
+ m->plp[id][m->n_plp[id]++] = *p;
+ }
+ }