X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bam_plcmd.c;h=e2601ecc3c26b41051d0148a9cf9696036d297fb;hb=a7dfcd368cda94d5847d3f2b42d4288399d77835;hp=96f5c459df9dd88aa35d69a361898739016875d8;hpb=9360df156b23158232dca10ffdb8a8ffdcaf0930;p=samtools.git diff --git a/bam_plcmd.c b/bam_plcmd.c index 96f5c45..e2601ec 100644 --- a/bam_plcmd.c +++ b/bam_plcmd.c @@ -17,6 +17,7 @@ KHASH_MAP_INIT_INT64(64, indel_list_t) #define BAM_PLF_INDEL_ONLY 0x04 #define BAM_PLF_GLF 0x08 #define BAM_PLF_VAR_ONLY 0x10 +#define BAM_PLF_2ND 0x20 typedef struct { bam_header_t *h; @@ -26,7 +27,7 @@ typedef struct { khash_t(64) *hash; uint32_t format; int tid, len, last_pos; - int mask, min_mapQ; + int mask; char *ref; glfFile fp_glf; // for glf output only } pu_data_t; @@ -181,7 +182,7 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p // call the consensus and indel if (d->format & BAM_PLF_CNS) // call consensus cns = bam_maqcns_call(n, pu, d->c); - if ((d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY)) && d->ref) { // call indels + if ((d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY)) && d->ref && pos < d->len) { // call indels if (proposed_indels) // the first element gives the size of the array r = bam_maqindel(n, pos, d->ido, pu, d->ref, proposed_indels[0], proposed_indels+1); else r = bam_maqindel(n, pos, d->ido, pu, d->ref, 0, 0); @@ -189,8 +190,10 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p // when only variant sites are asked for, test if the site is a variant if ((d->format & BAM_PLF_CNS) && (d->format & BAM_PLF_VAR_ONLY)) { if (!(bam_nt16_table[rb] != 15 && cns>>28 != bam_nt16_table[rb])) { // not a SNP - if (!(r && (r->gt == 2 || strcmp(r->s[r->gt], "*")))) // not an indel + if (!(r && (r->gt == 2 || strcmp(r->s[r->gt], "*")))) { // not an indel + if (r) bam_maqindel_ret_destroy(r); return 0; + } } } // print the first 3 columns @@ -246,6 +249,21 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p if (c > 126) c = 126; putchar(c); } + 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'); @@ -266,7 +284,8 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p 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\n", r->cnt1, r->cnt2, r->cnt_anti); + 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; @@ -274,14 +293,15 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p int bam_pileup(int argc, char *argv[]) { - int c; + 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->tid = -1; d->mask = BAM_DEF_MASK; d->c = bam_maqcns_init(); d->ido = bam_maqindel_opt_init(); - while ((c = getopt(argc, argv, "st:f:cT:N:r:l:im:gI:G:vM:q:")) >= 0) { + while ((c = getopt(argc, argv, "st:f:cT:N:r:l:im:gI:G:vM:S2a")) >= 0) { switch (c) { + case 'a': d->c->is_soap = 1; break; case 's': d->format |= BAM_PLF_SIMPLE; break; case 't': fn_list = strdup(optarg); break; case 'l': fn_pos = strdup(optarg); break; @@ -295,21 +315,25 @@ int bam_pileup(int argc, char *argv[]) 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 'I': d->ido->q_indel = atoi(optarg); break; case 'G': d->ido->r_indel = atof(optarg); break; - case 'q': d->min_mapQ = atoi(optarg); break; + case 'S': is_SAM = 1; break; default: fprintf(stderr, "Unrecognizd option '-%c'.\n", c); return 1; } } + if (fn_list) is_SAM = 1; if (optind == argc) { fprintf(stderr, "\n"); fprintf(stderr, "Usage: samtools pileup [options] |\n\n"); fprintf(stderr, "Option: -s simple (yet incomplete) pileup format\n"); + fprintf(stderr, " -S the input is in SAM\n"); + fprintf(stderr, " -a use the SOAPsnp model for SNP calling\n"); + fprintf(stderr, " -2 output the 2nd best call and quality\n"); fprintf(stderr, " -i only show lines/consensus with indels\n"); fprintf(stderr, " -m INT filtering reads with bits in INT [%d]\n", d->mask); - fprintf(stderr, " -q INT filtering reads with mapping quality smaller than INT [%d]\n", d->min_mapQ); fprintf(stderr, " -M INT cap mapping quality at INT [%d]\n", d->c->cap_mapQ); - fprintf(stderr, " -t FILE list of reference sequences (assume the input is in SAM)\n"); + fprintf(stderr, " -t FILE list of reference sequences (force -S)\n"); fprintf(stderr, " -l FILE list of sites at which pileup is output\n"); fprintf(stderr, " -f FILE reference sequence in the FASTA format\n\n"); fprintf(stderr, " -c output the maq consensus sequence\n"); @@ -335,16 +359,18 @@ int bam_pileup(int argc, char *argv[]) } if (d->fai == 0 && (d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY))) fprintf(stderr, "[bam_pileup] indels will not be called when -f is absent.\n"); + if (fn_fa && is_SAM && fn_list == 0) fn_list = samfaipath(fn_fa); + { samfile_t *fp; - fp = fn_list? samopen(argv[optind], "r", fn_list) : samopen(argv[optind], "rb", 0); + fp = is_SAM? samopen(argv[optind], "r", fn_list) : samopen(argv[optind], "rb", 0); if (fp == 0 || fp->header == 0) { fprintf(stderr, "[bam_pileup] fail to read the header: non-exisiting file or wrong format.\n"); return 1; } d->h = fp->header; if (fn_pos) d->hash = load_pos(fn_pos, d->h); - sampileup(fp, d->mask, d->min_mapQ, pileup_func, d); + sampileup(fp, d->mask, pileup_func, d); samclose(fp); // d->h will be destroyed here }