X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=misc%2Fbamcheck.c;h=c7db8abfb8c153fcc5a95c947a6013e258e6e95a;hb=1ef8ee37fa796dd58e4c18681939a2fafe84ee0c;hp=03cf3b233fec3b7d976dc07ac82fad94fd47a70c;hpb=8f9081fe5d62916987b928753d13e06c08322236;p=samtools.git diff --git a/misc/bamcheck.c b/misc/bamcheck.c index 03cf3b2..c7db8ab 100644 --- a/misc/bamcheck.c +++ b/misc/bamcheck.c @@ -8,14 +8,16 @@ - coverage distribution ignores softclips and deletions - some stats require sorted BAMs - GC content graph can have an untidy, step-like pattern when BAM contains multiple read lengths. - - The whole reads are used with -t, no splicing is done, no indels or soft clips are - considered, even small overlap is good enough to include the read in the stats. + - 'bases mapped' (stats->nbases_mapped) is calculated from read lengths given by BAM (core.l_qseq) + - With the -t option, the whole reads are used. Except for the number of mapped bases (cigar) + counts, no splicing is done, no indels or soft clips are considered, even small overlap is + good enough to include the read in the stats. + */ -#define BAMCHECK_VERSION "2012-03-29" +#define BAMCHECK_VERSION "2012-04-24" #define _ISOC99_SOURCE -#define _GNU_SOURCE #include #include #include @@ -146,7 +148,7 @@ typedef struct int filter_readlen; // Target regions - int nregions; + int nregions, reg_from,reg_to; regions_t *regions; // Auxiliary data @@ -159,6 +161,9 @@ typedef struct stats_t; void error(const char *format, ...); +void bam_init_header_hash(bam_header_t *header); +int is_in_regions(bam1_t *bam_line, stats_t *stats); + // Coverage distribution methods inline int coverage_idx(int min, int max, int n, int step, int depth) @@ -282,8 +287,9 @@ void count_indels(stats_t *stats,bam1_t *bam_line) if ( cig==1 ) { - int idx = is_fwd ? icycle : read_len-icycle-1; - if ( idx >= stats->nbases ) error("FIXME: %d vs %d\n", idx,stats->nbases); + int idx = is_fwd ? icycle : read_len-icycle; + if ( idx<0 ) error("FIXME: read_len=%d vs icycle=%d\n", read_len,icycle); + if ( idx >= stats->nbases || idx<0 ) error("FIXME: %d vs %d\n", idx,stats->nbases); stats->ins_cycles[idx]++; icycle += ncig; if ( ncig<=stats->nindels ) @@ -292,7 +298,8 @@ void count_indels(stats_t *stats,bam1_t *bam_line) } if ( cig==2 ) { - int idx = is_fwd ? icycle : read_len-icycle-1; + int idx = is_fwd ? icycle-1 : read_len-icycle-1; + if ( idx<0 ) continue; // discard meaningless deletions if ( idx >= stats->nbases ) error("FIXME: %d vs %d\n", idx,stats->nbases); stats->del_cycles[idx]++; if ( ncig<=stats->nindels ) @@ -515,15 +522,15 @@ void realloc_buffers(stats_t *stats, int seq_len) error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*sizeof(uint64_t)); memset(stats->deletions + stats->nbases, 0, (n-stats->nbases)*sizeof(uint64_t)); - stats->ins_cycles = realloc(stats->ins_cycles, n*sizeof(uint64_t)); + stats->ins_cycles = realloc(stats->ins_cycles, (n+1)*sizeof(uint64_t)); if ( !stats->ins_cycles ) - error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*sizeof(uint64_t)); - memset(stats->ins_cycles + stats->nbases, 0, (n-stats->nbases)*sizeof(uint64_t)); + error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,(n+1)*sizeof(uint64_t)); + memset(stats->ins_cycles + stats->nbases + 1, 0, (n+1-stats->nbases)*sizeof(uint64_t)); - stats->del_cycles = realloc(stats->del_cycles, n*sizeof(uint64_t)); + stats->del_cycles = realloc(stats->del_cycles, (n+1)*sizeof(uint64_t)); if ( !stats->del_cycles ) - error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*sizeof(uint64_t)); - memset(stats->del_cycles + stats->nbases, 0, (n-stats->nbases)*sizeof(uint64_t)); + error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,(n+1)*sizeof(uint64_t)); + memset(stats->del_cycles + stats->nbases + 1, 0, (n+1-stats->nbases)*sizeof(uint64_t)); stats->nbases = n; @@ -544,6 +551,9 @@ void collect_stats(bam1_t *bam_line, stats_t *stats) if ( stats->rmdup && IS_DUP(bam_line) ) return; + if ( !is_in_regions(bam_line,stats) ) + return; + int seq_len = bam_line->core.l_qseq; if ( !seq_len ) return; if ( stats->filter_readlen!=-1 && seq_len!=stats->filter_readlen ) return; @@ -665,19 +675,47 @@ void collect_stats(bam1_t *bam_line, stats_t *stats) stats->nmismatches += bam_aux2i(nm); // Number of mapped bases from cigar + // Conversion from uint32_t to MIDNSHP + // 012-4-- + // MIDNSHP if ( bam_line->core.n_cigar == 0) error("FIXME: mapped read with no cigar?\n"); - int readlen = seq_len; - for (i=0; icore.n_cigar; i++) + int readlen=seq_len; + if ( stats->regions ) + { + // Count only on-target bases + int iref = bam_line->core.pos + 1; + for (i=0; icore.n_cigar; i++) + { + int cig = bam1_cigar(bam_line)[i]&BAM_CIGAR_MASK; + int ncig = bam1_cigar(bam_line)[i]>>BAM_CIGAR_SHIFT; + if ( cig==2 ) readlen += ncig; + else if ( cig==0 ) + { + if ( iref < stats->reg_from ) ncig -= stats->reg_from-iref; + else if ( iref+ncig-1 > stats->reg_to ) ncig -= iref+ncig-1 - stats->reg_to; + if ( ncig<0 ) ncig = 0; + stats->nbases_mapped_cigar += ncig; + iref += bam1_cigar(bam_line)[i]>>BAM_CIGAR_SHIFT; + } + else if ( cig==1 ) + { + iref += ncig; + if ( iref>=stats->reg_from && iref<=stats->reg_to ) + stats->nbases_mapped_cigar += ncig; + } + } + } + else { - // Conversion from uint32_t to MIDNSHP - // 01--4-- - // MIDNSHP - if ( (bam1_cigar(bam_line)[i]&BAM_CIGAR_MASK)==0 || (bam1_cigar(bam_line)[i]&BAM_CIGAR_MASK)==1 ) - stats->nbases_mapped_cigar += bam1_cigar(bam_line)[i]>>BAM_CIGAR_SHIFT; - - if ( (bam1_cigar(bam_line)[i]&BAM_CIGAR_MASK)==2 ) - readlen += bam1_cigar(bam_line)[i]>>BAM_CIGAR_SHIFT; + // Count the whole read + for (i=0; icore.n_cigar; i++) + { + if ( (bam1_cigar(bam_line)[i]&BAM_CIGAR_MASK)==0 || (bam1_cigar(bam_line)[i]&BAM_CIGAR_MASK)==1 ) + stats->nbases_mapped_cigar += bam1_cigar(bam_line)[i]>>BAM_CIGAR_SHIFT; + if ( (bam1_cigar(bam_line)[i]&BAM_CIGAR_MASK)==2 ) + readlen += bam1_cigar(bam_line)[i]>>BAM_CIGAR_SHIFT; + } } stats->nbases_mapped += seq_len; @@ -819,23 +857,23 @@ void output_stats(stats_t *stats) printf(" %s",stats->argv[i]); printf("\n"); printf("# Summary Numbers. Use `grep ^SN | cut -f 2-` to extract this part.\n"); - printf("SN\tsequences:\t%ld\n", stats->nreads_1st+stats->nreads_2nd); + printf("SN\tsequences:\t%ld\n", (long)(stats->nreads_1st+stats->nreads_2nd)); printf("SN\tis paired:\t%d\n", stats->nreads_1st&&stats->nreads_2nd ? 1 : 0); printf("SN\tis sorted:\t%d\n", stats->is_sorted ? 1 : 0); - printf("SN\t1st fragments:\t%ld\n", stats->nreads_1st); - printf("SN\tlast fragments:\t%ld\n", stats->nreads_2nd); - printf("SN\treads mapped:\t%ld\n", stats->nreads_paired+stats->nreads_unpaired); - printf("SN\treads unmapped:\t%ld\n", stats->nreads_unmapped); - printf("SN\treads unpaired:\t%ld\n", stats->nreads_unpaired); - printf("SN\treads paired:\t%ld\n", stats->nreads_paired); - printf("SN\treads duplicated:\t%ld\n", stats->nreads_dup); - printf("SN\treads MQ0:\t%ld\n", stats->nreads_mq0); - printf("SN\ttotal length:\t%ld\n", stats->total_len); - printf("SN\tbases mapped:\t%ld\n", stats->nbases_mapped); - printf("SN\tbases mapped (cigar):\t%ld\n", stats->nbases_mapped_cigar); - printf("SN\tbases trimmed:\t%ld\n", stats->nbases_trimmed); - printf("SN\tbases duplicated:\t%ld\n", stats->total_len_dup); - printf("SN\tmismatches:\t%ld\n", stats->nmismatches); + printf("SN\t1st fragments:\t%ld\n", (long)stats->nreads_1st); + printf("SN\tlast fragments:\t%ld\n", (long)stats->nreads_2nd); + printf("SN\treads mapped:\t%ld\n", (long)(stats->nreads_paired+stats->nreads_unpaired)); + printf("SN\treads unmapped:\t%ld\n", (long)stats->nreads_unmapped); + printf("SN\treads unpaired:\t%ld\n", (long)stats->nreads_unpaired); + printf("SN\treads paired:\t%ld\n", (long)stats->nreads_paired); + printf("SN\treads duplicated:\t%ld\n", (long)stats->nreads_dup); + printf("SN\treads MQ0:\t%ld\n", (long)stats->nreads_mq0); + printf("SN\ttotal length:\t%ld\n", (long)stats->total_len); + printf("SN\tbases mapped:\t%ld\n", (long)stats->nbases_mapped); + printf("SN\tbases mapped (cigar):\t%ld\n", (long)stats->nbases_mapped_cigar); + printf("SN\tbases trimmed:\t%ld\n", (long)stats->nbases_trimmed); + printf("SN\tbases duplicated:\t%ld\n", (long)stats->total_len_dup); + printf("SN\tmismatches:\t%ld\n", (long)stats->nmismatches); printf("SN\terror rate:\t%e\n", (float)stats->nmismatches/stats->nbases_mapped_cigar); float avg_read_length = (stats->nreads_1st+stats->nreads_2nd)?stats->total_len/(stats->nreads_1st+stats->nreads_2nd):0; printf("SN\taverage length:\t%.0f\n", avg_read_length); @@ -843,9 +881,9 @@ void output_stats(stats_t *stats) printf("SN\taverage quality:\t%.1f\n", stats->total_len?stats->sum_qual/stats->total_len:0); printf("SN\tinsert size average:\t%.1f\n", avg_isize); printf("SN\tinsert size standard deviation:\t%.1f\n", sd_isize); - printf("SN\tinward oriented pairs:\t%ld\n", nisize_inward); - printf("SN\toutward oriented pairs:\t%ld\n", nisize_outward); - printf("SN\tpairs with other orientation:\t%ld\n", nisize_other); + printf("SN\tinward oriented pairs:\t%ld\n", (long)nisize_inward); + printf("SN\toutward oriented pairs:\t%ld\n", (long)nisize_outward); + printf("SN\tpairs with other orientation:\t%ld\n", (long)nisize_other); int ibase,iqual; if ( stats->max_lennbases ) stats->max_len++; @@ -857,7 +895,7 @@ void output_stats(stats_t *stats) printf("FFQ\t%d",ibase+1); for (iqual=0; iqual<=stats->max_qual; iqual++) { - printf("\t%ld", stats->quals_1st[ibase*stats->nquals+iqual]); + printf("\t%ld", (long)stats->quals_1st[ibase*stats->nquals+iqual]); } printf("\n"); } @@ -868,7 +906,7 @@ void output_stats(stats_t *stats) printf("LFQ\t%d",ibase+1); for (iqual=0; iqual<=stats->max_qual; iqual++) { - printf("\t%ld", stats->quals_2nd[ibase*stats->nquals+iqual]); + printf("\t%ld", (long)stats->quals_2nd[ibase*stats->nquals+iqual]); } printf("\n"); } @@ -882,7 +920,7 @@ void output_stats(stats_t *stats) printf("MPC\t%d",ibase+1); for (iqual=0; iqual<=stats->max_qual; iqual++) { - printf("\t%ld", stats->mpc_buf[ibase*stats->nquals+iqual]); + printf("\t%ld", (long)stats->mpc_buf[ibase*stats->nquals+iqual]); } printf("\n"); } @@ -892,7 +930,7 @@ void output_stats(stats_t *stats) for (ibase=0; ibasengc; ibase++) { if ( stats->gc_1st[ibase]==stats->gc_1st[ibase_prev] ) continue; - printf("GCF\t%.2f\t%ld\n", (ibase+ibase_prev)*0.5*100./(stats->ngc-1),stats->gc_1st[ibase_prev]); + printf("GCF\t%.2f\t%ld\n", (ibase+ibase_prev)*0.5*100./(stats->ngc-1), (long)stats->gc_1st[ibase_prev]); ibase_prev = ibase; } printf("# GC Content of last fragments. Use `grep ^GCL | cut -f 2-` to extract this part.\n"); @@ -900,7 +938,7 @@ void output_stats(stats_t *stats) for (ibase=0; ibasengc; ibase++) { if ( stats->gc_2nd[ibase]==stats->gc_2nd[ibase_prev] ) continue; - printf("GCL\t%.2f\t%ld\n", (ibase+ibase_prev)*0.5*100./(stats->ngc-1),stats->gc_2nd[ibase_prev]); + printf("GCL\t%.2f\t%ld\n", (ibase+ibase_prev)*0.5*100./(stats->ngc-1), (long)stats->gc_2nd[ibase_prev]); ibase_prev = ibase; } printf("# ACGT content per cycle. Use `grep ^GCC | cut -f 2-` to extract this part. The columns are: cycle, and A,C,G,T counts [%%]\n"); @@ -913,37 +951,37 @@ void output_stats(stats_t *stats) } printf("# Insert sizes. Use `grep ^IS | cut -f 2-` to extract this part. The columns are: pairs total, inward oriented pairs, outward oriented pairs, other pairs\n"); for (isize=1; isizeisize_inward[isize]+stats->isize_outward[isize]+stats->isize_other[isize]), - stats->isize_inward[isize],stats->isize_outward[isize],stats->isize_other[isize]); + printf("IS\t%d\t%ld\t%ld\t%ld\t%ld\n", isize, (long)(stats->isize_inward[isize]+stats->isize_outward[isize]+stats->isize_other[isize]), + (long)stats->isize_inward[isize], (long)stats->isize_outward[isize], (long)stats->isize_other[isize]); printf("# Read lengths. Use `grep ^RL | cut -f 2-` to extract this part. The columns are: read length, count\n"); int ilen; for (ilen=0; ilenmax_len; ilen++) { if ( stats->read_lengths[ilen]>0 ) - printf("RL\t%d\t%ld\n", ilen,stats->read_lengths[ilen]); + printf("RL\t%d\t%ld\n", ilen, (long)stats->read_lengths[ilen]); } printf("# Indel distribution. Use `grep ^ID | cut -f 2-` to extract this part. The columns are: length, number of insertions, number of deletions\n"); for (ilen=0; ilennindels; ilen++) { if ( stats->insertions[ilen]>0 || stats->deletions[ilen]>0 ) - printf("ID\t%d\t%ld\t%ld\n", ilen+1,stats->insertions[ilen],stats->deletions[ilen]); + printf("ID\t%d\t%ld\t%ld\n", ilen+1, (long)stats->insertions[ilen], (long)stats->deletions[ilen]); } printf("# Indels per cycle. Use `grep ^IC | cut -f 2-` to extract this part. The columns are: cycle, number of insertions, number of deletions\n"); - for (ilen=0; ilennbases; ilen++) + for (ilen=0; ilen<=stats->nbases; ilen++) { if ( stats->ins_cycles[ilen]>0 || stats->del_cycles[ilen]>0 ) - printf("IC\t%d\t%ld\t%ld\n", ilen+1,stats->ins_cycles[ilen],stats->del_cycles[ilen]); + printf("IC\t%d\t%ld\t%ld\n", ilen+1, (long)stats->ins_cycles[ilen], (long)stats->del_cycles[ilen]); } printf("# Coverage distribution. Use `grep ^COV | cut -f 2-` to extract this part.\n"); - printf("COV\t[<%d]\t%d\t%ld\n",stats->cov_min,stats->cov_min-1,stats->cov[0]); + printf("COV\t[<%d]\t%d\t%ld\n",stats->cov_min,stats->cov_min-1, (long)stats->cov[0]); int icov; for (icov=1; icovncov-1; icov++) - printf("COV\t[%d-%d]\t%d\t%ld\n",stats->cov_min + (icov-1)*stats->cov_step, stats->cov_min + icov*stats->cov_step-1,stats->cov_min + icov*stats->cov_step-1,stats->cov[icov]); - printf("COV\t[%d<]\t%d\t%ld\n",stats->cov_min + (stats->ncov-2)*stats->cov_step-1,stats->cov_min + (stats->ncov-2)*stats->cov_step-1,stats->cov[stats->ncov-1]); + printf("COV\t[%d-%d]\t%d\t%ld\n",stats->cov_min + (icov-1)*stats->cov_step, stats->cov_min + icov*stats->cov_step-1,stats->cov_min + icov*stats->cov_step-1, (long)stats->cov[icov]); + printf("COV\t[%d<]\t%d\t%ld\n",stats->cov_min + (stats->ncov-2)*stats->cov_step-1,stats->cov_min + (stats->ncov-2)*stats->cov_step-1, (long)stats->cov[stats->ncov-1]); // Calculate average GC content, then sort by GC and depth @@ -980,7 +1018,39 @@ void output_stats(stats_t *stats) } } -void bam_init_header_hash(bam_header_t *header); +size_t mygetline(char **line, size_t *n, FILE *fp) +{ + if (line == NULL || n == NULL || fp == NULL) + { + errno = EINVAL; + return -1; + } + if (*n==0 || !*line) + { + *line = NULL; + *n = 0; + } + + size_t nread=0; + int c; + while ((c=getc(fp))!= EOF && c!='\n') + { + if ( ++nread>=*n ) + { + *n += 255; + *line = realloc(*line, sizeof(char)*(*n)); + } + (*line)[nread-1] = c; + } + if ( nread>=*n ) + { + *n += 255; + *line = realloc(*line, sizeof(char)*(*n)); + } + (*line)[nread] = 0; + return nread>0 ? nread : -1; + +} void init_regions(stats_t *stats, char *file) { @@ -998,13 +1068,13 @@ void init_regions(stats_t *stats, char *file) ssize_t nread; int warned = 0; int prev_tid=-1, prev_pos=-1; - while ((nread = getline(&line, &len, fp)) != -1) + while ((nread = mygetline(&line, &len, fp)) != -1) { if ( line[0] == '#' ) continue; int i = 0; while ( i=nread ) error("Could not parse the file: %s\n", file); + if ( i>=nread ) error("Could not parse the file: %s [%s]\n", file,line); line[i] = 0; iter = kh_get(str, header_hash, line); @@ -1012,7 +1082,7 @@ void init_regions(stats_t *stats, char *file) if ( iter == kh_end(header_hash) ) { if ( !warned ) - fprintf(stderr,"Warning: Some sequences not present in the BAM (%s)\n", line); + fprintf(stderr,"Warning: Some sequences not present in the BAM, e.g. \"%s\". This message is printed only once.\n", line); warned = 1; continue; } @@ -1046,6 +1116,7 @@ void init_regions(stats_t *stats, char *file) stats->regions[tid].npos++; } if (line) free(line); + if ( !stats->regions ) error("Unable to map the -t sequences to the BAM sequences.\n"); fclose(fp); } @@ -1066,6 +1137,35 @@ static int fetch_read(const bam1_t *bam_line, void *data) return 1; } +void reset_regions(stats_t *stats) +{ + int i; + for (i=0; inregions; i++) + stats->regions[i].cpos = 0; +} + +int is_in_regions(bam1_t *bam_line, stats_t *stats) +{ + if ( !stats->regions ) return 1; + + if ( bam_line->core.tid >= stats->nregions || bam_line->core.tid<0 ) return 0; + if ( !stats->is_sorted ) error("The BAM must be sorted in order for -t to work.\n"); + + regions_t *reg = &stats->regions[bam_line->core.tid]; + if ( reg->cpos==reg->npos ) return 0; // done for this chr + + // Find a matching interval or skip this read. No splicing of reads is done, no indels or soft clips considered, + // even small overlap is enough to include the read in the stats. + int i = reg->cpos; + while ( inpos && reg->pos[i].to<=bam_line->core.pos ) i++; + if ( i>=reg->npos ) { reg->cpos = reg->npos; return 0; } + if ( bam_line->core.pos + bam_line->core.l_qseq + 1 < reg->pos[i].from ) return 0; + reg->cpos = i; + stats->reg_from = reg->pos[i].from; + stats->reg_to = reg->pos[i].to; + + return 1; +} void error(const char *format, ...) { @@ -1210,8 +1310,8 @@ int main(int argc, char *argv[]) stats->read_lengths = calloc(stats->nbases,sizeof(uint64_t)); stats->insertions = calloc(stats->nbases,sizeof(uint64_t)); stats->deletions = calloc(stats->nbases,sizeof(uint64_t)); - stats->ins_cycles = calloc(stats->nbases,sizeof(uint64_t)); - stats->del_cycles = calloc(stats->nbases,sizeof(uint64_t)); + stats->ins_cycles = calloc(stats->nbases+1,sizeof(uint64_t)); + stats->del_cycles = calloc(stats->nbases+1,sizeof(uint64_t)); if ( targets ) init_regions(stats, targets); @@ -1229,6 +1329,7 @@ int main(int argc, char *argv[]) int tid, beg, end; bam_parse_region(stats->sam->header, argv[i], &tid, &beg, &end); if ( tid < 0 ) continue; + reset_regions(stats); bam_fetch(stats->sam->x.bam, bam_idx, tid, beg, end, stats, fetch_read); } bam_index_destroy(bam_idx); @@ -1237,25 +1338,7 @@ int main(int argc, char *argv[]) { // Stream through the entire BAM ignoring off-target regions if -t is given while (samread(sam,bam_line) >= 0) - { - if ( stats->regions ) - { - if ( bam_line->core.tid >= stats->nregions || bam_line->core.tid<0 ) continue; - if ( !stats->is_sorted ) error("The BAM must be sorted in order for -t to work.\n"); - - regions_t *reg = &stats->regions[bam_line->core.tid]; - if ( reg->cpos==reg->npos ) continue; // done for this chr - - // Find a matching interval or skip this read. No splicing of reads is done, no indels or soft clips considered, - // even small overlap is enough to include the read in the stats. - int i = reg->cpos; - while ( inpos && reg->pos[i].to<=bam_line->core.pos ) i++; - if ( i>=reg->npos ) { reg->cpos = reg->npos; continue; } - if ( bam_line->core.pos + bam_line->core.l_qseq + 1 < reg->pos[i].from ) continue; - reg->cpos = i; - } collect_stats(bam_line,stats); - } } round_buffer_flush(stats,-1);