update

Author: cmaceves

src/call_consensus_clustering.cpp

@@ -1,6 +1,7 @@
 #include "call_consensus_clustering.h"
 #include "estimate_error.h"
 #include "gmm.h"
+#include "saga.h"
 #include <ostream>
 #include <iostream>
 #include <vector>
@@ -9,6 +10,33 @@
 #include <algorithm>
 #include <numeric>
 
+bool test_cluster_deviation(float nearest_cluster, float variant_cluster, float std_dev){
+  bool fluctuation = false;
+  //CLEANUP THIS CAN BE CALCULATED ONCE PER ALL CLUSTERS
+  //determine if the assigned and nearest cluster can be resolved based on variant fluctuation
+  std::vector<double> tmp = {(double) nearest_cluster, (double) variant_cluster};
+  double cluster_dev = calculate_standard_deviation(tmp);
+  if((double)std_dev > cluster_dev){
+    fluctuation = true;
+  }
+  return(fluctuation);
+}
+
+double find_neighboring_cluster(double freq, uint32_t cluster_assigned, std::vector<double> means){
+  //find closest cluster by mean value
+  double min_dist = 1;
+  uint32_t index = 0;
+  for(uint32_t i=0; i < means.size(); i++){
+    if(i == cluster_assigned) continue;
+    double dist = std::abs(means[i]-freq);
+    if(dist < min_dist){
+      min_dist = dist;
+      index = i;
+    }
+  }
+  return(means[index]);
+}
+
 void call_majority_consensus(std::vector<variant> variants, uint32_t max_position, std::string clustering_file, double default_threshold){
   //if we can't find a solution simply take the majority variant per position
   std::vector<std::string> nucs;
@@ -229,12 +257,15 @@ std::vector<std::vector<uint32_t>> find_combination_peaks(std::vector<float> sol
   //given a solution and the means, map each cluster to the cluster it contains
   for(uint32_t i=0; i < means.size(); i++){
     auto it = std::find(solution.begin(), solution.end(), means[i]);
+    //th mean is part of the solution
     if(it != solution.end()){
         cluster_indexes[i].push_back(i);
         float target = means[i];
         std::vector<float> distances(totals.size());
         std::transform(totals.begin(), totals.end(), distances.begin(), [target](float num) { return std::abs(target - num); }); 
         uint32_t count = 0;
+        
+        //this checks the distances from the mean to all other possible peaks
         for(uint32_t d=0; d < distances.size(); d++){
           if(distances[d] < 0.03) count += 1;
         }
@@ -313,9 +344,10 @@ std::vector<float> parse_clustering_results(std::string clustering_file){
 void cluster_consensus(std::vector<variant> variants, std::string clustering_file, std::string variants_file, double default_threshold){ 
   float depth_cutoff = 10; 
   double error = 0.10; 
-  float solution_error = 0.10;
- 
-  std::vector<float> error_rate = cluster_error(variants_file);
+  float solution_error = 0.05;
+  double quality_threshold = 20; 
+
+  std::vector<float> error_rate = cluster_error(variants_file, quality_threshold);
   float freq_lower_bound = error_rate[0];
   float freq_upper_bound = error_rate[1];
 
@@ -340,7 +372,7 @@ void cluster_consensus(std::vector<variant> variants, std::string clustering_fil
   //find position wise frequency pairs
   std::vector<std::vector<float>> pairs = frequency_pair_finder(variants, freq_lower_bound, freq_upper_bound, means); 
   std::vector<std::vector<float>> solutions = find_solutions(means, error);  
-
+  
   //find peaks that can't be a subset of other peaks
   std::vector<float> non_subset_means;
   for(uint32_t i=0; i < means.size(); i++){
@@ -358,7 +390,6 @@ void cluster_consensus(std::vector<variant> variants, std::string clustering_fil
         realistic_solutions.push_back(solutions[i]);
      }
   }
-
   //check each solution that every possible peak is accounted for
   std::vector<std::vector<float>> solution_sets;
   for(uint32_t i=0; i < realistic_solutions.size(); i++){
@@ -367,6 +398,15 @@ void cluster_consensus(std::vector<variant> variants, std::string clustering_fil
       solution_sets.push_back(realistic_solutions[i]);
     }
   }
+
+  for(auto sol : solution_sets){
+    std::cerr << "\nsolution" << std::endl;
+    for(auto s : sol){
+      std::cerr << s << " ";
+    }
+  }
+  std::cerr << "\n" << std::endl;
+
   std::vector<float> solution;
   bool traditional_majority= false; //if we can't find a solution call a traditional majority consensus
   if(solution_sets.size() == 0){
@@ -389,11 +429,10 @@ void cluster_consensus(std::vector<variant> variants, std::string clustering_fil
   }
   std::vector<float> unresolved;
   std::vector<std::vector<uint32_t>> cluster_groups = find_combination_peaks(solution, means, unresolved);
-  for(auto u : unresolved) std::cerr << "unresolved " << u << std::endl;
+  
   std::vector<std::vector<uint32_t>> inverse_groups(means.size());
   for(uint32_t i=0; i < cluster_groups.size(); i++){
     for(uint32_t j=0; j < cluster_groups[i].size(); j++){
-      //std::cerr << cluster_groups[i][j] << " i " << i << " j " << j << std::endl; 
       inverse_groups[cluster_groups[i][j]].push_back(i);
     }
   }
@@ -465,7 +504,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 == 0){
+    if(variants[i].position == 13572){
       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;
@@ -482,11 +521,33 @@ void cluster_consensus(std::vector<variant> variants, std::string clustering_fil
       continue;
     }
     //if this position is experiencing fluctuation across amplicons, call ambiguity
-    if(variants[i].amplicon_flux && variants[i].freq < freq_upper_bound){
+    if(variants[i].amplicon_flux && variants[i].freq < freq_upper_bound && variants[i].freq > freq_lower_bound){
+
+      //find all clusters not part of the same assignment
+      std::vector<double> other_population_clusters;
+      for(uint32_t j=0; j < inverse_groups.size(); j++){
+        //check to make sure you're lookin at a group that's part of the solution
+        auto mit = std::find(solution.begin(), solution.end(), means[j]);
+        if(mit == solution.end()) continue;
+        auto it = std::find(inverse_groups[j].begin(), inverse_groups[j].end(), variants[i].cluster_assigned);      
+        //assigned cluster is not apart of the population
+        if(it == inverse_groups[j].end())
+        for(auto ig : inverse_groups[j]){
+          //CLEAN UP this will push redundant things back
+          other_population_clusters.push_back(means[ig]);
+        }
+      }
+
+      //find the second closest cluster index
+      double closest_mean = find_neighboring_cluster((double)variants[i].gapped_freq, variants[i].cluster_assigned, other_population_clusters);
+      //check if the cluster is within the standard dev of the variant
+      bool fluctuating = test_cluster_deviation(closest_mean, means[variants[i].cluster_assigned], variants[i].std_dev);       
       if(print){
-        std::cerr << "position is experiencing fluctuation" << std::endl;
+        std::cerr << "fluctuating " << fluctuating << std::endl;
+      }
+      if(fluctuating){
+        continue;
       }
-      continue;
     }
 
     //if this amplicon is experiencing fluctuation across amplicons, call ambiguity
@@ -541,6 +602,7 @@ void cluster_consensus(std::vector<variant> variants, std::string clustering_fil
       }
       continue;
     }
+
     for(uint32_t j=0; j < inverse_groups.size(); j++){
       //check to make sure you're lookin at a group that's part of the solution
       auto mit = std::find(solution.begin(), solution.end(), means[j]);
@@ -566,6 +628,7 @@ 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/estimate_error.cpp

@@ -1,95 +1,122 @@
 #include "estimate_error.h"
 #include "gmm.h"
+#include "saga.h"
 
-std::vector<float> cluster_error(std::string filename){
+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);
+    double stddev = std::sqrt(sq_sum / data.size() - mean * mean);
+
+    std::vector<double> z_scores;
+    for (double x : data)
+        z_scores.push_back((x - mean) / stddev);
+    return z_scores;
+}
+
+std::vector<std::vector<uint32_t>> determine_outlier_points(std::vector<double> data, std::vector<double> means){
+    std::vector<std::vector<double>> clusters(means.size());
+    //assign points to clusters
+    for(auto tmp : data){
+      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(tmp);
+    }
+
+    std::vector<std::vector<uint32_t>> removal_points(means.size());
+    //calculate cluster specific percentiles
+    for(uint32_t j=0; j < clusters.size(); j++){
+      std::vector<double> z_scores = z_score(clusters[j]);
+      for(uint32_t i=0; i < z_scores.size(); i++){
+        double abs = std::abs(z_scores[i]);
+        if(abs >= 10){
+          removal_points[j].push_back(i);
+        }
+      }
+    }
+    return(removal_points);
+}
+
+std::vector<float> cluster_error(std::string filename, double quality_threshold){
   /*
     Here we use clustering to determine the value of the noise.
   */
 
-  float lower_bound = 0.01;
-  float upper_bound = 0.10;
+  float lower_bound = 0.50;
+  float upper_bound = 0.99;
   uint32_t depth_cutoff = 10;
-  float quality_threshold = 20;
   uint32_t round_val = 4;
 
   std::vector<uint32_t> deletion_positions = find_deletion_positions(filename, depth_cutoff, lower_bound, upper_bound, round_val);
   std::vector<variant> base_variants;
-  parse_internal_variants(filename, base_variants, depth_cutoff, lower_bound, upper_bound, deletion_positions, round_val);
+  parse_internal_variants(filename, base_variants, depth_cutoff, lower_bound, upper_bound, deletion_positions, round_val, quality_threshold);
   if(base_variants[0].version_1_var){
     calculate_reference_frequency(base_variants, filename, depth_cutoff, lower_bound, upper_bound, deletion_positions);
   }
 
-  gaussian_mixture_model retrained_original;
   std::vector<variant> variants_original;
   uint32_t useful_count_original = 0;
   uint32_t max_pos = 0;
+  std::vector<double> frequencies;
    for(uint32_t i=0; i < base_variants.size(); i++){
     if(base_variants[i].position > max_pos) max_pos = base_variants[i].position;
-    //changed this for TEST
-    if(!base_variants[i].amplicon_flux && !base_variants[i].depth_flag && !base_variants[i].outside_freq_range && !base_variants[i].qual_flag && !base_variants[i].del_flag && !base_variants[i].amplicon_masked && !base_variants[i].primer_masked){
-      
+    if(!base_variants[i].amplicon_flux && !base_variants[i].depth_flag && !base_variants[i].outside_freq_range && !base_variants[i].qual_flag && !base_variants[i].del_flag && !base_variants[i].amplicon_masked && !base_variants[i].primer_masked){      
       useful_count_original++;
       variants_original.push_back(base_variants[i]);
+      frequencies.push_back(base_variants[i].freq);
     }
   }
+
   arma::mat data_original(1, useful_count_original, arma::fill::zeros);
-  //std::cerr << useful_count_original << std::endl;
   uint32_t count_original=0;
   for(uint32_t i = 0; i < variants_original.size(); i++){
-    //check if variant should be filtered for first pass model
-    double tmp = static_cast<double>(variants_original[i].freq); //transform
+    double tmp = static_cast<double>(variants_original[i].gapped_freq);
     data_original.col(count_original) = tmp;
     count_original += 1;
   }
-  uint32_t n = 3;
-  retrained_original = retrain_model(n, data_original, variants_original, 2, 0.0001);
-  assign_clusters(variants_original, retrained_original, 2);
+
+  uint32_t n = 2;
+  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);
+
   for(auto var : variants_original){
-    int cluster = var.cluster_assigned;
-    clusters[cluster].push_back(var.freq);
-    //std::cerr << var.freq << " " << cluster << std::endl;
+    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);
   }
   uint32_t j = 0;
   uint32_t largest=0;
-  std::vector<float> cluster_percents;
   for(uint32_t i=0; i < clusters.size(); i++){
-    std::cerr << clusters[i].size() << std::endl;
-    cluster_percents.push_back((float)clusters[i].size() / (float)variants_original.size());
     if(clusters[i].size() > largest){
       j = i;
       largest = clusters[i].size();
     }
   }
-  for(auto me : retrained_original.means){
-    std::cerr << me << std::endl;
-  }  
-  for(auto per : cluster_percents){
-    std::cerr << per << std::endl;
-  }
-  //get the upper edge of the noise cluster
-  auto max_it = std::max_element(clusters[j].begin(), clusters[j].end());
-  std::cerr << *max_it << std::endl;
-  std::vector<float> max_error_pos(max_pos+1);
-  for(uint32_t i=0; i < base_variants.size(); i++){
-    if(base_variants[i].freq <= *max_it){
-      max_error_pos[base_variants[i].position] += base_variants[i].freq;
+  std::vector<double> universal_cluster = 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++){
+    auto it = std::find(outliers.begin(), outliers.end(), i);
+    if(it == outliers.end()){
+      cleaned_cluster.push_back(universal_cluster[i]);
     }
   }
 
-  float max_error = 0;  
-  for(auto x : max_error_pos){
-    if(x >= *max_it && x > max_error){
-      max_error = x;
-    }
-  }
-  //std::cerr << "max error " << max_error << std::endl;
+  //get the upper edge of the noise cluster
+  auto min_it = std::min_element(cleaned_cluster.begin(), cleaned_cluster.end());
+  std::cerr << "min it " << *min_it << std::endl;
   std::vector<float> tmp;
-  //tmp.push_back(*max_it);
-  //tmp.push_back(1-*max_it);
-
-  tmp.push_back(max_error);
-  tmp.push_back(1-max_error);
-  //exit(0);
+  tmp.push_back(1-(*min_it));
+  tmp.push_back(*min_it);
   return tmp;
 }

src/estimate_error.h

@@ -3,5 +3,5 @@
 #include "./include/armadillo"
 #ifndef estimate_error
 #define estimate_error
-std::vector<float> cluster_error(std::string filename);
+std::vector<float> cluster_error(std::string filename, double quality_threshold);
 #endif