reducing excess code, changing noise calculation

Author: cmaceves

src/call_consensus_clustering.cpp

@@ -209,9 +209,6 @@ std::vector<std::vector<float>> find_solutions(std::vector<float> means, float e
   std::sort(results.begin(), results.end());
   results.erase(std::unique(results.begin(), results.end()), results.end());
 
-  auto max_iter = std::max_element(means.begin(), means.end());
-  auto min_iter = std::min_element(means.begin(), means.end());
-
   std::vector<std::vector<float>> final_results;
   //constrain that the solutions must add to 1
   for(uint32_t i=0; i < results.size(); i++){
@@ -341,13 +338,11 @@ std::vector<float> parse_clustering_results(std::string clustering_file){
   }  
   return(numbers);
 }
-void cluster_consensus(std::vector<variant> variants, std::string clustering_file, std::string variants_file, double default_threshold){ 
-  float depth_cutoff = 10; 
+void cluster_consensus(std::vector<variant> variants, std::string clustering_file, std::string variants_file, double default_threshold, uint32_t min_depth, uint8_t min_qual){ 
   double error = 0.10; 
-  float solution_error = 0.05;
-  double quality_threshold = 20; 
+  float solution_error = 0.10;
 
-  double error_rate = cluster_error(variants_file, quality_threshold, depth_cutoff);
+  double error_rate = cluster_error(variants_file, min_qual, min_depth);
   float freq_lower_bound = 1-error_rate;
   float freq_upper_bound = error_rate;
 
@@ -504,7 +499,7 @@ void cluster_consensus(std::vector<variant> variants, std::string clustering_fil
     //TESTLINES
     if(variants[i].nuc.find('+') != std::string::npos) continue;
     //TESTLINES
-    if(variants[i].position == 13572){
+    if(variants[i].position == 22029){
       print = true;
       std::cerr << "\ntop freq " << variants[i].freq << " " << variants[i].nuc << " cluster " << variants[i].cluster_assigned << " " << variants[i].gapped_freq << std::endl;
       std::cerr << "vague assignment " << variants[i].vague_assignment << " del pos " << variants[i].pos_del_flag << " depth flag " << variants[i].depth_flag << std::endl;
@@ -628,7 +623,6 @@ void cluster_consensus(std::vector<variant> variants, std::string clustering_fil
     auto it = std::find(solution.begin(), solution.end(), tmp_mean);
     if(it == solution.end()) continue;
     file << ">"+clustering_file+"_cluster_"+ std::to_string(means[i]) << "\n";
-    std::cerr << all_sequences[i][13571] << std::endl;
     file << all_sequences[i] << "\n";
   }
   file.close(); 

src/call_consensus_clustering.h

@@ -2,6 +2,6 @@
 #ifndef call_consensus_clustering
 #define call_consensus_clustering
 
-void cluster_consensus(std::vector<variant> variants, std::string clustering_file, std::string variants_filename, double default_threshold);
+void cluster_consensus(std::vector<variant> variants, std::string clustering_file, std::string variants_filename, double default_threshold, uint32_t min_depth, uint8_t min_qual);
 void find_combinations(std::vector<float> means, uint32_t index, std::vector<float> &current, std::vector<std::vector<float>> &results);
 #endif

src/estimate_error.cpp

@@ -2,6 +2,13 @@
 #include "gmm.h"
 #include "saga.h"
 
+std::vector<std::vector<double>> assign_nearest_cluster(){
+  /*
+    Given a set of means and variants assign each variant to the nearest cluster.
+  */
+  
+}
+
 std::vector<double> z_score(std::vector<double> data) {
     double mean = std::accumulate(data.begin(), data.end(), 0.0) / data.size();
     double sq_sum = std::inner_product(data.begin(), data.end(), data.begin(), 0.0);
@@ -39,7 +46,7 @@ std::vector<std::vector<uint32_t>> determine_outlier_points(std::vector<double>
     return(removal_points);
 }
 
-double cluster_error(std::string filename, double quality_threshold, uint32_t depth_cutoff){
+double cluster_error(std::string filename, uint8_t quality_threshold, uint32_t depth_cutoff){
   /*
     Here we use clustering to determine the value of the noise.
   */
@@ -76,32 +83,42 @@ double cluster_error(std::string filename, double quality_threshold, uint32_t de
     count_original += 1;
   }
 
-  uint32_t n = 3;
-  gaussian_mixture_model model = retrain_model(n, data_original, variants_original, 2, 0.00001);
-  std::vector<double> means = model.means;
-  
-  //for each cluster this describes the points which are outliers
-  std::vector<std::vector<uint32_t>> removal_points = determine_outlier_points(frequencies, means);
-  std::vector<std::vector<double>> clusters(n);
+  //start with a small n value and if we don't find two major clusters we increase the number of clusters
+  uint32_t n = 2;
+  kmeans_model model;
 
-  for(auto var : variants_original){
-    double tmp = (double)var.gapped_freq;
-    auto it = std::min_element(means.begin(), means.end(), 
-        [tmp](double a, double b) {
-            return std::fabs(a - tmp) < std::fabs(b - tmp);
-        });
-    uint32_t index = std::distance(means.begin(), it);
-    clusters[index].push_back(var.gapped_freq);
+  while(n <= 5){
+    model = train_model(n, data_original);
+    for(auto m : model.means){
+      std::cerr << m << std::endl;
+    }
+    std::vector<double> means = model.means;
+    bool stop=true;
+    for(uint32_t i=0; i < model.clusters.size(); i++){
+      float percent = (float)model.clusters[i].size() / (float)data_original.size();
+      if(percent > 0.85) {
+        n++;
+        stop = false;
+      }
+      std::cerr << "n " << n << " percent " << percent << std::endl;
+    }
+    if(stop) break;
   }
+  //for each cluster this describes the points which are outliers
+  std::vector<std::vector<uint32_t>> removal_points = determine_outlier_points(frequencies, model.means);
+
   uint32_t j = 0;
   uint32_t largest=0;
-  for(uint32_t i=0; i < clusters.size(); i++){
-    if(clusters[i].size() > largest){
+  for(uint32_t i=0; i < model.clusters.size(); i++){
+    //std::cerr << "percent " << (float)model.clusters[i].size() / (float)data_original.size() << std::endl;
+    double mean = std::accumulate(model.clusters[i].begin(), model.clusters[i].end(), 0.0) / model.clusters[i].size();
+    //std::cerr << "mean " << mean << std::endl;
+    if(model.clusters[i].size() > largest){
       j = i;
-      largest = clusters[i].size();
+      largest = model.clusters[i].size();
     }
   }
-  std::vector<double> universal_cluster = clusters[j];
+  std::vector<double> universal_cluster = model.clusters[j];
   std::vector<uint32_t> outliers = removal_points[j];
   std::vector<double> cleaned_cluster;
   for(uint32_t i=0; i < universal_cluster.size(); i++){
@@ -113,5 +130,7 @@ double cluster_error(std::string filename, double quality_threshold, uint32_t de
 
   //get the upper edge of the noise cluster
   auto min_it = std::min_element(cleaned_cluster.begin(), cleaned_cluster.end());
-  return *min_it;
+  std::cerr << *min_it << std::endl;
+  //exit(0);
+  return(*min_it);
 }

src/estimate_error.h

@@ -3,5 +3,5 @@
 #include "./include/armadillo"
 #ifndef estimate_error
 #define estimate_error
-double cluster_error(std::string filename, double quality_threshold, uint32_t depth_cutoff);
+double cluster_error(std::string filename, uint8_t quality_threshold, uint32_t depth_cutoff);
 #endif

src/interval_tree.cpp

@@ -138,8 +138,6 @@ void IntervalTree::detect_primer_issues(ITNode *root, uint32_t find_position){
 
 void IntervalTree::detect_abberations(ITNode *root, uint32_t find_position){
   if (root==NULL) return;
-  //unsure why this line is wrogn, but it is
-  //if (find_position < (uint32_t)root->data->low) return;
     for(uint32_t i=0; i < root->amp_positions.size(); i++){
       if(find_position == root->amp_positions[i].pos){
         test_flux.push_back(root->amp_positions[i]);

src/ivar.cpp

@@ -331,6 +331,8 @@ int main(int argc, char *argv[]) {
   //ivar contam
   if (cmd.compare("contam") == 0) {
     g_args.min_threshold = 0;
+    g_args.min_depth = 10;
+    g_args.min_qual = 20;
     opt = getopt(argc, argv, contam_opt_str);
     while (opt != -1) {
       switch (opt) {
@@ -352,8 +354,8 @@ int main(int argc, char *argv[]) {
       opt = getopt(argc, argv, contam_opt_str);
     }
     if (!g_args.variants.empty() && !g_args.prefix.empty()) {
-      std::vector<variant> variants = gmm_model(g_args.variants, g_args.prefix);
-      cluster_consensus(variants, g_args.prefix, g_args.variants, g_args.min_threshold);
+      std::vector<variant> variants = gmm_model(g_args.variants, g_args.prefix, g_args.min_depth, g_args.min_qual);
+      cluster_consensus(variants, g_args.prefix, g_args.variants, g_args.min_threshold, g_args.min_depth, g_args.min_qual);
     }
     res = 0; 
     g_args.prefix = get_filename_without_extension(g_args.prefix, ".bam");

src/primer_bed.cpp

@@ -1,17 +1,12 @@
 #include "primer_bed.h"
 #include "allele_functions.h"
 #include "ref_seq.h"
-#include <chrono>
 std::string primer::get_name() { return name; }
 
 std::string primer::get_region() { return region; }
 
 int primer::get_score() { return score; }
 
-std::vector<position> primer::get_positions() { return positions; }
-
-void primer::set_positions(position pos) { positions.push_back(pos); }
-
 uint32_t primer::get_start() const { return start; }
 
 uint32_t primer::get_end() const { return end; }
@@ -75,22 +70,6 @@ float weighted_similarity(std::string string_1, std::string string_2){
   return similarity_score / total_weight;
 }
 
-void primer::populate_positions(){
-  for(uint32_t i=0; i < this->hardcoded_length; i++){
-    position tmp;
-    tmp.alleles = populate_basic_alleles();
-    if(this->get_strand() == '+'){
-      tmp.pos = i + this->get_start();
-    } else {
-      if(this->get_end() > i){
-        break;
-      }
-      tmp.pos = this->get_end() - i;
-    }
-    positions.push_back(tmp);
-  }
-}
-
 std::vector<primer> populate_from_file(std::string path, int32_t offset = 0) {
   std::ifstream data(path.c_str());
   std::string line;

src/primer_bed.h

@@ -20,14 +20,10 @@ class primer {
   int16_t pair_indice;
   int16_t indice;
   uint32_t read_count = 0;
-  std::vector<position> positions; //allele counts for this primer
 
  public:
-  uint32_t hardcoded_length=300; //CHANGE
   std::string get_name();
   std::string get_region();
-  std::vector<position> get_positions();
-  void set_positions(position pos);
   int get_score();
   uint32_t get_start() const;
   uint32_t get_end() const;
@@ -36,7 +32,6 @@ class primer {
   int16_t get_pair_indice();
   int16_t get_indice() const;
   uint32_t get_read_count() const;
-  void populate_positions();
   void set_start(uint32_t s);
   void set_end(uint32_t e);
   void set_strand(char s);
@@ -64,8 +59,4 @@ primer get_max_end(std::vector<primer> primers);
 std::vector<uint32_t> get_nlengths();
 std::vector<uint32_t> get_start_positions();
 std::vector<bool> get_direction();
-std::vector<std::vector<std::string>> get_qnames();
-void set_positions(position pos);
-std::vector<position> get_positions();
-void populate_positions();
 #endif

src/saga.cpp

@@ -191,6 +191,9 @@ void merge_reads(const bam1_t* read1, const bam1_t* read2, IntervalTree &amplico
 
   amplicons.find_read_amplicon(start_forward, end_reverse, found_amplicon, bam_get_qname(read1), amp_start, amp_dist);   
   if(!found_amplicon){
+    if(start_forward < 23063 && end_reverse > 23063){
+    std::cerr << start_forward << " " << end_reverse << " " << bam_get_qname(read1) << " " << bam_get_qname(read2) << std::endl;
+    }
     amplicons.add_read_variants(final_positions, final_bases, final_qualities, min_qual);
   } else {
     amplicons.assign_read_amplicon(amp_start, final_positions, final_bases, final_qualities, min_qual);
@@ -402,6 +405,7 @@ int preprocess_reads(std::string bam, std::string bed, std::string bam_out, std:
       start_pos = aln->core.pos;
     }
     //TEST LINES
+    if(start_pos < 22063 || start_pos > 24063) continue;
     //if(start_pos < 15176 || start_pos > 17176) continue;
     bam1_t *r = aln;
     //get the md tag
@@ -519,7 +523,7 @@ int preprocess_reads(std::string bam, std::string bed, std::string bam_out, std:
   std::vector<uint32_t> flagged_positions; 
   std::vector<float> std_deviations;
   std::vector<std::string> pos_nuc;
-  uint32_t test_pos = 0;
+  uint32_t test_pos = 23064;
   //detect fluctuating variants across amplicons
   for(uint32_t i=0; i < amplicons.max_pos; i++){
     amplicons.test_flux.clear();