Beta 18 (#60)

* init commit

* Added an FDR estimator (--fdr) and improved the behvaior of --tune_indel

Author: DanielOlson

src/cli.cpp

@@ -100,6 +100,10 @@ void Settings::prepare_settings() {
                "Do not output BED or JSON annotation")
       ->group("Output");
 
+  app.add_flag("--fdr", this->estimate_fdr,
+               "Estimate the False Discovery rate (runtime will be twice as long)")
+      ->group("Output");
+
   // *************
   // Mask options
   // *************
@@ -480,7 +484,7 @@ bool Settings::parse_input(int argc, const char **argv) {
     passed = false;
   }
 
-  if (this->tune_only || this->tune_medium || this->tune_large) {
+  if (this->tune_only || this->tune_medium || this->tune_large || this->tune_indels) {
     this->tune = true;
   }
 

src/cli.hpp

@@ -5,7 +5,7 @@
 #ifndef ULTRA_CLI_HPP
 #define ULTRA_CLI_HPP
 
-#define ULTRA_VERSION_STRING "1.0.0 (beta 17)"
+#define ULTRA_VERSION_STRING "1.0.0 (beta 18)"
 
 
 #include "../lib/CLI11.hpp"
@@ -26,6 +26,7 @@ struct Settings {
   float p_value_loc = 4.27294921875;
   float p_value_scale = 1.8913602828979492;
   float p_value_freq = 0.000330256;
+  bool estimate_fdr = false;
 
   bool json = false;
   bool hide_seq = false;

src/main.cpp

@@ -133,10 +133,13 @@ int main(int argc, const char *argv[]) {
   }
 
   // Perform the actual run
+
+
   auto *ultra = new Ultra(settings);
   ultra->AnalyzeFile();
   ultra->OutputRepeats(true);
-
+  unsigned long long true_coverage = ultra->Coverage();
+  unsigned long long shuff_coverage = 0;
   // Check to see if we are making a masked file
   if (settings->produce_mask) {
     FILE *f = fopen(settings->mask_file.c_str(), "w");
@@ -146,6 +149,21 @@ int main(int argc, const char *argv[]) {
   }
 
   delete ultra;
+  if (settings->estimate_fdr) {
+    settings->suppress_out = true;
+    settings->hide_settings = true;
+    settings->produce_mask = true;
+    settings->no_split = true;
+    settings->max_split = 0;
+    settings->json = false;
+    ultra = new Ultra(settings);
+    ultra->shuffleSequence = true;
+    ultra->AnalyzeFile();
+    ultra->OutputRepeats(true);
+    shuff_coverage = ultra->Coverage();
+    float fdr = (float)shuff_coverage / (float)true_coverage;
+    printf("Estimated false discovery rate: %g = (%llu / %llu)\n", fdr, shuff_coverage, true_coverage);
+  }
 
   return 0;
 }

src/ultra.cpp

@@ -297,6 +297,20 @@ void Ultra::OutputULTRASettings() {
 }
 
 void Ultra::OutputRepeat(RepeatRegion *r, bool isSubRep) {
+
+  if (r->readID > last_read_id) {
+    last_rep_end = 0;
+  }
+
+  auto rep_start = r->sequenceStart;
+  if (rep_start < last_rep_end) {
+    rep_start = last_rep_end;
+  }
+
+  if (rep_start < r->sequenceStart + r->repeatLength) {
+    total_coverage += (r->sequenceStart + r->repeatLength) - rep_start;
+  }
+
   if (!settings->suppress_out)
     writer->WriteRepeat(r);
   if (settings->produce_mask) {
@@ -316,23 +330,7 @@ void Ultra::SortRepeatRegions() {
 }
 
 unsigned long long Ultra::Coverage() {
-  unsigned long long coverage = 0;
-
-  // Itterate through each sequence and find the coverage for that sequence
-  for (auto tup : this->masks_for_seq) {
-    auto seq_id = tup.first;
-    auto regions = tup.second;
-    auto cleaned_regions = CleanedMasks(regions);
-
-    for (auto region : *cleaned_regions) {
-      coverage += region.end - region.start;
-    }
-
-    cleaned_regions->clear();
-    delete cleaned_regions;
-
-  }
-  return coverage;
+    return total_coverage;
 }
 
 Ultra::Ultra(Settings *s) {

src/ultra.hpp

@@ -67,6 +67,10 @@ class Ultra {
   bool firstRepeat = true;
   unsigned long repeatBuffer = 2000;
 
+  unsigned long long total_coverage;
+  int last_read_id = 0;
+  unsigned long long last_rep_end = 0;
+
   bool storeTraceAndSequence = false;