Merge branch 'master' of github.com:smithlabcode/dnmtools into selectsites-input-file-type

Author: andrewdavidsmith

src/analysis/bsrate.cpp

@@ -17,6 +17,8 @@
  * General Public License for more details.
  */
 
+#include <bamxx.hpp>
+
 #include <algorithm>
 #include <fstream>
 #include <iostream>
@@ -32,14 +34,13 @@
 #include "dnmt_error.hpp"
 #include "smithlab_utils.hpp"
 
-#include <bamxx.hpp>
-
 using std::accumulate;
 using std::cerr;
 using std::cout;
 using std::endl;
 using std::max;
 using std::numeric_limits;
+using std::pair;
 using std::runtime_error;
 using std::string;
 using std::unordered_map;
@@ -48,10 +49,12 @@ using std::vector;
 
 using bamxx::bam_rec;
 
-static void
+static pair<uint32_t, uint32_t>
 count_states_pos(const bool INCLUDE_CPGS, const string &chrom,
                  const bam_rec &aln, vector<size_t> &unconv,
                  vector<size_t> &conv, vector<size_t> &err, size_t &hanging) {
+  uint32_t n_conv = 0, n_uconv = 0;
+
   /* iterate through reference, query/read and fragment */
   const auto seq = bam_get_seq(aln);
   const auto beg_cig = bam_get_cigar(aln);
@@ -73,10 +76,14 @@ count_states_pos(const bool INCLUDE_CPGS, const string &chrom,
             (rpos >= chrom_lim || !is_guanine(chrom[rpos + 1]) ||
              INCLUDE_CPGS)) {
           const auto qc = seq_nt16_str[bam_seqi(seq, qpos)];
-          if (qc == 'C')
+          if (qc == 'C') {
             ++unconv[fpos];
-          else if (qc == 'T')
+            ++n_uconv;
+          }
+          else if (qc == 'T') {
             ++conv[fpos];
+            ++n_conv;
+          }
           else if (qc != 'N')
             ++err[fpos];
         }
@@ -90,12 +97,15 @@ count_states_pos(const bool INCLUDE_CPGS, const string &chrom,
   }
 
   assert(qpos == get_l_qseq(aln));
+  return {n_conv, n_conv + n_uconv};
 }
 
-static void
+static pair<uint32_t, uint32_t>
 count_states_neg(const bool INCLUDE_CPGS, const string &chrom,
                  const bam_rec &aln, vector<size_t> &unconv,
                  vector<size_t> &conv, vector<size_t> &err, size_t &hanging) {
+  uint32_t n_conv = 0, n_uconv = 0;
+
   /* iterate backward over query/read positions but forward over
      reference and fragment positions */
   const auto seq = bam_get_seq(aln);
@@ -117,10 +127,14 @@ count_states_neg(const bool INCLUDE_CPGS, const string &chrom,
         if (is_guanine(chrom[rpos]) &&
             (rpos == 0 || !is_cytosine(chrom[rpos - 1]) || INCLUDE_CPGS)) {
           const auto qc = seq_nt16_str[bam_seqi(seq, qpos - 1)];
-          if (qc == 'C')
+          if (qc == 'C') {
             ++unconv[fpos];
-          else if (qc == 'T')
+            ++n_uconv;
+          }
+          else if (qc == 'T') {
             ++conv[fpos];
+            ++n_conv;
+          }
           else if (qc != 'N')
             ++err[fpos];
         }
@@ -133,6 +147,7 @@ count_states_neg(const bool INCLUDE_CPGS, const string &chrom,
     }
   }
   assert(qpos == 0);
+  return {n_conv, n_conv + n_uconv};
 }
 
 static void
@@ -223,16 +238,65 @@ write_output(const string &outfile, const vector<size_t> &ucvt_count_p,
   }
 }
 
+template<typename T> static inline void
+update_per_read_stats(const pair<T, T> &x, vector<vector<T>> &tab) {
+  if (x.second < tab.size()) ++tab[x.second][x.first];
+}
+
+static inline vector<double>
+format_histogram(const vector<vector<uint32_t>> &tab,
+                 const size_t n_hist_bins) {
+  constexpr auto epsilon = 1e-6;
+  vector<double> hist(n_hist_bins, 0.0);
+  for (size_t i = 1; i < tab.size(); ++i) {
+    const double denom = i + epsilon;
+    for (size_t j = 0; j <= i; ++j) {
+      const double frac = j / denom;
+      const auto bin_id = std::floor(frac * n_hist_bins);
+      assert(bin_id < hist.size());
+      hist[bin_id] += tab[i][j];
+    }
+  }
+  return hist;
+}
+
+static void
+write_hanging_read_message(std::ostream &out, const size_t n_hanging) {
+  out << "Warning: hanging reads detected at chromosome ends "
+      << "(N=" << n_hanging << ")." << endl
+      << "High numbers of hanging reads suggest inconsistent " << endl
+      << "reference genomes between stages of analysis." << endl
+      << "This is likely to result in analysis errors." << endl;
+}
+
+template<typename T> static void
+write_per_read_histogram(const vector<vector<T>> &tab, const size_t n_hist_bins,
+                         std::ostream &out) {
+  const auto hist = format_histogram(tab, n_hist_bins);
+  out << std::fixed;
+  for (auto i = 0.0; i < hist.size(); ++i)
+    out << std::setprecision(3) << i / hist.size() << '\t'
+        << std::setprecision(3) << (i + 1) / hist.size() << '\t'
+        << std::setprecision(0) << hist[i] << endl;
+}
+
 int
 main_bsrate(int argc, const char **argv) {
   try {
-    // ASSUMED MAXIMUM LENGTH OF A FRAGMENT
-    static const size_t output_size = 10000;
+    // assumed maximum length of a fragment
+    constexpr const size_t output_size = 10000;
+
+    // Assumed maximum cytosines per fragment. Currently the per-read
+    // information is collected as counts in a 2D array, so that later
+    // features may be able to fit distributions based these counts.
+    constexpr const size_t max_cytosine_per_frag = 1000;
 
     bool VERBOSE = false;
     bool INCLUDE_CPGS = false;
     bool reads_are_a_rich = false;
+    bool report_per_read = false;
     size_t n_threads = 1;
+    size_t n_hist_bins = 20;
 
     string chroms_file;
     string outfile;
@@ -254,6 +318,10 @@ main_bsrate(int argc, const char **argv) {
                       seq_to_use);
     opt_parse.add_opt("a-rich", 'A', "reads are A-rich", false,
                       reads_are_a_rich);
+    opt_parse.add_opt("per-read", 'p', "report per-read conversion to terminal",
+                      false, report_per_read);
+    opt_parse.add_opt("bins", 'b', "number of bins for per-read", false,
+                      n_hist_bins);
     opt_parse.add_opt("threads", 't', "number of threads", false, n_threads);
     opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE);
     vector<string> leftover_args;
@@ -295,8 +363,7 @@ main_bsrate(int argc, const char **argv) {
     bamxx::bam_header hdr(hts);
     if (!hdr) throw dnmt_error("failed to read header");
 
-    if (n_threads > 1)
-      tp.set_io(hts);
+    if (n_threads > 1) tp.set_io(hts);
 
     // map the bam header index for each "target" to a sequence in the
     // reference genome
@@ -307,6 +374,9 @@ main_bsrate(int argc, const char **argv) {
       chrom_lookup.insert({sam_hdr_name2tid(hdr.h, names[i].data()), i});
     }
 
+    vector<vector<uint32_t>> per_read_counts(
+      max_cytosine_per_frag, vector<uint32_t>(max_cytosine_per_frag, 0));
+
     vector<size_t> unconv_count_pos(output_size, 0ul);
     vector<size_t> conv_count_pos(output_size, 0ul);
     vector<size_t> unconv_count_neg(output_size, 0ul);
@@ -324,7 +394,6 @@ main_bsrate(int argc, const char **argv) {
     unordered_set<int32_t> chroms_seen;
 
     while (hts.read(hdr, aln)) {
-
       if (reads_are_a_rich) flip_conversion(aln);
 
       // get the correct chrom if it has changed
@@ -351,24 +420,25 @@ main_bsrate(int argc, const char **argv) {
       }
 
       if (use_this_chrom) {
-        // do the work for this mapped read
-        if (bam_is_rev(aln))
-          count_states_neg(INCLUDE_CPGS, chroms[chrom_idx], aln,
-                           unconv_count_neg, conv_count_neg, err_neg, hanging);
-        else
-          count_states_pos(INCLUDE_CPGS, chroms[chrom_idx], aln,
-                           unconv_count_pos, conv_count_pos, err_pos, hanging);
+        // do the work for the current mapped read
+        const auto conv_result =
+          bam_is_rev(aln) ? count_states_neg(INCLUDE_CPGS, chroms[chrom_idx],
+                                             aln, unconv_count_neg,
+                                             conv_count_neg, err_neg, hanging)
+                          : count_states_pos(INCLUDE_CPGS, chroms[chrom_idx],
+                                             aln, unconv_count_pos,
+                                             conv_count_pos, err_pos, hanging);
+        if (report_per_read)
+          update_per_read_stats(conv_result, per_read_counts);
       }
     }
     write_output(outfile, unconv_count_pos, conv_count_pos, unconv_count_neg,
                  conv_count_neg, err_pos, err_neg);
 
-    if (hanging > 0)  // some overhanging reads
-      cerr << "Warning: hanging reads detected at chrom ends "
-           << "(N=" << hanging << ")" << endl
-           << "High numbers of hanging reads suggest mismatch "
-           << "between assembly provided here and that used for mapping"
-           << endl;
+    if (hanging > 0) write_hanging_read_message(cerr, hanging);
+
+    if (report_per_read)
+      write_per_read_histogram(per_read_counts, n_hist_bins, cout);
   }
   catch (const runtime_error &e) {
     cerr << e.what() << endl;

src/analysis/hmr-rep.cpp

@@ -415,6 +415,10 @@ main_hmr_rep(int argc, const char **argv) {
     const vector<string> cpgs_files(leftover_args);
     /****************** END COMMAND LINE OPTIONS *****************/
 
+    for (auto &filename : cpgs_files)
+      if (!is_msite_file(filename))
+        throw runtime_error("malformed counts file: " + filename);
+
     vector<string> params_in_file;
     if (!params_in_files.empty()) {
       params_in_file = split_comma(params_in_files);

src/analysis/hmr.cpp

@@ -442,6 +442,9 @@ main_hmr(int argc, const char **argv) {
     const string cpgs_file = leftover_args.front();
     /****************** END COMMAND LINE OPTIONS *****************/
 
+    if (!is_msite_file(cpgs_file))
+      throw runtime_error("malformed counts file: " + cpgs_file);
+
     // separate the regions by chrom and by desert
     vector<MSite> cpgs;
     vector<pair<double, double> > meth;

src/analysis/hypermr.cpp

@@ -278,6 +278,9 @@ main_hypermr(int argc, const char **argv) {
     const string cpgs_file = leftover_args.front();
     /****************** END COMMAND LINE OPTIONS *****************/
 
+    if (!is_msite_file(cpgs_file))
+      throw runtime_error("malformed counts file: " + cpgs_file);
+
     if (VERBOSE) cerr << "[loading_data]" << endl;
     vector<MSite> cpgs;
     vector<pair<double, double>> meth;

src/analysis/levels.cpp

@@ -120,6 +120,9 @@ main_levels(int argc, const char **argv) {
     const string meth_file = leftover_args.front();
     /****************** END COMMAND LINE OPTIONS *****************/
 
+    if (!is_msite_file(meth_file))
+      throw runtime_error("malformed counts file: " + meth_file);
+
     const counts_file_format guessed_format =
       guess_counts_file_format(meth_file);
     if (guessed_format != counts_file_format::ordinary) {