Merge pull request yatisht#334 from yceh/master

automatically disable optimization if ref contain ambiguous nucleotide
Sicheng-Shu · Mar 29, 2023 · 308e929 · 308e929
2 parents 5dc4aae + 3b57a84
commit 308e929
Showing 1 changed file with 18 additions and 2 deletions.
diff --git a/src/usher-sampled/driver/main.cpp b/src/usher-sampled/driver/main.cpp
@@ -155,8 +155,11 @@ void leader_thread_optimization(MAT::Tree& tree,std::vector<mutated_t>& position
                 timeout=true;
                 break;
             }
-            if (new_parsimony_score==last_parsimony_score&&optimization_radius>tree.max_level) {
-                timeout=true;
+            if (new_parsimony_score==last_parsimony_score) {
+                if(optimization_radius>tree.max_level){
+                    timeout=true;
+                }
+                break;
             }
             last_parsimony_score=new_parsimony_score;
         }
@@ -243,6 +246,19 @@ static int leader_thread(
     tree.condense_leaves();
     fix_parent(tree);
     tree.check_leaves();
+    bool have_ambiguous_ref=false;
+    for (int position=1; position<MAT::Mutation::refs.size(); position++) {
+        auto nuc=MAT::Mutation::refs[position];
+        if (nuc&(nuc-1)) {
+            fprintf(stderr, "\nWARNING: Ref nuc @ %d : %c is ambiguous\n", position,MAT::get_nt(MAT::Mutation::refs[position]));
+            have_ambiguous_ref=true;
+        }
+    }
+    if (have_ambiguous_ref) {
+        fprintf(stderr, "WARNING: Reference contain ambiguous nucleotide, optimization is disabled\n");
+        options.initial_optimization_radius=0;
+        optimization_radius=0;
+    }
     if (options.initial_optimization_radius>0) {
         for (auto& pos : position_wise_out) {
             pos.erase(std::remove_if(pos.begin(), pos.end(), [sample_start_idx](const std::pair<long, nuc_one_hot>& in) {