Merge pull request #122 from smithlabcode/sym-check-sorted

sym: checking that the input is sorted

Author: andrewdavidsmith

src/utils/symmetric-cpgs.cpp

@@ -22,7 +22,8 @@
 #include <iostream>
 #include <fstream>
 #include <stdexcept>
-
+#include <unordered_set>
+#include <filesystem>
 #include <bamxx.hpp>
 
 // from smithlab_cpp
@@ -36,56 +37,70 @@ using std::string;
 using std::cout;
 using std::cerr;
 using std::endl;
+using std::unordered_set;
 
 using bamxx::bgzf_file;
 
-inline bool
+static inline bool
 found_symmetric(const MSite &prev_cpg, const MSite &curr_cpg) {
   // assumes check for CpG already done
-  return (prev_cpg.strand == '+' &&
-          curr_cpg.strand == '-' &&
+  return (prev_cpg.strand == '+' && curr_cpg.strand == '-' &&
           prev_cpg.pos + 1 == curr_cpg.pos);
 }
 
-template <class T>
-static void
+template<class T> static bool
 process_sites(bgzf_file &in, T &out) {
-
   MSite prev_site, curr_site;
   bool prev_is_cpg = false;
   if (read_site(in, prev_site))
-    if (prev_site.is_cpg())
-      prev_is_cpg = true;
+    if (prev_site.is_cpg()) prev_is_cpg = true;
+
+  unordered_set<string> chroms_seen;
+  bool sites_are_sorted = true;
 
   while (read_site(in, curr_site)) {
     if (curr_site.is_cpg()) {
+      if (curr_site.chrom != prev_site.chrom) {
+        const auto chrom_itr = chroms_seen.find(curr_site.chrom);
+        if (chrom_itr != cend(chroms_seen)) {
+          sites_are_sorted = false;
+          break;
+        }
+        else
+          chroms_seen.insert(curr_site.chrom);
+      }
+      else if (curr_site.pos <= prev_site.pos) {
+        sites_are_sorted = false;
+        break;
+      }
       if (prev_is_cpg && found_symmetric(prev_site, curr_site)) {
         prev_site.add(curr_site);
         write_site(out, prev_site);
       }
       prev_is_cpg = true;
     }
-    else prev_is_cpg = false;
+    else
+      prev_is_cpg = false;
     std::swap(prev_site, curr_site);
   }
+
+  return sites_are_sorted;
 }
 
 int
 main_symmetric_cpgs(int argc, const char **argv) {
-
   try {
-
     string outfile;
     // (not used) bool VERBOSE = false;
 
     const string description =
       "Get CpG sites and make methylation levels symmetric.";
 
     /****************** COMMAND LINE OPTIONS ********************/
-    OptionParser opt_parse(strip_path(argv[0]),
-                           description, "<methcounts-file>");
-    opt_parse.add_opt("output", 'o', "output file (default: stdout)",
-                      false, outfile);
+    OptionParser opt_parse(strip_path(argv[0]), description,
+                           "<methcounts-file>");
+    opt_parse.add_opt("output", 'o', "output file (default: stdout)", false,
+                      outfile);
     // opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE);
     std::vector<string> leftover_args;
     opt_parse.parse(argc, argv, leftover_args);
@@ -112,18 +127,27 @@ main_symmetric_cpgs(int argc, const char **argv) {
     bgzf_file in(filename, "r");
     if (!in) throw std::runtime_error("could not open file: " + filename);
 
+    bool sites_are_sorted = true;
     if (outfile.empty() || !has_gz_ext(outfile)) {
       std::ofstream of;
       if (!outfile.empty()) of.open(outfile.c_str());
       std::ostream out(outfile.empty() ? std::cout.rdbuf() : of.rdbuf());
-      process_sites(in, out);
+      sites_are_sorted = process_sites(in, out);
     }
     else {
       bgzf_file out(outfile, "w");
-      process_sites(in, out);
+      sites_are_sorted = process_sites(in, out);
+    }
+
+    if (!sites_are_sorted) {
+      cerr << "sites are not sorted in: " << filename << endl;
+      namespace fs = std::filesystem;
+      const fs::path outpath{outfile};
+      if (fs::exists(outpath)) fs::remove(outpath);
+      return EXIT_FAILURE;
     }
   }
-  catch (const std::runtime_error &e)  {
+  catch (const std::runtime_error &e) {
     cerr << e.what() << endl;
     return EXIT_FAILURE;
   }