@@ -1,5 +1,5 @@

 Package: seqbias

-Version: 1.1.1

+Version: 1.1.3

 Date: 25-12-2010

 Title: Estimation of per-position bias in high-throughput sequencing data

 Description: This package implements a model of per-position sequencing bias in

@@ -109,21 +109,27 @@ setClass( "seqbias",

-"count.reads" <- function( reads_fn, I, binary = TRUE )

+"count.reads" <- function( reads_fn, I, sb = NULL, binary = FALSE, sum.counts = FALSE )

 {

     require(GenomicRanges)

     stopifnot( is( I, "GRanges" ) )

+    if (!is.null(sb) & class(sb) != "seqbias") {

+        stop( "The 'sb' parameter of count.reads must be a seqbias class." )

+    }

+

     bam_ptr <- .Call( "seqbias_open_bam", path.expand(reads_fn),

                       PACKAGE = "seqbias" )

     counts <- tapply( I,

                       INDEX = 1:length(I),

                       FUN   = function(x) .Call( "seqbias_count_reads",

+                                                 sb,

                                                  bam_ptr,

                                                  as.character(seqnames(x)),

                                                  start(x), end(x),

                                                  as.character(strand(x)),

+                                                 sum.counts,

                       PACKAGE = "seqbias" ) )

     if( binary ) lapply( counts, FUN = function(c) as.integer( c > 0 ) )

@@ -126,7 +126,7 @@ seqs[neg_idx] <- reverseComplement( seqs[neg_idx] )

 Finally, we count the number of reads mapping to each position in our sampled

 intervals.

 <<>>=

-counts <- count.reads( reads_fn, I )

+counts <- count.reads( reads_fn, I, binary = T )

 @

 Unless the \texttt{binary} argument is FALSE, this function returns a 0-1 vector,

@@ -192,7 +192,7 @@ more reads are available a more accurate model can be trained at the expense of

 the training procedure taking up to several minutes.

 <<results=hide>>=

-sb <- seqbias.fit( ref_fn, reads_fn, L = 5, R = 20 )

+sb <- seqbias.fit( ref_fn, reads_fn, L = 5, R = 15 )

 @

 The \texttt{L} and \texttt{R} arguments control the maximum number of positions

@@ -3,21 +3,25 @@

 \title{Counting reads across intervals}

 \description{Counts the number of reads starting at each position across given

 genomic intervals}

-\usage{count.reads(reads_fn, I, binary=TRUE)}

+\usage{count.reads(reads_fn, I, sb=NULL, binary=FALSE, sum.counts=TRUE)}

 \arguments{

   \item{reads_fn}{filename of aligned reads in BAM format}

   \item{I}{a GRanges object giving valid genomic intervals}

+  \item{sb}{a seqbias object}

   \item{binary}{if \code{TRUE}, return a 0-1 vector, otherwise return a vector

   counting the number of reads mapped to each position}

+  \item{sum.counts}{if \code{TRUE} return the total read count for each interval}

 }

 \details{

   Given an indexed BAM file, this function counts the number of reads starting

-  at each position of each provided interval. These counts can then be

-  normalized by dividing by the predicted bias obtained from 'bias.predict'.

-

-  By default, a 0-1 vector is returned, where positions at which no reads are

-  mapped are 0, and those with one or more are 1. If \code{binary} is

-  \code{FALSE}, the number of reads mapping to each position is returned.

+  at each position of each provided interval. If a seqbias object is provided

+  through the \code{sb} attribute, counts are corrected for sequence bias.

+  The total read count for each interval is returned if \code{sum.counts} is

+  \code{TRUE}.

+  

+  If \code{binary} is \code{TRUE} a 0-1 vector is returned instead, where

+  positions at which no reads are mapped are 0, and those with one or more are

+  1.

 }

 \value{

   A list of numeric vectors is returned, one for each interval provided. Each

@@ -36,7 +36,6 @@ kmer_matrix::kmer_matrix(const YAML::Node& node)

     A = new double[size1 * size2];

     const YAML::Node& node_A = node["A"];

-

     size_t i;

     for (i = 0; i < size1 * size2; ++i) {

         node_A[i] >> A[i];

@@ -186,7 +186,8 @@ class motif_trainer

             compute_reachability();

             logger::print(col_base);

-            snprintf(msg, sizeof(msg), "round %4zu (ic = %0.4e)", round_num, ic);

+            snprintf(msg, sizeof(msg), "round %4lu (ic = %0.4e)",

+                                       (unsigned long)round_num, ic);

             logger::print(msg);

             logger::print(col_base);

             col = 0;

@@ -442,7 +443,6 @@ class motif_trainer

             else i_start = i_end = j;

             for (i = i_start; i <= i_end; ++i) {

-

                 /* skip existing edges */

                 if (M.has_edge(i, j)) continue;

@@ -800,7 +800,7 @@ motif::motif(const YAML::Node& node)

     /* m */

     unsigned int m_;

     node["m"] >> m_;

-    m = (size_t) m;

+    m = (size_t) m_;

     /* parents */

@@ -809,8 +809,10 @@ motif::motif(const YAML::Node& node)

     const YAML::Node& parents_node = node["parents"];

     size_t i;

-    for (i = 0; i < m; ++i) {

-        parents_node[i] >> parents[i];

+    int b;

+    for (i = 0; i < m * m; ++i) {

+        parents_node[i] >> b;

+        parents[i] = (bool) b;

     }

@@ -840,7 +842,7 @@ void motif::to_yaml(YAML::Emitter& out) const

     out << YAML::Flow << YAML::BeginSeq;

     size_t i;

-    for (i = 0; i < m * m; ++i) out << parents[i];

+    for (i = 0; i < m * m; ++i) out << (parents[i] ? 1 : 0);

     out << YAML::EndSeq;

@@ -860,6 +862,45 @@ void motif::to_yaml(YAML::Emitter& out) const

 }

+string motif::model_graph(int offset) const

+{

+    string graph_str;

+    char strbuf[512];

+

+    graph_str += "digraph {\n";

+    graph_str += "splines=\"true\";\n";

+    graph_str += "node [shape=\"box\"];\n";

+

+    /* print nodes */

+    size_t i, j;

+    for (j = 0; j < m; j++) {

+        snprintf(strbuf, sizeof(strbuf), 

+                 "n%d [label=\"%d\",pos=\"%d,0\",style=\"%s\"];\n",

+                 (int) j, (int) j - offset, (int) j * 100,

+                 parents[j * m + j] ? "solid" : "dotted");

+

+        graph_str += strbuf;

+    }

+

+    /* print edges */

+    for (j = 0; j < m; ++j) {

+        if (!parents[j * m + j]) continue;

+

+        for (i = 0; i < m; ++i) {

+            if (i == j) continue;

+            if (parents[j * m + i]) {

+                snprintf(strbuf, sizeof(strbuf),

+                         "n%lu -> n%lu;\n", (unsigned long) i, (unsigned long) j);

+                graph_str += strbuf;

+            }

+        }

+    }

+

+    graph_str += "}\n";

+    return graph_str;

+}

+

+

 double motif::eval(const twobitseq& seq, size_t offset)

 {

     double ll0 = 0.0;

@@ -45,6 +45,10 @@ class motif

         /** emit yaml serialization */

         void to_yaml(YAML::Emitter& out) const;

+        /** print the model graph in GraphViz (dot) format,

+         *  adjusting position labels by the given offset. */

+        std::string model_graph(int offset) const;

+

     private:

         motif();

@@ -153,8 +153,13 @@ SEXP seqbias_predict( SEXP seqbias,

     pos_t c_start, c_end;

     strand_t c_strand;

+

     coerce_genomic_coords( seqname, start, end, strand,

                            &c_seqname, &c_start, &c_end, &c_strand );

+

+    /* we expect 1-based coordinates, but work with 0-based */

+    c_start -= 1;

+    c_end   -= 1;

     if( c_strand != 0 && c_strand != 1 ) {

         warning( "strand should be '+' or '-'" );

@@ -217,11 +222,13 @@ SEXP seqbias_open_bam( SEXP reads_fn )

 }

-SEXP seqbias_count_reads( SEXP bam_ptr,

+SEXP seqbias_count_reads( SEXP seqbias,

+                          SEXP bam_ptr,

                           SEXP seqname,

                           SEXP start,

                           SEXP end,

-                          SEXP strand )

+                          SEXP strand,

+                          SEXP sum_counts )

 {

     if( TYPEOF(bam_ptr) != EXTPTRSXP ) error( "argument is not a indexed bam pointer" );

     indexed_bam_f* c_bam_ptr = (indexed_bam_f*)EXTPTR_PTR( bam_ptr );

@@ -232,13 +239,52 @@ SEXP seqbias_count_reads( SEXP bam_ptr,

     coerce_genomic_coords( seqname, start, end, strand,

                            &c_seqname, &c_start, &c_end, &c_strand );

+

+

+    /* we expect 1-based coordinates, but use 0-based */

+    c_start -= 1;

+    c_end   -= 1;

+    double* bs[2] = {NULL, NULL};

+

+    if( !isNull(seqbias) ) {

+

+        SEXP idx;

+        PROTECT( idx = allocVector( STRSXP, 1 ) );

+        SET_STRING_ELT(idx, 0, mkChar("ptr"));

+

+        sequencing_bias* c_seqbias = NULL;

+        if( TYPEOF(GET_SLOT(seqbias, idx)) != EXTPTRSXP ||

+            !(c_seqbias = (sequencing_bias*)EXTPTR_PTR(GET_SLOT(seqbias, idx))) ) {

+            error( "first argument is not a proper seqbias class." );

+        }

+        

+        if (c_strand == strand_na || c_strand == strand_pos ) {

+            bs[0] = c_seqbias->get_bias(c_seqname, c_start, c_end, strand_pos);

+        }

+

+        if (c_strand == strand_na || c_strand == strand_neg ) {

+            bs[1] = c_seqbias->get_bias(c_seqname, c_start, c_end, strand_neg);

+            std::reverse(bs[1], bs[1] + (c_end - c_start + 1));

+        }

+

+        UNPROTECT(1);

+    }

+

+    bool c_sum_counts = asLogical(sum_counts) == TRUE;

+

     /* init vector */

     SEXP v;

-    PROTECT( v = allocVector( REALSXP, c_end - c_start + 1 ) );

-    pos_t i;

-    for( i = 0; i < c_end - c_start + 1; i++ ) REAL(v)[i] = 0;

+    if (c_sum_counts) {

+        PROTECT( v = allocVector( REALSXP, 1 ) );

+        REAL(v)[0] = 0;

+    }

+    else {

+        PROTECT( v = allocVector( REALSXP, c_end - c_start + 1 ) );

+        pos_t i;

+        for( i = 0; i < c_end - c_start + 1; i++ ) REAL(v)[i] = 0;

+    }

     const size_t region_len = 1024;

@@ -261,26 +307,41 @@ SEXP seqbias_count_reads( SEXP bam_ptr,

     bam_iter_t it = bam_iter_query( c_bam_ptr->idx, bam_ref_id, bam_start, bam_end );

     bam1_t* b = bam_init1();

     pos_t x;

+    strand_t s;

+

     while( bam_iter_read( c_bam_ptr->f->x.bam, it, b ) > 0 ) {

-        if( bam1_strand(b) !=  c_strand ) continue;

+        s = (strand_t) bam1_strand(b);

+        if( c_strand != strand_na && s !=  c_strand ) continue;

+

         x = bam1_strand(b) == 1 ? bam_calend( &b->core,  bam1_cigar(b) ) - 1 : b->core.pos;

-        if( c_start <= x && x <= c_end ) REAL(v)[x - c_start]++;

+        if (x < c_start || x > c_end) continue;

+

+        if (c_sum_counts) {

+            REAL(v)[0] += bs[s] == NULL ? 1.0 : (1.0 / bs[s][x - c_start]);

+        }

+        else {

+            REAL(v)[x - c_start] += bs[s] == NULL ? 1.0 : (1.0 / bs[s][x - c_start]);

+        }

     }

-    if( c_strand == strand_neg ) {

+    if( c_strand == strand_neg && !c_sum_counts ) {

         std::reverse(REAL(v), REAL(v) + (c_end - c_start + 1));

     }

     bam_iter_destroy( it );

     bam_destroy1(b);

+    delete [] bs[0];

+    delete [] bs[1];

+

     UNPROTECT(1);

     return v;

 }

+

 void dealloc_kmer_matrix( SEXP M )

 {

     if( TYPEOF(M) != EXTPTRSXP ) error( "argument is not a kmer_matrix pointer" );

@@ -433,7 +494,7 @@ void R_init_seqbias( DllInfo* info )

         { "seqbias_load",        (DL_FUNC) &seqbias_load, 2 },

         { "seqbias_save",        (DL_FUNC) &seqbias_save, 2 },

         { "seqbias_open_bam",       (DL_FUNC) &seqbias_open_bam, 1 },

-        { "seqbias_count_reads",    (DL_FUNC) &seqbias_count_reads, 5 },

+        { "seqbias_count_reads",    (DL_FUNC) &seqbias_count_reads, 7 },

         { "seqbias_alloc_kmer_matrix",          (DL_FUNC) &seqbias_alloc_kmer_matrix, 2 },

         { "seqbias_tally_kmers",                (DL_FUNC) &seqbias_tally_kmers, 4 },

         { "sebqias_dataframe_from_kmer_matrix", (DL_FUNC) &seqbias_dataframe_from_kmer_matrix, 2 },

@@ -11,6 +11,7 @@

 #include "miscmath.h"

 #include "samtools_extra.h"

 #include "samtools/faidx.h"

+#include <climits>

 #include <cmath>

 #include <cctype>

 #include <ctime>

@@ -55,17 +56,6 @@ static double rand_gauss(double sigma)

 }

-static double rand_trunc_gauss(double sigma, double a, double b)

-{

-    double x;

-    do {

-        x = rand_gauss(sigma);

-    } while (x < a || x > b);

-

-    return x;

-}

-

-

 double gauss_pdf (const double x, const double sigma)

 {

   double u = x / fabs (sigma);

@@ -86,25 +76,30 @@ static int read_pos_count_compare(const void* p1, const void* p2)

 }

-/*

-static int read_pos_tid_count_compare(const void* p1, const void* p2)

-{

-    int c = (int)(((read_pos*)p1)->tid - ((read_pos*)p2)->tid);

-

-    if (c == 0) {

-        return (int)((read_pos*)p2)->count - (int)((read_pos*)p1)->count;

-    }

-    else return c;

-}

-*/

-

-

 sequencing_bias::sequencing_bias()

     : ref_f(NULL)

     , M(NULL)

 {}

+

+sequencing_bias::sequencing_bias(const char* model_fn)

+    : ref_f(NULL)

+    , M(NULL)

+{

+    std::ifstream fin;

+    fin.open(model_fn);

+

+    YAML::Parser parser(fin);

+    YAML::Node doc;

+    parser.GetNextDocument(doc);

+

+    doc["L"] >> L;

+    doc["R"] >> R;

+

+    M = new motif(doc["motif"]);

+}

+

 sequencing_bias::sequencing_bias(const char* ref_fn,

                                  const char* model_fn)

     : ref_f(NULL)

@@ -122,12 +117,15 @@ sequencing_bias::sequencing_bias(const char* ref_fn,

     M = new motif(doc["motif"]);

-    this->ref_fn = ref_fn;

-    ref_f = fai_load(ref_fn);

-    if (ref_f == NULL) {

-        logger::abort("Can't open indexed FASTA file %s\n", ref_fn);

+    if (ref_fn != NULL) {

+        this->ref_fn = ref_fn;

+        ref_f = fai_load(ref_fn);

+        if (ref_f == NULL) {

+            logger::abort("Can't open indexed FASTA file %s\n", ref_fn);

+        }

     }

+    else ref_f = NULL;

 }

@@ -233,7 +231,6 @@ void sequencing_bias::build(const char* ref_fn,

     while (samread(reads_f, read) > 0) {

         if (read->core.n_cigar != 1) continue;

-        if (read->core.flag & BAM_FREAD2) continue;

         k++;

         if (k % 1000000 == 0) {

             logger::info("hashed %zu reads.", k);

@@ -265,7 +262,7 @@ void sequencing_bias::build(const char* ref_fn,

     this->L = L;

     this->R = R;

-    unsigned int i;

+    unsigned int i = 0;

     read_pos* S_tmp;

     read_pos* S;

@@ -273,34 +270,12 @@ void sequencing_bias::build(const char* ref_fn,

     const size_t max_dump = 10000000;

     pos_table_dump(T, &S_tmp, &N, max_dump);

-    /* sort by count */

-    qsort(S_tmp, N, sizeof(read_pos), read_pos_count_compare);

-

-    /* consider only reads with at least one duplicate */

-    for (i = 0; i < N; ++i) {

-        if (S_tmp[i].count <= 1) break;

-    }

-

-    /* (unless there are very few of these reads */

-    if (i > 10000) {

-        max_reads = std::min<size_t>(max_reads, i);

-        logger::info("%zu reads with duplicates.", i);

-    }

-    else {

-        i = N;

-    }

-

-

-    /* ignore the top 1%, as they tend to be vastly higher than anything else,

-     * and thus can throw things off when training with small numbers of reads

-     * */

-    S = S_tmp + i/100;

-    max_reads = min<size_t>(max_reads, 99*i/100); 

+    S = S_tmp;

     /* sort by tid (so we can load one chromosome at a time) */

+    random_shuffle(S, S + N);

     qsort(S, std::min<size_t>(max_reads, N), sizeof(read_pos), read_pos_tid_compare);

-

     /* sample foreground and background kmer frequencies */

     ref_f = fai_load(ref_fn);

     if (ref_f == NULL) {

@@ -337,7 +312,7 @@ void sequencing_bias::build(const char* ref_fn,

             seqname = T->seq_names[ S[i].tid ];

             if (seq) free(seq); 

-            seq = faidx_fetch_seq(ref_f, seqname, 0, 1000000000, &seqlen);

+            seq = faidx_fetch_seq(ref_f, seqname, 0, INT_MAX, &seqlen);

             logger::info("read sequence %s.", seqname);

@@ -374,7 +349,7 @@ void sequencing_bias::build(const char* ref_fn,

         /* adjust the current read position randomly, and sample */

         for (bg_sample_num = 0; bg_sample_num < bg_samples;) {

-            bg_pos = S[i].pos + (pos_t)round_away(rand_trunc_gauss(12, -200, 200));

+            bg_pos = S[i].pos + (pos_t)round_away(rand_gauss(500));

             if (S[i].strand == strand_neg) {

                 if (bg_pos < R || bg_pos >= seqlen - L) continue;

@@ -464,17 +439,201 @@ double* sequencing_bias::get_bias(const char* seqname,

         bs[i] = M->eval(seq, i);

     }

-

     free(seqstr);

     return bs;

 }

-string sequencing_bias::print_model_graph()

+string sequencing_bias::model_graph() const

+{

+    return M->model_graph((int) L);

+}

+

+

+kmer_matrix tabulate_bias(double* kl,

+                          pos_t L, pos_t R, int k,

+                          const char* ref_fn,

+                          const char* reads_fn,

+                          const char* model_fn)

 {

-    // TODO

-    //return M->print_model_graph(L);

-    return string();

+    /* This function procedes very much like the sequencing_bias::build

+     * function, but does not finally train a motif, but rather tabulates

+     * k-mer frequencies. */

+    size_t max_reads = 250000;

+

+    kmer_matrix dest((size_t) (L + 1 + R), (size_t) k);

+    dest.set_all(0.0);

+

+    faidx_t* ref_f = fai_load(ref_fn);

+    if (ref_f == NULL) {

+        logger::abort("Can't open fasta file '%s'.", ref_fn);

+    }

+

+    samfile_t* reads_f = samopen(reads_fn, "rb", NULL);

+    if (reads_f == NULL) {

+        logger::abort("Can't open bam file '%s'.", reads_fn);

+    }

+

+    bam_index_t* reads_index = bam_index_load(reads_fn);

+    if (reads_index == NULL) {

+        logger::abort("Can't open bam index '%s.bai'.", reads_fn);

+    }

+

+    sequencing_bias* sb = NULL;

+    if (model_fn != NULL) {

+        sb = new sequencing_bias(ref_fn, model_fn);

+    }

+

+

+    bam_init_header_hash(reads_f->header);

+    bam1_t* read = bam_init1();

+

+    pos_table T;

+    pos_table_create(&T, reads_f->header->n_targets);

+    T.seq_names = reads_f->header->target_name;

+

+    size_t hashed_count = 0;

+    while (samread(reads_f, read) > 0) {

+        if (read->core.n_cigar != 1) continue;

+        if (++hashed_count % 1000000 == 0) {

+            logger::info("hashed %zu reads.", hashed_count);

+        }

+        pos_table_inc(&T, read);

+    }

+    logger::info("hashed %zu reads.", hashed_count);

+

+    read_pos* S_tmp;

+    read_pos* S;

+    size_t N;

+    const size_t max_dump = 10000000;

+    pos_table_dump(&T, &S_tmp, &N, max_dump);

+

+    /* sort by count */

+    qsort(S_tmp, N, sizeof(read_pos), read_pos_count_compare);

+

+    /* consider only reads with at least one duplicate */

+    size_t i;

+    for (i = 0; i < N; ++i) {

+        if (S_tmp[i].count <= 1) break;

+    }

+

+    /* (unless there are very few of these reads */

+    if (i > 10000) {

+        max_reads = std::min<size_t>(max_reads, i);

+        logger::info("%zu reads with duplicates.", i);

+    }

+    else {

+        i = N;

+    }

+

+    /* ignore the top 1%, as they tend to be vastly higher than anything else,

+     * and thus can throw things off when training with small numbers of reads

+     * */

+    S = S_tmp + i/100;

+    max_reads = min<size_t>(max_reads, 99*i/100); 

+

+    /* sort by tid (so we can load one chromosome at a time) */

+    qsort(S, std::min<size_t>(max_reads, N), sizeof(read_pos), read_pos_tid_compare);

+

+    twobitseq tbs;

+

+    char* local_seq;

+    local_seq = new char[ (k - 1) + L + 1 + R + 1 ];

+    local_seq[(k - 1) + L + 1 + R] = '\0';

+

+    double         w;

+    char*          seq       = NULL;

+    int            seqlen    = 0;

+    int            curr_tid  = -1;

+    char*          seqname   = NULL;

+

+    for (i = 0; i < N && i < max_reads; i++) {

+        if (S[i].tid != curr_tid) {

+            seqname = T.seq_names[S[i].tid];

+            free(seq);

+            seq = faidx_fetch_seq(ref_f, seqname, 0, INT_MAX, &seqlen);

+            logger::info("read sequence %s.", seqname);

+            curr_tid = S[i].tid;

+

+            if (seq == NULL) {

+                logger::warn("warning: reference sequence not found, skipping.");

+            }

+        }

+

+        if (seq == NULL) continue;

+

+        if (S[i].strand == strand_neg) {

+            if (S[i].pos < R || S[i].pos >= seqlen - L - (k-1)) continue;

+            memcpy(local_seq, seq + S[i].pos - R, ((k-1)+L+1+R)*sizeof(char));

+            seqrc(local_seq, L+1+R);

+        }

+        else {

+            if (S[i].pos < L + (k-1) || S[i].pos >= seqlen - R) continue;

+            memcpy(local_seq, seq + (S[i].pos - L - (k-1)), ((k-1)+L+1+R)*sizeof(char));

+        }

+

+        if (sb) {

+            // TODO: get_bias

+        }

+        else w = 1.0;

+

+        tbs = local_seq;

+        kmer K;

+        for (pos_t pos = (k-1); pos < (k-1) + L + 1 + R; ++pos) {

+            K = tbs.get_kmer(k, pos);

+            dest(pos - (k-1), K) += w;

+        }

+    }

+

+    /* compute KL divergence */

+

+    /* estimate a background distribution by averaging across the entire window

+     * */

+    int four_to_k = 1 << (2 * k);

+    double* bg = new double[four_to_k];

+    memset(bg, 0, four_to_k * sizeof(double));

+    for (pos_t pos = 0; pos < L + 1 + R; ++pos) {

+        for (kmer K = 0; K < (kmer) four_to_k; ++K) {

+            bg[K] += dest(pos, K);

+        }

+    }

+

+    kmer_matrix norm_dest(dest);

+    norm_dest.make_distribution();

+

+    double z = 0.0;

+    for (kmer K = 0; K < (kmer) four_to_k; ++K) z += bg[K];

+    for (kmer K = 0; K < (kmer) four_to_k; ++K) bg[K] /= z;

+

+

+    /* Compute the (symetric) kullback-leibler divegnnce */

+    memset(kl, 0, (L + 1 + R) * sizeof(double));

+    for (pos_t pos = 0; pos < L + 1 + R; ++pos) {

+        kl[pos] = 0.0;

+        for (kmer K = 0; K < (kmer) four_to_k; ++K) {

+            if (norm_dest(pos, K) > 0.0) {

+                kl[pos] += norm_dest(pos, K) * (log2(norm_dest(pos, K)) - log2(bg[K]));

+            }

+

+            if (bg[K] > 0.0) {

+                kl[pos] += bg[K] * (log2(bg[K]) - log2(norm_dest(pos, K)));

+            }

+        }

+    }

+

+    delete [] bg;

+

+

+    free(seq);

+    free(local_seq);

+    free(S_tmp);

+    bam_destroy1(read);

+    pos_table_destroy(&T);

+    delete sb;

+    bam_index_destroy(reads_index);

+    samclose(reads_f);

+

+    return dest;

 }

@@ -21,7 +21,13 @@

 class sequencing_bias

 {

     public:

-        /** Load a model that has bee previously trained. */

+        /** Load a model that has been previously trained, without a reference

+         * sequence. (It can not be actually used this way.)

+         */

+        sequencing_bias(const char* model_fn);

+

+

+        /** Load a model that has been previously trained. */

         sequencing_bias(const char* ref_fn,

                         const char* model_fn);

@@ -74,7 +80,7 @@ class sequencing_bias

         void to_yaml(YAML::Emitter&) const;

         /** Return a string of the model graph in dot format. */

-        std::string print_model_graph();

+        std::string model_graph() const;

     private:

@@ -107,4 +113,13 @@ class sequencing_bias

 };

+/** Tabulate the bias in a given dataset, optionally correcting for bias using

+ *  the given model. */

+kmer_matrix tabulate_bias(double* kl,

+                          pos_t L, pos_t R, int k,

+                          const char* ref_fn,

+                          const char* reads_fn,

+                          const char* model_fn = NULL);

+

+

 #endif

@@ -8,10 +8,8 @@

 #include "twobitseq.hpp"

 #include <algorithm>

 #include <cstring>

+#include <cstdlib>

-/* Don't expose this function. No one should have to know how we are mapping

- * letters to numbers.

- */

 kmer nuc_to_num(char c)

 {

     switch( c ) {

@@ -26,36 +24,48 @@ kmer nuc_to_num(char c)

     }

 }

-void num_to_nuc(char* dest, int n, int k)

+

+

+void num_to_nuc(char* dest, kmer K, int k)

 {

     int i;

     /* read backwards, then reverse */

-    for( i = 0; i < k; i++ ) {

-        switch( n & 0x3 ) {

+    for (i = 0; i < k; ++i) {

+        switch (K & 0x3) {

             case 0: dest[i] = 'a'; break;

             case 1: dest[i] = 'c'; break;

             case 2: dest[i] = 'g'; break;

             case 3: dest[i] = 't'; break;

         }

-        n >>= 2;

+

+        K >>= 2;

     }

-    dest[i] = '\0';

+    dest[i] = '\0';

     std::reverse(dest, dest + i);

 }

 const size_t twobitseq::max_kmer = 4 * sizeof(kmer);

+twobitseq::twobitseq()

+    : xs(NULL)

+    , n(0)

+{

+

+}

+

 twobitseq::twobitseq(const char* seq)

     : xs(NULL)

     , n(0)

 {

+    if (seq == NULL) return;

     n = strlen(seq);

     if (n == 0) return;

-    xs = new kmer[n / max_kmer + 1];

+    xs = reinterpret_cast<kmer*>(

+            malloc_or_die((n / max_kmer + 1) * sizeof(kmer)));

     memset(xs, 0, (n / max_kmer + 1) * sizeof(kmer));

     size_t i;

@@ -77,7 +87,8 @@ twobitseq::twobitseq(const twobitseq& other)

     n = other.n;

     if (n == 0) return;

-    xs = new kmer[n / max_kmer + 1];

+    xs = reinterpret_cast<kmer*>(

+            malloc_or_die((n / max_kmer + 1) * sizeof(kmer)));

     memcpy(xs, other.xs, (n / max_kmer + 1) * sizeof(kmer));

 }

@@ -85,20 +96,60 @@ twobitseq::twobitseq(const twobitseq& other)

 twobitseq::~twobitseq()

 {

-    delete [] xs;

+    free(xs);

 }

 void twobitseq::operator = (const twobitseq& other)

 {

     n = other.n;

-

-    delete [] xs;

-    xs = new kmer[n / max_kmer + 1];

+    xs = reinterpret_cast<kmer*>(

+            realloc_or_die(xs, (n / max_kmer + 1) * sizeof(kmer)));

     memcpy(xs, other.xs, (n / max_kmer + 1) * sizeof(kmer));

 }

+void twobitseq::operator = (const char* seq)

+{

+    if (seq == NULL) {

+        n = 0;

+        free(xs);

+        xs = NULL;

+        return;

+    }

+

+    n = strlen(seq);

+    xs = reinterpret_cast<kmer*>(

+            realloc_or_die(xs, (n / max_kmer + 1) * sizeof(kmer)));

+    memset(xs, 0, (n / max_kmer + 1) * sizeof(kmer));

+

+    size_t i;

+    size_t block, offset;

+    for (i = 0; i < n; ++i) {

+        block  = i / max_kmer;

+        offset = i % max_kmer;

+

+        xs[block] = xs[block] | (nuc_to_num(seq[i]) << (2*offset));

+    }

+}

+

+

+kmer twobitseq::get_kmer(int k, pos_t pos)

+{

+    size_t block, offset;

+    size_t i;

+    kmer K = 0;

+

+    for (i = 0; i < (size_t) k; ++i) {

+        block  = (i + (size_t) (pos - (k - 1))) / max_kmer;

+        offset = (i + (size_t) (pos - (k - 1))) % max_kmer;