+#include <assert.h>
+#include "bam2bcf.h"
+#include "sample.h"
+
+#define MPLP_GLF 0x10
+#define MPLP_NO_COMP 0x20
+#define MPLP_NO_ORPHAN 0x40
+#define MPLP_REALN 0x80
+#define MPLP_FMT_DP 0x100
+#define MPLP_FMT_SP 0x200
+
+typedef struct {
+ int max_mq, min_mq, flag, min_baseQ, capQ_thres, max_depth;
+ char *reg, *fn_pos;
+ faidx_t *fai;
+ kh_64_t *hash;
+} mplp_conf_t;
+
+typedef struct {
+ bamFile fp;
+ bam_iter_t iter;
+ int min_mq, flag, ref_id, capQ_thres;
+ char *ref;
+} mplp_aux_t;
+
+typedef struct {
+ int n;
+ int *n_plp, *m_plp;
+ bam_pileup1_t **plp;
+} mplp_pileup_t;
+
+static int mplp_func(void *data, bam1_t *b)
+{
+ extern int bam_realn(bam1_t *b, const char *ref);
+ extern int bam_prob_realn(bam1_t *b, const char *ref);
+ extern int bam_cap_mapQ(bam1_t *b, char *ref, int thres);
+ mplp_aux_t *ma = (mplp_aux_t*)data;
+ int ret, skip = 0;
+ do {
+ int has_ref = (ma->ref && ma->ref_id == b->core.tid)? 1 : 0;
+ ret = ma->iter? bam_iter_read(ma->fp, ma->iter, b) : bam_read1(ma->fp, b);
+ if (ret < 0) break;
+ skip = 0;
+ if (has_ref && (ma->flag&MPLP_REALN)) bam_prob_realn(b, ma->ref);
+ if (has_ref && ma->capQ_thres > 10) {
+ int q = bam_cap_mapQ(b, ma->ref, ma->capQ_thres);
+ if (q < 0) skip = 1;
+ else if (b->core.qual > q) b->core.qual = q;
+ } else if (b->core.flag&BAM_FUNMAP) skip = 1;
+ else if (b->core.qual < ma->min_mq) skip = 1;
+ else if ((ma->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);
+ assert(id >= 0 && id < m->n);
+ 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;
+ }
+ }
+}
+
+static int mpileup(mplp_conf_t *conf, int n, char **fn)