- }
- if (d->format & BAM_PLF_2ND) { // print 2nd calls and qualities
- const unsigned char *q;
- putchar('\t');
- for (i = 0; i < n; ++i) {
- const bam_pileup1_t *p = pu + i;
- q = bam_aux_get(p->b, "E2");
- putchar(q? q[p->qpos + 1] : 'N');
- }
- putchar('\t');
- for (i = 0; i < n; ++i) {
- const bam_pileup1_t *p = pu + i;
- q = bam_aux_get(p->b, "U2");
- putchar(q? q[p->qpos + 1] : '!');
- }
- }
- // print mapping quality if -s is flagged on the command line
- if (d->format & BAM_PLF_SIMPLE) {
- putchar('\t');
- for (i = 0; i < n; ++i) {
- int c = pu[i].b->core.qual + 33;
- if (c > 126) c = 126;
- putchar(c);
- }
- }
- // print read position
- if (d->format & BAM_PLF_READPOS) {
- putchar('\t');
- for (i = 0; i < n; ++i) {
- int x = pu[i].qpos;
- int l = pu[i].b->core.l_qseq;
- printf("%d,", x < l/2? x+1 : -((l-1)-x+1));
- }
- }
- putchar('\n');
- // print the indel line if r has been calculated. This only happens if:
- // a) -c or -i are flagged, AND b) the reference sequence is available
- if (r) {
- printf("%s\t%d\t*\t", d->h->target_name[tid], pos + 1);
- if (r->gt < 2) printf("%s/%s\t", r->s[r->gt], r->s[r->gt]);
- else printf("%s/%s\t", r->s[0], r->s[1]);
- printf("%d\t%d\t", r->q_cns, r->q_ref);
- printf("%d\t%d\t", rms_mapq, n);
- printf("%s\t%s\t", r->s[0], r->s[1]);
- //printf("%d\t%d\t", r->gl[0], r->gl[1]);
- printf("%d\t%d\t%d\t", r->cnt1, r->cnt2, r->cnt_anti);
- printf("%d\t%d\n", r->cnt_ref, r->cnt_ambi);
- bam_maqindel_ret_destroy(r);
- }
- return 0;
-}
-
-int bam_pileup(int argc, char *argv[])
-{
- int c, is_SAM = 0;
- char *fn_list = 0, *fn_fa = 0, *fn_pos = 0;
- pu_data_t *d = (pu_data_t*)calloc(1, sizeof(pu_data_t));
- d->max_depth = 0;
- d->tid = -1; d->mask = BAM_DEF_MASK;
- d->c = bam_maqcns_init();
- d->c->is_soap = 1; // change the default model
- d->ido = bam_maqindel_opt_init();
- while ((c = getopt(argc, argv, "st:f:cT:N:r:l:d:im:gI:G:vM:S2aR:PA")) >= 0) {
- switch (c) {
- case 'a': d->c->is_soap = 1; break;
- case 'A': d->c->is_soap = 0; break;
- case 's': d->format |= BAM_PLF_SIMPLE; break;
- case 't': fn_list = strdup(optarg); break;
- case 'l': fn_pos = strdup(optarg); break;
- case 'f': fn_fa = strdup(optarg); break;
- case 'T': d->c->theta = atof(optarg); break;
- case 'N': d->c->n_hap = atoi(optarg); break;
- case 'r': d->c->het_rate = atof(optarg); d->ido->r_snp = d->c->het_rate; break;
- case 'M': d->c->cap_mapQ = atoi(optarg); break;
- case 'd': d->max_depth = atoi(optarg); break;
- case 'c': d->format |= BAM_PLF_CNS; break;
- case 'i': d->format |= BAM_PLF_INDEL_ONLY; break;
- case 'v': d->format |= BAM_PLF_VAR_ONLY; break;
- case 'm': d->mask = strtol(optarg, 0, 0); break;
- case 'g': d->format |= BAM_PLF_GLF; break;
- case '2': d->format |= BAM_PLF_2ND; break;
- case 'P': d->format |= BAM_PLF_READPOS; break;
- case 'I': d->ido->q_indel = atoi(optarg); break;
- case 'G': d->ido->r_indel = atof(optarg); break;
- case 'S': is_SAM = 1; break;
- case 'R':
- if (strcmp(optarg, "random") == 0) d->format |= BAM_PLF_RANBASE;
- else if (strcmp(optarg, "first") == 0) d->format |= BAM_PLF_1STBASE;
- else if (strcmp(optarg, "all") == 0) d->format |= BAM_PLF_ALLBASE;
- else fprintf(stderr, "[bam_pileup] unrecognized -R\n");
- break;
- default: fprintf(stderr, "Unrecognizd option '-%c'.\n", c); return 1;