Merge pull request #133 from smithlabcode/pmd-summary

Pmd summary

Author: masarunakajima

docs/README.md

@@ -19,15 +19,15 @@ pip install -U mkdocs
 
 Build the HTML documentation by running
 ```
-mkdocs build
+mkdocs build -f docs/mkdocs.yml
 ```
 which will create a `site` directory where markdown files are
 converted to HTML
 
 Create a local host for the HTML documentation by running
 
 ```
-mkdocs serve
+mkdocs serve -f docs/mkdocs.yml
 ```
 
 This will create the documentation, usually at http://localhost:8000 .

docs/content/hmr.md

@@ -170,3 +170,9 @@ File in which to write parameters learned during the current run.
 A random number seed. Randomization is used in a shuffling step prior
 to filering candidate HMRs. This parameter is typically only used for
 testing (default: 408).
+
+```txt
+ -S, -summary
+```
+Write the analysis summary to this file. The summary is not
+reported unless a file is specified here. This option is correct as of v1.4.0.

docs/content/pmd.md

@@ -186,3 +186,10 @@ Write the trained parameters to this file.
  -s, -seed
 ```
 Specify a random seed value.
+
+```txt
+ -S, -summary
+```
+Write the analysis summary to this file. The summary is not
+reported unless a file is specified here. This option is correct as of v1.4.0.
+

docs/content/uniq.md

@@ -49,10 +49,10 @@ most helpful when the input and output are both BAM, where the threads
 can really speed things up.
 
 ```txt
- -S, -stats
+ -S, -summary
 ```
 Save statistics on duplication rates to this file. The statistics are not
-reported unless a file is specified here.
+reported unless a file is specified here. This option is correct as of v1.4.0.
 
 ```txt
  -hist

src/analysis/hmr.cpp

@@ -50,6 +50,34 @@ using std::unordered_set;
 
 using bamxx::bgzf_file;
 
+struct hmr_summary {
+  hmr_summary(const vector<GenomicRegion> &hmrs) {
+    hmr_count = hmrs.size();
+    hmr_total_size = accumulate(cbegin(hmrs), cend(hmrs), 0,
+                                [](const uint64_t t, const GenomicRegion &p) {
+                                  return t + p.get_width(); });
+    hmr_mean_size = 
+      static_cast<double>(hmr_total_size)/std::max(1ul, hmr_count);
+  }
+  // hmr_count is the number of identified HMRs.
+  uint64_t hmr_count{};
+  // total_hmr_size is the sum of the sizes of the identified HMRs
+  uint64_t hmr_total_size{};
+  // mean_hmr_size is the mean size of the identified HMRs
+  double hmr_mean_size{};
+
+  string tostring() {
+    std::ostringstream oss;
+    oss << "hmr_count: " << hmr_count << endl
+        << "hmr_total_size: " << hmr_total_size << endl
+        << "hmr_mean_size: "
+        << std::fixed << std::setprecision(2) << hmr_mean_size;
+
+    return oss.str();
+  }
+};
+
+
 static GenomicRegion
 as_gen_rgn(const MSite &s) {
   return GenomicRegion(s.chrom, s.pos, s.pos + 1);
@@ -386,6 +414,8 @@ main_hmr(int argc, const char **argv) {
     // corrections for small values
     const double tolerance = 1e-10;
 
+    string summary_file;
+
     string params_in_file;
     string params_out_file;
 
@@ -414,11 +444,16 @@ main_hmr(int argc, const char **argv) {
     opt_parse.add_opt("post-meth", '\0', "output file for single-CpG posteiror "
                       "methylation probability (default: none)",
                       false, meth_post_outfile);
+    opt_parse.add_opt("post-meth", '\0', "output file for single-CpG posteiror "
+                      "methylation probability (default: none)",
+                      false, meth_post_outfile);
     opt_parse.add_opt("params-in", 'P', "HMM parameter file "
                       "(override training)", false, params_in_file);
     opt_parse.add_opt("params-out", 'p', "write HMM parameters to this "
                       "file (default: none)", false, params_out_file);
     opt_parse.add_opt("seed", 's', "specify random seed", false, rng_seed);
+    opt_parse.add_opt("summary", 'S', "write summary output here", false, 
+                      summary_file);
     opt_parse.set_show_defaults();
     vector<string> leftover_args;
     opt_parse.parse(argc, argv, leftover_args);
@@ -557,6 +592,11 @@ main_hmr(int argc, const char **argv) {
         out << cpg << '\n';
       }
     }
+    if (!summary_file.empty()) {
+      std::ofstream summary_out(summary_file);
+      if (!summary_out) throw runtime_error("failed to open: " + summary_file);
+      summary_out << hmr_summary(domains).tostring() << endl;
+    }
   }
   catch (runtime_error &e) {
     cerr << e.what() << endl;

src/analysis/pmd.cpp

@@ -26,6 +26,7 @@
 #include <stdexcept>
 #include <unordered_set>
 #include <random>
+#include <sstream>
 
 #include "smithlab_utils.hpp"
 #include "smithlab_os.hpp"
@@ -54,6 +55,33 @@ using std::to_string;
 
 using bamxx::bgzf_file;
 
+struct pmd_summary {
+  pmd_summary(const vector<GenomicRegion> &pmds) {
+    pmd_count = pmds.size();
+    pmd_total_size = accumulate(cbegin(pmds), cend(pmds), 0,
+                                [](const uint64_t t, const GenomicRegion &p) {
+                                  return t + p.get_width(); });
+    pmd_mean_size = 
+      static_cast<double>(pmd_total_size)/std::max(1ul, pmd_count);
+  }
+  // pmd_count is the number of identified PMDs.
+  uint64_t pmd_count{};
+  // total_pmd_size is the sum of the sizes of the identified PMDs
+  uint64_t pmd_total_size{};
+  // mean_pmd_size is the mean size of the identified PMDs
+  double pmd_mean_size{};
+
+  string tostring() {
+    std::ostringstream oss;
+    oss << "pmd_count: " << pmd_count << endl
+        << "pmd_total_size: " << pmd_total_size << endl
+        << "pmd_mean_size: "
+        << std::fixed << std::setprecision(2) << pmd_mean_size;
+
+    return oss.str();
+  }
+};
+
 static void
 get_adjacent_distances(const vector<GenomicRegion> &pmds,
                        vector<size_t> &dists) {
@@ -1105,6 +1133,8 @@ main_pmd(int argc, const char **argv) {
     static const double tolerance = 1e-5;
     static const double min_prob  = 1e-10;
 
+    string summary_file;
+
     string params_in_files;
     string params_out_file;
     string posteriors_out_prefix;
@@ -1138,6 +1168,9 @@ main_pmd(int argc, const char **argv) {
     opt_parse.add_opt("posteriors-out", 'r',
                       "write out posterior probabilities in methcounts format",
                       false, posteriors_out_prefix);
+    opt_parse.add_opt("summary", 'S',
+                      "write summary output here",
+                      false, summary_file);
     opt_parse.add_opt("params-out", 'p', "write HMM parameters to this file",
                       false, params_out_file);
     opt_parse.add_opt("seed", 's', "specify random seed",
@@ -1385,6 +1418,12 @@ main_pmd(int argc, const char **argv) {
                         reps_fg_alpha, reps_fg_beta,
                         reps_bg_alpha, reps_bg_beta);
 
+    if (!summary_file.empty()) {
+      ofstream summary_out(summary_file);
+      if (!summary_out) throw runtime_error("failed to open: " + summary_file);
+      summary_out << pmd_summary(good_domains).tostring() << endl;
+    }
+
     ofstream of;
     if (!outfile.empty()) of.open(outfile);
     ostream out(outfile.empty() ? std::cout.rdbuf() : of.rdbuf());

src/utils/uniq.cpp

@@ -77,25 +77,66 @@ struct rd_stats {  // keep track of good bases/reads in and out
   }
 };
 
+
+struct uniq_summary {
+  uniq_summary(const rd_stats &rs_in, const rd_stats &rs_out,
+               const size_t reads_duped) {
+    total_reads = rs_in.reads;
+    total_bases = rs_in.bases;
+    unique_reads = rs_out.reads;
+    unique_read_bases = rs_out.bases;
+    reads_removed = rs_in.reads - rs_out.reads;
+    non_duplicate_fraction = static_cast<double>(rs_out.reads - reads_duped) / 
+                             std::max(1ul, rs_in.reads);
+    duplication_rate = static_cast<double>(reads_removed + reads_duped) / 
+                       std::max(1ul, reads_duped);
+    duplicate_reads = reads_duped;
+  }
+
+  // total_reads is the number of input reads
+  size_t total_reads{};
+  // total_bases is the total number of input bases
+  size_t total_bases{};
+  // unique_reads is the number of unique reads 
+  size_t unique_reads{};
+  // unique_read_bases is the total number of bases for the unique reads
+  size_t unique_read_bases{};
+  // non_duplicate_fraction is the ratio of the number of unique reads with 
+  // no duplicates to that of the input reads 
+  double non_duplicate_fraction{};
+  // duplicate_reads is the number of unique reads with at least one duplicate
+  size_t duplicate_reads{};
+  // reads_removed is the number of duplicate reads that have been removed
+  size_t reads_removed{};
+  // duplication_rate is the average number of duplicates for the reads with 
+  // at least one duplicate (>1 by definition)
+  double duplication_rate{};
+
+  string tostring() {
+    std::ostringstream oss;
+    oss << "total_reads: " << total_reads << endl
+        << "total_bases: " << total_bases << endl
+        << "unique_reads: " << unique_reads << endl
+        << "unique_read_bases: " << unique_read_bases << endl
+        << "non_duplicate_fraction: " << non_duplicate_fraction << endl
+        << "duplicate_reads: " << duplicate_reads << endl
+        << "reads_removed: " << reads_removed << endl
+        << "duplication_rate: " << duplication_rate;
+
+    return oss.str();
+  }
+};
+
+
+
 static void
 write_stats_output(const rd_stats &rs_in, const rd_stats &rs_out,
                    const size_t reads_duped, const string &statfile) {
   if (!statfile.empty()) {
-    const size_t reads_removed = rs_in.reads - rs_out.reads;
-    const double non_dup_frac =
-      (rs_out.reads - reads_duped) / static_cast<double>(rs_in.reads);
-    const double dup_rate =
-      (reads_removed + reads_duped) / static_cast<double>(reads_duped);
+    uniq_summary summary(rs_in, rs_out, reads_duped);
     ofstream out_stat(statfile);
     if (!out_stat) throw runtime_error("bad stats output file");
-    out_stat << "total_reads: " << rs_in.reads << endl
-             << "total_bases: " << rs_in.bases << endl
-             << "unique_reads: " << rs_out.reads << endl
-             << "unique_read_bases: " << rs_out.bases << endl
-             << "non_duplicate_fraction: " << non_dup_frac << endl
-             << "duplicate_reads: " << reads_duped << endl
-             << "reads_removed: " << reads_removed << endl
-             << "duplication_rate: " << dup_rate << endl;
+    out_stat << summary.tostring() << endl;
   }
 }
 
@@ -245,7 +286,7 @@ main_uniq(int argc, const char **argv) {
                            "sorted mapped reads",
                            "<in-file> [out-file]", 2);
     opt_parse.add_opt("threads", 't', "number of threads", false, n_threads);
-    opt_parse.add_opt("stats", 'S', "statistics output file", false, statfile);
+    opt_parse.add_opt("summary", 'S', "statistics output file", false, statfile);
     opt_parse.add_opt("add-count", 'a', "add duplicate counts to reads", false,
                       add_dup_count);
     opt_parse.add_opt("hist", '\0',