X-Git-Url: https://git.donarmstrong.com/?p=samtools.git;a=blobdiff_plain;f=bam2depth.c;h=02311efc35db9155a5fc4d3174ad5e38506bea82;hp=95588f96b5a4ebc1a531bc11bb5c2c43183bf08d;hb=43e736c6325dd64edf0f5549a22b51380f7bba25;hpb=50706bea83d8a1518485283876333279f4ae7137 diff --git a/bam2depth.c b/bam2depth.c index 95588f9..02311ef 100644 --- a/bam2depth.c +++ b/bam2depth.c @@ -13,6 +13,7 @@ typedef struct { // auxiliary data structure bamFile fp; // the file handler bam_iter_t iter; // NULL if a region not specified + int min_mapQ, min_len; // mapQ filter; length filter } aux_t; void *bed_read(const char *fn); // read a BED or position list file @@ -20,50 +21,80 @@ void bed_destroy(void *_h); // destroy the BED data structure int bed_overlap(const void *_h, const char *chr, int beg, int end); // test if chr:beg-end overlaps // This function reads a BAM alignment from one BAM file. -static int read_bam(void *data, bam1_t *b) +static int read_bam(void *data, bam1_t *b) // read level filters better go here to avoid pileup { aux_t *aux = (aux_t*)data; // data in fact is a pointer to an auxiliary structure - return aux->iter? bam_iter_read(aux->fp, aux->iter, b) : bam_read1(aux->fp, b); + int ret = aux->iter? bam_iter_read(aux->fp, aux->iter, b) : bam_read1(aux->fp, b); + if (!(b->core.flag&BAM_FUNMAP)) { + if ((int)b->core.qual < aux->min_mapQ) b->core.flag |= BAM_FUNMAP; + else if (aux->min_len && bam_cigar2qlen(&b->core, bam1_cigar(b)) < aux->min_len) b->core.flag |= BAM_FUNMAP; + } + return ret; } +int read_file_list(const char *file_list,int *n,char **argv[]); + #ifdef _MAIN_BAM2DEPTH int main(int argc, char *argv[]) #else int main_depth(int argc, char *argv[]) #endif { - int i, n, tid, beg, end, pos, *n_plp, baseQ = 0; + int i, n, tid, beg, end, pos, *n_plp, baseQ = 0, mapQ = 0, min_len = 0, nfiles; const bam_pileup1_t **plp; char *reg = 0; // specified region void *bed = 0; // BED data structure + char *file_list = NULL, **fn = NULL; bam_header_t *h = 0; // BAM header of the 1st input aux_t **data; bam_mplp_t mplp; // parse the command line - while ((n = getopt(argc, argv, "r:b:q:")) >= 0) { + while ((n = getopt(argc, argv, "r:b:q:Q:l:f:")) >= 0) { switch (n) { + case 'l': min_len = atoi(optarg); break; // minimum query length case 'r': reg = strdup(optarg); break; // parsing a region requires a BAM header case 'b': bed = bed_read(optarg); break; // BED or position list file can be parsed now case 'q': baseQ = atoi(optarg); break; // base quality threshold + case 'Q': mapQ = atoi(optarg); break; // mapping quality threshold + case 'f': file_list = optarg; break; } } - if (optind == argc) { - fprintf(stderr, "Usage: bam2depth [-r reg] [-q baseQthres] [-b in.bed] [...]\n"); + if (optind == argc && !file_list) { + fprintf(stderr, "\n"); + fprintf(stderr, "Usage: samtools depth [options] in1.bam [in2.bam [...]]\n"); + fprintf(stderr, "Options:\n"); + fprintf(stderr, " -b list of positions or regions\n"); + fprintf(stderr, " -f list of input BAM filenames, one per line [null]\n"); + fprintf(stderr, " -l minQLen\n"); + fprintf(stderr, " -q base quality threshold\n"); + fprintf(stderr, " -Q mapping quality threshold\n"); + fprintf(stderr, " -r region\n"); + fprintf(stderr, "\n"); return 1; } // initialize the auxiliary data structures - n = argc - optind; // the number of BAMs on the command line - data = calloc(n, sizeof(void*)); - beg = 0; end = 1<<30; tid = -1; + if (file_list) + { + if ( read_file_list(file_list,&nfiles,&fn) ) return 1; + n = nfiles; + argv = fn; + optind = 0; + } + else + n = argc - optind; // the number of BAMs on the command line + data = calloc(n, sizeof(void*)); // data[i] for the i-th input + beg = 0; end = 1<<30; tid = -1; // set the default region for (i = 0; i < n; ++i) { bam_header_t *htmp; data[i] = calloc(1, sizeof(aux_t)); data[i]->fp = bam_open(argv[optind+i], "r"); // open BAM + data[i]->min_mapQ = mapQ; // set the mapQ filter + data[i]->min_len = min_len; // set the qlen filter htmp = bam_header_read(data[i]->fp); // read the BAM header if (i == 0) { - h = htmp; // keep the header for the 1st BAM + h = htmp; // keep the header of the 1st BAM if (reg) bam_parse_region(h, reg, &tid, &beg, &end); // also parse the region } else bam_header_destroy(htmp); // if not the 1st BAM, trash the header if (tid >= 0) { // if a region is specified and parsed successfully @@ -80,15 +111,15 @@ int main_depth(int argc, char *argv[]) while (bam_mplp_auto(mplp, &tid, &pos, n_plp, plp) > 0) { // come to the next covered position if (pos < beg || pos >= end) continue; // out of range; skip if (bed && bed_overlap(bed, h->target_name[tid], pos, pos + 1) == 0) continue; // not in BED; skip - fputs(h->target_name[tid], stdout); printf("\t%d", pos+1); - for (i = 0; i < n; ++i) { + fputs(h->target_name[tid], stdout); printf("\t%d", pos+1); // a customized printf() would be faster + for (i = 0; i < n; ++i) { // base level filters have to go here int j, m = 0; for (j = 0; j < n_plp[i]; ++j) { const bam_pileup1_t *p = plp[i] + j; // DON'T modfity plp[][] unless you really know if (p->is_del || p->is_refskip) ++m; // having dels or refskips at tid:pos else if (bam1_qual(p->b)[p->qpos] < baseQ) ++m; // low base quality } - printf("\t%d", n_plp[i] - m); + printf("\t%d", n_plp[i] - m); // this the depth to output } putchar('\n'); } @@ -103,5 +134,10 @@ int main_depth(int argc, char *argv[]) } free(data); free(reg); if (bed) bed_destroy(bed); + if ( file_list ) + { + for (i=0; i