- // pileup strings
- printf("%d\t", n);
- for (i = 0; i < n; ++i) {
- const bam_pileup1_t *p = pu + i;
- if (max_mapq < p->b->core.qual) max_mapq = p->b->core.qual;
- if (p->is_head) printf("^%c", p->b->core.qual > 93? 126 : p->b->core.qual + 33);
- if (!p->is_del) {
- int c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos)];
- if (toupper(c) == toupper(rb)) c = bam1_strand(p->b)? ',' : '.';
- else c = bam1_strand(p->b)? tolower(c) : toupper(c);
- putchar(c);
- if (p->indel > 0) {
- printf("+%d", p->indel);
- for (j = 1; j <= p->indel; ++j) {
- c = bam_nt16_rev_table[bam1_seqi(bam1_seq(p->b), p->qpos + j)];
- putchar(bam1_strand(p->b)? tolower(c) : toupper(c));
+ max_indel_depth = conf->max_indel_depth * sm->n;
+ bam_mplp_set_maxcnt(iter, max_depth);
+ while (bam_mplp_auto(iter, &tid, &pos, n_plp, plp) > 0) {
+ if (conf->reg && (pos < beg0 || pos >= end0)) continue; // out of the region requested
+ if (conf->bed && tid >= 0 && !bed_overlap(conf->bed, h->target_name[tid], pos, pos+1)) continue;
+ if (tid != ref_tid) {
+ free(ref); ref = 0;
+ if (conf->fai) ref = faidx_fetch_seq(conf->fai, h->target_name[tid], 0, 0x7fffffff, &ref_len);
+ for (i = 0; i < n; ++i) data[i]->ref = ref, data[i]->ref_id = tid;
+ ref_tid = tid;
+ }
+ if (conf->flag & MPLP_GLF) {
+ int total_depth, _ref0, ref16;
+ bcf1_t *b = calloc(1, sizeof(bcf1_t));
+ for (i = total_depth = 0; i < n; ++i) total_depth += n_plp[i];
+ group_smpl(&gplp, sm, &buf, n, fn, n_plp, plp, conf->flag & MPLP_IGNORE_RG);
+ _ref0 = (ref && pos < ref_len)? ref[pos] : 'N';
+ ref16 = bam_nt16_table[_ref0];
+ for (i = 0; i < gplp.n; ++i)
+ bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], ref16, bca, bcr + i);
+ bcf_call_combine(gplp.n, bcr, ref16, &bc);
+ bcf_call2bcf(tid, pos, &bc, b, bcr, conf->fmt_flag, 0, 0);
+ bcf_write(bp, bh, b);
+ bcf_destroy(b);
+ // call indels
+ if (!(conf->flag&MPLP_NO_INDEL) && total_depth < max_indel_depth && 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);
+ if (bcf_call_combine(gplp.n, bcr, -1, &bc) >= 0) {
+ b = calloc(1, sizeof(bcf1_t));
+ bcf_call2bcf(tid, pos, &bc, b, bcr, conf->fmt_flag, 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');
+ for (i = 0; i < n; ++i) {
+ int j, cnt;
+ for (j = cnt = 0; j < n_plp[i]; ++j) {
+ const bam_pileup1_t *p = plp[i] + j;
+ if (bam1_qual(p->b)[p->qpos] >= conf->min_baseQ) ++cnt;