new file mode 100644

@@ -0,0 +1,1521 @@

+/* 

+    Author: petr.danecek@sanger

+    gcc -Wall -Winline -g -O2 -I ~/git/samtools bamcheck.c -o bamcheck -lm -lz -L ~/git/samtools -lbam -lpthread

+

+    Assumptions, approximations and other issues:

+        - GC-depth graph does not split reads, the starting position determines which bin is incremented.

+            There are small overlaps between bins (max readlen-1). However, the bins are big (20k).

+        - 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.

+        - '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-09-04"

+

+#define _ISOC99_SOURCE

+#include <stdio.h>

+#include <stdlib.h>

+#include <stdarg.h>

+#include <string.h>

+#include <math.h>

+#include <ctype.h>

+#include <getopt.h>

+#include <errno.h>

+#include <assert.h>

+#include "faidx.h"

+#include "khash.h"

+#include "sam.h"

+#include "sam_header.h"

+#include "razf.h"

+

+#define BWA_MIN_RDLEN 35

+#define IS_PAIRED(bam) ((bam)->core.flag&BAM_FPAIRED && !((bam)->core.flag&BAM_FUNMAP) && !((bam)->core.flag&BAM_FMUNMAP))

+#define IS_UNMAPPED(bam) ((bam)->core.flag&BAM_FUNMAP)

+#define IS_REVERSE(bam) ((bam)->core.flag&BAM_FREVERSE)

+#define IS_MATE_REVERSE(bam) ((bam)->core.flag&BAM_FMREVERSE)

+#define IS_READ1(bam) ((bam)->core.flag&BAM_FREAD1)

+#define IS_READ2(bam) ((bam)->core.flag&BAM_FREAD2)

+#define IS_DUP(bam) ((bam)->core.flag&BAM_FDUP)

+

+typedef struct 

+{

+    int32_t line_len, line_blen;

+    int64_t len;

+    uint64_t offset;

+} 

+faidx1_t;

+KHASH_MAP_INIT_STR(kh_faidx, faidx1_t)

+KHASH_MAP_INIT_STR(kh_bam_tid, int)

+KHASH_MAP_INIT_STR(kh_rg, const char *)

+struct __faidx_t {

+    RAZF *rz;

+    int n, m;

+    char **name;

+    khash_t(kh_faidx) *hash;

+};

+

+typedef struct

+{

+    float gc;

+    uint32_t depth;

+}

+gc_depth_t;

+

+// For coverage distribution, a simple pileup

+typedef struct 

+{

+    int64_t pos;

+    int size, start;

+    int *buffer;

+} 

+round_buffer_t;

+

+typedef struct { uint32_t from, to; } pos_t;

+typedef struct

+{

+    int npos,mpos,cpos;

+    pos_t *pos;

+}

+regions_t;

+

+typedef struct

+{

+    // Parameters

+    int trim_qual;      // bwa trim quality

+

+    // Dimensions of the quality histogram holder (quals_1st,quals_2nd), GC content holder (gc_1st,gc_2nd),

+    //  insert size histogram holder

+    int nquals;         // The number of quality bins 

+    int nbases;         // The maximum sequence length the allocated array can hold

+    int nisize;         // The maximum insert size that the allocated array can hold

+    int ngc;            // The size of gc_1st and gc_2nd

+    int nindels;        // The maximum indel length for indel distribution

+

+    // Arrays for the histogram data

+    uint64_t *quals_1st, *quals_2nd;

+    uint64_t *gc_1st, *gc_2nd;

+    uint64_t *isize_inward, *isize_outward, *isize_other;

+    uint64_t *acgt_cycles;

+    uint64_t *read_lengths;

+    uint64_t *insertions, *deletions;

+    uint64_t *ins_cycles_1st, *ins_cycles_2nd, *del_cycles_1st, *del_cycles_2nd;

+

+    // The extremes encountered

+    int max_len;            // Maximum read length

+    int max_qual;           // Maximum quality

+    float isize_main_bulk;  // There are always some unrealistically big insert sizes, report only the main part

+    int is_sorted;

+

+    // Summary numbers

+    uint64_t total_len;

+    uint64_t total_len_dup;

+    uint64_t nreads_1st;

+    uint64_t nreads_2nd;

+    uint64_t nreads_filtered;

+    uint64_t nreads_dup;

+    uint64_t nreads_unmapped;

+    uint64_t nreads_unpaired;

+    uint64_t nreads_paired;

+    uint64_t nreads_anomalous;

+    uint64_t nreads_mq0;

+    uint64_t nbases_mapped;

+    uint64_t nbases_mapped_cigar;

+    uint64_t nbases_trimmed;  // bwa trimmed bases

+    uint64_t nmismatches;

+    uint64_t nreads_QCfailed, nreads_secondary;

+

+    // GC-depth related data

+    uint32_t ngcd, igcd;        // The maximum number of GC depth bins and index of the current bin

+    gc_depth_t *gcd;            // The GC-depth bins holder

+    int gcd_bin_size;           // The size of GC-depth bin

+    uint32_t gcd_ref_size;      // The approximate size of the genome

+    int32_t tid, gcd_pos;       // Position of the current bin

+    int32_t pos;                // Position of the last read

+

+    // Coverage distribution related data

+    int ncov;                       // The number of coverage bins

+    uint64_t *cov;                  // The coverage frequencies

+    int cov_min,cov_max,cov_step;   // Minimum, maximum coverage and size of the coverage bins

+    round_buffer_t cov_rbuf;        // Pileup round buffer

+

+    // Mismatches by read cycle 

+    uint8_t *rseq_buf;              // A buffer for reference sequence to check the mismatches against

+    int mrseq_buf;                  // The size of the buffer

+    int32_t rseq_pos;               // The coordinate of the first base in the buffer

+    int32_t nrseq_buf;              // The used part of the buffer

+    uint64_t *mpc_buf;              // Mismatches per cycle

+

+    // Filters

+    int filter_readlen;

+

+    // Target regions

+    int nregions, reg_from,reg_to;

+    regions_t *regions;

+

+    // Auxiliary data

+    int flag_require, flag_filter;

+    double sum_qual;                // For calculating average quality value 

+    samfile_t *sam;             

+    khash_t(kh_rg) *rg_hash;        // Read groups to include, the array is null-terminated

+    faidx_t *fai;                   // Reference sequence for GC-depth graph

+    int argc;                       // Command line arguments to be printed on the output

+    char **argv;

+}

+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)

+{

+    if ( depth < min )

+        return 0;

+

+    if ( depth > max )

+        return n-1;

+

+    return 1 + (depth - min) / step;

+}

+

+inline int round_buffer_lidx2ridx(int offset, int size, int64_t refpos, int64_t pos)

+{

+    return (offset + (pos-refpos) % size) % size;

+}

+

+void round_buffer_flush(stats_t *stats, int64_t pos)

+{

+    int ibuf,idp;

+

+    if ( pos==stats->cov_rbuf.pos ) 

+        return;

+

+    int64_t new_pos = pos;

+    if ( pos==-1 || pos - stats->cov_rbuf.pos >= stats->cov_rbuf.size )

+    {

+        // Flush the whole buffer, but in sequential order, 

+        pos = stats->cov_rbuf.pos + stats->cov_rbuf.size - 1;

+    }

+

+    if ( pos < stats->cov_rbuf.pos ) 

+        error("Expected coordinates in ascending order, got %ld after %ld\n", pos,stats->cov_rbuf.pos);

+

+    int ifrom = stats->cov_rbuf.start;

+    int ito = round_buffer_lidx2ridx(stats->cov_rbuf.start,stats->cov_rbuf.size,stats->cov_rbuf.pos,pos-1);

+    if ( ifrom>ito )

+    {

+        for (ibuf=ifrom; ibuf<stats->cov_rbuf.size; ibuf++)

+        {

+            if ( !stats->cov_rbuf.buffer[ibuf] )

+                continue;

+            idp = coverage_idx(stats->cov_min,stats->cov_max,stats->ncov,stats->cov_step,stats->cov_rbuf.buffer[ibuf]);

+            stats->cov[idp]++;

+            stats->cov_rbuf.buffer[ibuf] = 0;

+        }

+        ifrom = 0;

+    }

+    for (ibuf=ifrom; ibuf<=ito; ibuf++)

+    {

+        if ( !stats->cov_rbuf.buffer[ibuf] )

+            continue;

+        idp = coverage_idx(stats->cov_min,stats->cov_max,stats->ncov,stats->cov_step,stats->cov_rbuf.buffer[ibuf]);

+        stats->cov[idp]++;

+        stats->cov_rbuf.buffer[ibuf] = 0;

+    }

+    stats->cov_rbuf.start = (new_pos==-1) ? 0 : round_buffer_lidx2ridx(stats->cov_rbuf.start,stats->cov_rbuf.size,stats->cov_rbuf.pos,pos);

+    stats->cov_rbuf.pos   = new_pos;

+}

+

+void round_buffer_insert_read(round_buffer_t *rbuf, int64_t from, int64_t to)

+{

+    if ( to-from >= rbuf->size )

+        error("The read length too big (%d), please increase the buffer length (currently %d)\n", to-from+1,rbuf->size);

+    if ( from < rbuf->pos )

+        error("The reads are not sorted (%ld comes after %ld).\n", from,rbuf->pos);

+

+    int ifrom,ito,ibuf;

+    ifrom = round_buffer_lidx2ridx(rbuf->start,rbuf->size,rbuf->pos,from);

+    ito   = round_buffer_lidx2ridx(rbuf->start,rbuf->size,rbuf->pos,to);

+    if ( ifrom>ito )

+    {

+        for (ibuf=ifrom; ibuf<rbuf->size; ibuf++)

+            rbuf->buffer[ibuf]++;

+        ifrom = 0;

+    }

+    for (ibuf=ifrom; ibuf<=ito; ibuf++)

+        rbuf->buffer[ibuf]++;

+}

+

+// Calculate the number of bases in the read trimmed by BWA

+int bwa_trim_read(int trim_qual, uint8_t *quals, int len, int reverse) 

+{

+    if ( len<BWA_MIN_RDLEN ) return 0;

+

+    // Although the name implies that the read cannot be trimmed to more than BWA_MIN_RDLEN,

+    //  the calculation can in fact trim it to (BWA_MIN_RDLEN-1). (bwa_trim_read in bwa/bwaseqio.c).

+    int max_trimmed = len - BWA_MIN_RDLEN + 1;

+    int l, sum=0, max_sum=0, max_l=0;

+

+    for (l=0; l<max_trimmed; l++)

+    {

+        sum += trim_qual - quals[ reverse ? l : len-1-l ];

+        if ( sum<0 ) break;

+        if ( sum>max_sum )

+        {

+            max_sum = sum;

+            // This is the correct way, but bwa clips from some reason one base less

+            // max_l   = l+1;

+            max_l   = l;

+        }

+    }

+    return max_l;

+}

+

+

+void count_indels(stats_t *stats,bam1_t *bam_line) 

+{

+    int is_fwd = IS_REVERSE(bam_line) ? 0 : 1;

+    int is_1st = IS_READ1(bam_line) ? 1 : 0;

+    int icig;

+    int icycle = 0;

+    int read_len = bam_line->core.l_qseq;

+    for (icig=0; icig<bam_line->core.n_cigar; icig++) 

+    {

+        // Conversion from uint32_t to MIDNSHP

+        //  0123456

+        //  MIDNSHP

+        int cig  = bam1_cigar(bam_line)[icig] & BAM_CIGAR_MASK;

+        int ncig = bam1_cigar(bam_line)[icig] >> BAM_CIGAR_SHIFT;

+

+        if ( cig==1 )

+        {

+            int idx = is_fwd ? icycle : read_len-icycle-ncig;

+            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, %s:%d %s\n", idx,stats->nbases, stats->sam->header->target_name[bam_line->core.tid],bam_line->core.pos+1,bam1_qname(bam_line));

+            if ( is_1st ) 

+                stats->ins_cycles_1st[idx]++;

+            else

+                stats->ins_cycles_2nd[idx]++;

+            icycle += ncig;

+            if ( ncig<=stats->nindels )

+                stats->insertions[ncig-1]++;

+            continue;

+        }

+        if ( cig==2 )

+        {

+            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);

+            if ( is_1st )

+                stats->del_cycles_1st[idx]++;

+            else

+                stats->del_cycles_2nd[idx]++;

+            if ( ncig<=stats->nindels )

+                stats->deletions[ncig-1]++;

+            continue;

+        }

+        if ( cig!=3 && cig!=5 )

+            icycle += ncig;

+    }

+}

+

+void count_mismatches_per_cycle(stats_t *stats,bam1_t *bam_line) 

+{

+    int is_fwd = IS_REVERSE(bam_line) ? 0 : 1;

+    int icig,iread=0,icycle=0;

+    int iref = bam_line->core.pos - stats->rseq_pos;

+    int read_len   = bam_line->core.l_qseq;

+    uint8_t *read  = bam1_seq(bam_line);

+    uint8_t *quals = bam1_qual(bam_line);

+    uint64_t *mpc_buf = stats->mpc_buf;

+    for (icig=0; icig<bam_line->core.n_cigar; icig++) 

+    {

+        // Conversion from uint32_t to MIDNSHP

+        //  0123456

+        //  MIDNSHP

+        int cig  = bam1_cigar(bam_line)[icig] & BAM_CIGAR_MASK;

+        int ncig = bam1_cigar(bam_line)[icig] >> BAM_CIGAR_SHIFT;

+        if ( cig==1 )

+        {

+            iread  += ncig;

+            icycle += ncig;

+            continue;

+        }

+        if ( cig==2 )

+        {

+            iref += ncig;

+            continue;

+        }

+        if ( cig==4 )

+        {

+            icycle += ncig;

+            // Soft-clips are present in the sequence, but the position of the read marks a start of non-clipped sequence

+            //   iref += ncig;

+            iread  += ncig;

+            continue;

+        }

+        if ( cig==5 )

+        {

+            icycle += ncig;

+            continue;

+        }

+        // Ignore H and N CIGARs. The letter are inserted e.g. by TopHat and often require very large

+        //  chunk of refseq in memory. Not very frequent and not noticable in the stats.

+        if ( cig==3 || cig==5 ) continue;

+        if ( cig!=0 )

+            error("TODO: cigar %d, %s:%d %s\n", cig,stats->sam->header->target_name[bam_line->core.tid],bam_line->core.pos+1,bam1_qname(bam_line));

+       

+        if ( ncig+iref > stats->nrseq_buf )

+            error("FIXME: %d+%d > %d, %s, %s:%d\n",ncig,iref,stats->nrseq_buf, bam1_qname(bam_line),stats->sam->header->target_name[bam_line->core.tid],bam_line->core.pos+1);

+

+        int im;

+        for (im=0; im<ncig; im++)

+        {

+            uint8_t cread = bam1_seqi(read,iread);

+            uint8_t cref  = stats->rseq_buf[iref];

+

+            // ---------------15

+            // =ACMGRSVTWYHKDBN

+            if ( cread==15 )

+            {

+                int idx = is_fwd ? icycle : read_len-icycle-1;

+                if ( idx>stats->max_len )

+                    error("mpc: %d>%d\n",idx,stats->max_len);

+                idx = idx*stats->nquals;

+                if ( idx>=stats->nquals*stats->nbases )

+                    error("FIXME: mpc_buf overflow\n");

+                mpc_buf[idx]++;

+            }

+            else if ( cref && cread && cref!=cread )

+            {

+                uint8_t qual = quals[iread] + 1;

+                if ( qual>=stats->nquals )

+                    error("TODO: quality too high %d>=%d (%s %d %s)\n", qual,stats->nquals, stats->sam->header->target_name[bam_line->core.tid],bam_line->core.pos+1,bam1_qname(bam_line));

+

+                int idx = is_fwd ? icycle : read_len-icycle-1;

+                if ( idx>stats->max_len )

+                    error("mpc: %d>%d\n",idx,stats->max_len);

+

+                idx = idx*stats->nquals + qual;

+                if ( idx>=stats->nquals*stats->nbases )

+                    error("FIXME: mpc_buf overflow\n");

+                mpc_buf[idx]++;

+            }

+

+            iref++;

+            iread++;

+            icycle++;

+        }

+    }

+}

+

+void read_ref_seq(stats_t *stats,int32_t tid,int32_t pos)

+{

+    khash_t(kh_faidx) *h;

+    khiter_t iter;

+    faidx1_t val;

+    char *chr, c;

+    faidx_t *fai = stats->fai;

+

+    h = fai->hash;

+    chr = stats->sam->header->target_name[tid];

+

+    // ID of the sequence name

+    iter = kh_get(kh_faidx, h, chr);

+    if (iter == kh_end(h)) 

+        error("No such reference sequence [%s]?\n", chr);

+    val = kh_value(h, iter);

+

+    // Check the boundaries

+    if (pos >= val.len)

+        error("Was the bam file mapped with the reference sequence supplied?"

+              " A read mapped beyond the end of the chromosome (%s:%d, chromosome length %d).\n", chr,pos,val.len);

+    int size = stats->mrseq_buf;

+    // The buffer extends beyond the chromosome end. Later the rest will be filled with N's.

+    if (size+pos > val.len) size = val.len-pos;

+

+    // Position the razf reader

+    razf_seek(fai->rz, val.offset + pos / val.line_blen * val.line_len + pos % val.line_blen, SEEK_SET);

+

+    uint8_t *ptr = stats->rseq_buf;

+    int nread = 0;

+    while ( nread<size && razf_read(fai->rz,&c,1) && !fai->rz->z_err )

+    {

+        if ( !isgraph(c) )

+            continue;

+

+        // Conversion between uint8_t coding and ACGT

+        //      -12-4---8-------

+        //      =ACMGRSVTWYHKDBN

+        if ( c=='A' || c=='a' )

+            *ptr = 1;

+        else if ( c=='C' || c=='c' )

+            *ptr = 2;

+        else if ( c=='G' || c=='g' )

+            *ptr = 4;

+        else if ( c=='T' || c=='t' )

+            *ptr = 8;

+        else

+            *ptr = 0;

+        ptr++;

+        nread++;

+    }

+    if ( nread < stats->mrseq_buf )

+    {

+        memset(ptr,0, stats->mrseq_buf - nread);

+        nread = stats->mrseq_buf;

+    }

+    stats->nrseq_buf = nread;

+    stats->rseq_pos  = pos;

+    stats->tid       = tid;

+}

+

+float fai_gc_content(stats_t *stats, int pos, int len)

+{

+    uint32_t gc,count,c;

+    int i = pos - stats->rseq_pos, ito = i + len;

+    assert( i>=0 && ito<=stats->nrseq_buf );

+

+    // Count GC content

+    gc = count = 0;

+    for (; i<ito; i++)

+    {

+        c = stats->rseq_buf[i];

+        if ( c==2 || c==4 )

+        {

+            gc++;

+            count++;

+        }

+        else if ( c==1 || c==8 )

+            count++;

+    }

+    return count ? (float)gc/count : 0;

+}

+

+void realloc_rseq_buffer(stats_t *stats)

+{

+    int n = stats->nbases*10;

+    if ( stats->gcd_bin_size > n ) n = stats->gcd_bin_size;

+    if ( stats->mrseq_buf<n )

+    {

+        stats->rseq_buf = realloc(stats->rseq_buf,sizeof(uint8_t)*n);

+        stats->mrseq_buf = n;

+    }

+}

+

+void realloc_gcd_buffer(stats_t *stats, int seq_len)

+{

+    if ( seq_len >= stats->gcd_bin_size )

+        error("The --GC-depth bin size (%d) is set too low for the read length %d\n", stats->gcd_bin_size, seq_len);

+

+    int n = 1 + stats->gcd_ref_size / (stats->gcd_bin_size - seq_len);

+    if ( n <= stats->igcd )

+        error("The --GC-depth bin size is too small or reference genome too big; please decrease the bin size or increase the reference length\n");

+        

+    if ( n > stats->ngcd )

+    {

+        stats->gcd = realloc(stats->gcd, n*sizeof(gc_depth_t));

+        if ( !stats->gcd )

+            error("Could not realloc GCD buffer, too many chromosomes or the genome too long?? [%u %u]\n", stats->ngcd,n);

+        memset(&(stats->gcd[stats->ngcd]),0,(n-stats->ngcd)*sizeof(gc_depth_t));

+        stats->ngcd = n;

+    }

+

+    realloc_rseq_buffer(stats);

+}

+

+void realloc_buffers(stats_t *stats, int seq_len)

+{

+    int n = 2*(1 + seq_len - stats->nbases) + stats->nbases;

+

+    stats->quals_1st = realloc(stats->quals_1st, n*stats->nquals*sizeof(uint64_t));

+    if ( !stats->quals_1st )

+        error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*stats->nquals*sizeof(uint64_t));

+    memset(stats->quals_1st + stats->nbases*stats->nquals, 0, (n-stats->nbases)*stats->nquals*sizeof(uint64_t));

+

+    stats->quals_2nd = realloc(stats->quals_2nd, n*stats->nquals*sizeof(uint64_t));

+    if ( !stats->quals_2nd )

+        error("Could not realloc buffers, the sequence too long: %d (2x%ld)\n", seq_len,n*stats->nquals*sizeof(uint64_t));

+    memset(stats->quals_2nd + stats->nbases*stats->nquals, 0, (n-stats->nbases)*stats->nquals*sizeof(uint64_t));

+

+    if ( stats->mpc_buf )

+    {

+        stats->mpc_buf = realloc(stats->mpc_buf, n*stats->nquals*sizeof(uint64_t));

+        if ( !stats->mpc_buf )

+            error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*stats->nquals*sizeof(uint64_t));

+        memset(stats->mpc_buf + stats->nbases*stats->nquals, 0, (n-stats->nbases)*stats->nquals*sizeof(uint64_t));

+    }

+

+    stats->acgt_cycles = realloc(stats->acgt_cycles, n*4*sizeof(uint64_t));

+    if ( !stats->acgt_cycles )

+        error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*4*sizeof(uint64_t));

+    memset(stats->acgt_cycles + stats->nbases*4, 0, (n-stats->nbases)*4*sizeof(uint64_t));

+

+    stats->read_lengths = realloc(stats->read_lengths, n*sizeof(uint64_t));

+    if ( !stats->read_lengths )

+        error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*sizeof(uint64_t));

+    memset(stats->read_lengths + stats->nbases, 0, (n-stats->nbases)*sizeof(uint64_t));

+

+    stats->insertions = realloc(stats->insertions, n*sizeof(uint64_t));

+    if ( !stats->insertions )

+        error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,n*sizeof(uint64_t));

+    memset(stats->insertions + stats->nbases, 0, (n-stats->nbases)*sizeof(uint64_t));

+

+    stats->deletions = realloc(stats->deletions, n*sizeof(uint64_t));

+    if ( !stats->deletions )

+        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_1st = realloc(stats->ins_cycles_1st, (n+1)*sizeof(uint64_t));

+    if ( !stats->ins_cycles_1st )

+        error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,(n+1)*sizeof(uint64_t));

+    memset(stats->ins_cycles_1st + stats->nbases + 1, 0, (n-stats->nbases)*sizeof(uint64_t));

+

+    stats->ins_cycles_2nd = realloc(stats->ins_cycles_2nd, (n+1)*sizeof(uint64_t));

+    if ( !stats->ins_cycles_2nd )

+        error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,(n+1)*sizeof(uint64_t));

+    memset(stats->ins_cycles_2nd + stats->nbases + 1, 0, (n-stats->nbases)*sizeof(uint64_t));

+

+    stats->del_cycles_1st = realloc(stats->del_cycles_1st, (n+1)*sizeof(uint64_t));

+    if ( !stats->del_cycles_1st )

+        error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,(n+1)*sizeof(uint64_t));

+    memset(stats->del_cycles_1st + stats->nbases + 1, 0, (n-stats->nbases)*sizeof(uint64_t));

+

+    stats->del_cycles_2nd = realloc(stats->del_cycles_2nd, (n+1)*sizeof(uint64_t));

+    if ( !stats->del_cycles_2nd )

+        error("Could not realloc buffers, the sequence too long: %d (%ld)\n", seq_len,(n+1)*sizeof(uint64_t));

+    memset(stats->del_cycles_2nd + stats->nbases + 1, 0, (n-stats->nbases)*sizeof(uint64_t));

+

+    stats->nbases = n;

+

+    // Realloc the coverage distribution buffer

+    int *rbuffer = calloc(sizeof(int),seq_len*5);

+    n = stats->cov_rbuf.size-stats->cov_rbuf.start;

+    memcpy(rbuffer,stats->cov_rbuf.buffer+stats->cov_rbuf.start,n);

+    if ( stats->cov_rbuf.start>1 )

+        memcpy(rbuffer+n,stats->cov_rbuf.buffer,stats->cov_rbuf.start);

+    stats->cov_rbuf.start = 0;

+    free(stats->cov_rbuf.buffer);

+    stats->cov_rbuf.buffer = rbuffer;

+    stats->cov_rbuf.size = seq_len*5;

+

+    realloc_rseq_buffer(stats);

+}

+

+void collect_stats(bam1_t *bam_line, stats_t *stats)

+{

+    if ( stats->rg_hash )

+    {

+        const uint8_t *rg = bam_aux_get(bam_line, "RG");

+        if ( !rg ) return; 

+        khiter_t k = kh_get(kh_rg, stats->rg_hash, (const char*)(rg + 1));

+        if ( k == kh_end(stats->rg_hash) ) return;

+    }

+    if ( stats->flag_require && (bam_line->core.flag & stats->flag_require)!=stats->flag_require )

+    {

+        stats->nreads_filtered++;

+        return;

+    }

+    if ( stats->flag_filter && (bam_line->core.flag & stats->flag_filter) )

+    {

+        stats->nreads_filtered++;

+        return;

+    }

+    if ( !is_in_regions(bam_line,stats) )

+        return;

+    if ( stats->filter_readlen!=-1 && bam_line->core.l_qseq!=stats->filter_readlen ) 

+        return;

+

+    if ( bam_line->core.flag & BAM_FQCFAIL ) stats->nreads_QCfailed++;

+    if ( bam_line->core.flag & BAM_FSECONDARY ) stats->nreads_secondary++;

+

+    int seq_len = bam_line->core.l_qseq;

+    if ( !seq_len ) return;

+

+    if ( seq_len >= stats->nbases )

+        realloc_buffers(stats,seq_len);

+    if ( stats->max_len<seq_len )

+        stats->max_len = seq_len;

+

+    stats->read_lengths[seq_len]++;

+

+    // Count GC and ACGT per cycle

+    uint8_t base, *seq  = bam1_seq(bam_line);

+    int gc_count  = 0;

+    int i;

+    int reverse = IS_REVERSE(bam_line);

+    for (i=0; i<seq_len; i++)

+    {

+        // Conversion from uint8_t coding to ACGT

+        //      -12-4---8-------

+        //      =ACMGRSVTWYHKDBN

+        //       01 2   3

+        base = bam1_seqi(seq,i);

+        base /= 2;

+        if ( base==1 || base==2 ) gc_count++;

+        else if ( base>2 ) base=3;

+        if ( 4*(reverse ? seq_len-i-1 : i) + base >= stats->nbases*4 ) 

+            error("FIXME: acgt_cycles\n");

+        stats->acgt_cycles[ 4*(reverse ? seq_len-i-1 : i) + base ]++;

+    }

+    int gc_idx_min = gc_count*(stats->ngc-1)/seq_len;

+    int gc_idx_max = (gc_count+1)*(stats->ngc-1)/seq_len;

+    if ( gc_idx_max >= stats->ngc ) gc_idx_max = stats->ngc - 1;

+

+    // Determine which array (1st or 2nd read) will these stats go to,

+    //  trim low quality bases from end the same way BWA does, 

+    //  fill GC histogram

+    uint64_t *quals;

+    uint8_t *bam_quals = bam1_qual(bam_line);

+    if ( bam_line->core.flag&BAM_FREAD2 )

+    {

+        quals  = stats->quals_2nd;

+        stats->nreads_2nd++;

+        for (i=gc_idx_min; i<gc_idx_max; i++)

+            stats->gc_2nd[i]++;

+    }

+    else

+    {

+        quals = stats->quals_1st;

+        stats->nreads_1st++;

+        for (i=gc_idx_min; i<gc_idx_max; i++)

+            stats->gc_1st[i]++;

+    }

+    if ( stats->trim_qual>0 ) 

+        stats->nbases_trimmed += bwa_trim_read(stats->trim_qual, bam_quals, seq_len, reverse);

+

+    // Quality histogram and average quality

+    for (i=0; i<seq_len; i++)

+    {

+        uint8_t qual = bam_quals[ reverse ? seq_len-i-1 : i];

+        if ( qual>=stats->nquals )

+            error("TODO: quality too high %d>=%d (%s %d %s)\n", qual,stats->nquals,stats->sam->header->target_name[bam_line->core.tid],bam_line->core.pos+1,bam1_qname(bam_line));

+        if ( qual>stats->max_qual )

+            stats->max_qual = qual;

+

+        quals[ i*stats->nquals+qual ]++;

+        stats->sum_qual += qual;

+    }

+

+    // Look at the flags and increment appropriate counters (mapped, paired, etc)

+    if ( IS_UNMAPPED(bam_line) )

+        stats->nreads_unmapped++;

+    else

+    {

+        if ( !bam_line->core.qual )

+            stats->nreads_mq0++;

+

+        count_indels(stats,bam_line);

+

+        if ( !IS_PAIRED(bam_line) )

+            stats->nreads_unpaired++;

+        else 

+        {

+            stats->nreads_paired++;

+

+            if ( bam_line->core.tid!=bam_line->core.mtid )

+                stats->nreads_anomalous++;

+

+            // The insert size is tricky, because for long inserts the libraries are

+            // prepared differently and the pairs point in other direction. BWA does

+            // not set the paired flag for them.  Similar thing is true also for 454

+            // reads. Mates mapped to different chromosomes have isize==0.

+            int32_t isize = bam_line->core.isize;

+            if ( isize<0 ) isize = -isize;

+            if ( isize >= stats->nisize )

+                isize = stats->nisize-1;

+            if ( isize>0 || bam_line->core.tid==bam_line->core.mtid )

+            {

+                int pos_fst = bam_line->core.mpos - bam_line->core.pos;

+                int is_fst  = IS_READ1(bam_line) ? 1 : -1;

+                int is_fwd  = IS_REVERSE(bam_line) ? -1 : 1;

+                int is_mfwd = IS_MATE_REVERSE(bam_line) ? -1 : 1;

+

+                if ( is_fwd*is_mfwd>0 )

+                    stats->isize_other[isize]++;

+                else if ( is_fst*pos_fst>0 )

+                {

+                    if ( is_fst*is_fwd>0 )

+                        stats->isize_inward[isize]++;

+                    else

+                        stats->isize_outward[isize]++;

+                }

+                else if ( is_fst*pos_fst<0 )

+                {

+                    if ( is_fst*is_fwd>0 )

+                        stats->isize_outward[isize]++;

+                    else

+                        stats->isize_inward[isize]++;

+                }

+            }

+        }

+

+        // Number of mismatches 

+        uint8_t *nm = bam_aux_get(bam_line,"NM");

+        if (nm) 

+            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;

+        if ( stats->regions )

+        {

+            // Count only on-target bases

+            int iref = bam_line->core.pos + 1;

+            for (i=0; i<bam_line->core.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

+        {

+            // Count the whole read

+            for (i=0; i<bam_line->core.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;

+

+        if ( stats->tid==bam_line->core.tid && bam_line->core.pos<stats->pos )

+            stats->is_sorted = 0;

+        stats->pos = bam_line->core.pos;

+

+        if ( stats->is_sorted )

+        {

+            if ( stats->tid==-1 || stats->tid!=bam_line->core.tid )

+                round_buffer_flush(stats,-1);

+

+            // Mismatches per cycle and GC-depth graph. For simplicity, reads overlapping GCD bins

+            //  are not splitted which results in up to seq_len-1 overlaps. The default bin size is

+            //  20kbp, so the effect is negligible.

+            if ( stats->fai )

+            {

+                int inc_ref = 0, inc_gcd = 0;

+                // First pass or new chromosome

+                if ( stats->rseq_pos==-1 || stats->tid != bam_line->core.tid ) { inc_ref=1; inc_gcd=1; }

+                // Read goes beyond the end of the rseq buffer

+                else if ( stats->rseq_pos+stats->nrseq_buf < bam_line->core.pos+readlen ) { inc_ref=1; inc_gcd=1; }

+                // Read overlaps the next gcd bin

+                else if ( stats->gcd_pos+stats->gcd_bin_size < bam_line->core.pos+readlen ) 

+                { 

+                    inc_gcd = 1;

+                    if ( stats->rseq_pos+stats->nrseq_buf < bam_line->core.pos+stats->gcd_bin_size ) inc_ref = 1;

+                }

+                if ( inc_gcd )

+                {

+                    stats->igcd++;

+                    if ( stats->igcd >= stats->ngcd )

+                        realloc_gcd_buffer(stats, readlen);

+                    if ( inc_ref )

+                        read_ref_seq(stats,bam_line->core.tid,bam_line->core.pos);

+                    stats->gcd_pos = bam_line->core.pos;

+                    stats->gcd[ stats->igcd ].gc = fai_gc_content(stats, stats->gcd_pos, stats->gcd_bin_size);

+                }

+

+                count_mismatches_per_cycle(stats,bam_line);

+            }

+            // No reference and first pass, new chromosome or sequence going beyond the end of the gcd bin

+            else if ( stats->gcd_pos==-1 || stats->tid != bam_line->core.tid || bam_line->core.pos - stats->gcd_pos > stats->gcd_bin_size )

+            {

+                // First pass or a new chromosome

+                stats->tid     = bam_line->core.tid;

+                stats->gcd_pos = bam_line->core.pos;

+                stats->igcd++;

+                if ( stats->igcd >= stats->ngcd )

+                    realloc_gcd_buffer(stats, readlen);

+            }

+            stats->gcd[ stats->igcd ].depth++;

+            // When no reference sequence is given, approximate the GC from the read (much shorter window, but otherwise OK)

+            if ( !stats->fai )

+                stats->gcd[ stats->igcd ].gc += (float) gc_count / seq_len;

+

+            // Coverage distribution graph

+            round_buffer_flush(stats,bam_line->core.pos);

+            round_buffer_insert_read(&(stats->cov_rbuf),bam_line->core.pos,bam_line->core.pos+seq_len-1);

+        }

+    }

+

+    stats->total_len += seq_len;

+    if ( IS_DUP(bam_line) )

+    {

+        stats->total_len_dup += seq_len;

+        stats->nreads_dup++;

+    }

+}

+

+// Sort by GC and depth

+#define GCD_t(x) ((gc_depth_t *)x)

+static int gcd_cmp(const void *a, const void *b)

+{