Merge branch 'v2.7' into 'master'

v2.7.1 See merge request stat/graphtyper-github!2
DecodeGenetics · Mar 31, 2021 · ffb1ab2 · ffb1ab2
2 parents b1a69e8 + 63a3b38
commit ffb1ab2
Show file tree

Hide file tree

Showing 13 changed files with 61 additions and 33 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -6,7 +6,7 @@ include(ExternalProject)
 # The version number
 set (graphtyper_VERSION_MAJOR 2)
 set (graphtyper_VERSION_MINOR 7)
-set (graphtyper_VERSION_PATCH 0)
+set (graphtyper_VERSION_PATCH 1)
 set(STATIC_DIR "" CACHE STRING "If set, GraphTyper will be built as a static binary using libraries from the given STATIC_DIR.")
 
 # Get the current working branch

diff --git a/include/graphtyper/utilities/options.hpp b/include/graphtyper/utilities/options.hpp
@@ -84,7 +84,7 @@ class Options
   bool is_csi{false};
   bool force_align_both_orientations{false};
   int sam_flag_filter{3840};
-  long max_files_open{576}; // Maximum amount of SAM/BAM/CRAM files can be opened at the same time
+  long max_files_open{864}; // Maximum amount of SAM/BAM/CRAM files can be opened at the same time
   long soft_cap_of_variants_in_100_bp_window{22};
   bool get_sample_names_from_filename{false};
   bool output_all_variants{false};

diff --git a/src/graph/constructor.cpp b/src/graph/constructor.cpp
@@ -1692,7 +1692,7 @@ add_var_record(std::vector<VarRecord> & var_records,
             assert(key == "GT_HAPLOTYPE");
             auto const val_ul = std::stoul(val_str);
             ref.anti_events.insert(val_ul);
-            alt.anti_events.insert(-val_ul);
+            //alt.anti_events.insert(-val_ul);
           }
         } // for (auto const & val : values)
       } // for (auto const & info : infos)

diff --git a/src/typer/caller.cpp b/src/typer/caller.cpp
@@ -47,9 +47,9 @@ namespace
 
 #ifndef NDEBUG
 std::string const debug_read_name = "HISEQ1:33:H9YY4ADXX:1:2110:2792:58362/2";
-long debug_event_pos{699214};
-char debug_event_type{'D'};
-std::size_t debug_event_size{3};
+long debug_event_pos{15001};
+char debug_event_type{'X'};
+std::size_t debug_event_size{1};
 #endif // NDEBUG
 
 
@@ -1017,13 +1017,13 @@ run_first_pass(bam1_t * hts_rec,
         }
         else
         {
-//#ifndef NDEBUG
-//          if (debug_event_type == 'X' && snp.pos == debug_event_pos)
-//          {
-//            BOOST_LOG_TRIVIAL(info) << __HERE__ << " debug snp has bad support " << snp.to_string() << " "
-//                                    << info.to_string() << " file_i=" << file_i;
-//          }
-//#endif // DEBUG
+#ifndef NDEBUG
+          if (debug_event_type == 'X' && snp.pos == debug_event_pos)
+          {
+            BOOST_LOG_TRIVIAL(info) << __HERE__ << " debug snp has bad support " << snp.to_string() << " "
+                                    << info.to_string() << " file_i=" << file_i;
+          }
+#endif // DEBUG
 
           snp_it = bucket.events.erase(snp_it);
         }
@@ -1233,9 +1233,10 @@ run_first_pass(bam1_t * hts_rec,
       if (support_ratio < 0.3)
         support_ratio = 0.3;
 
-      // 0 low coverage or ambigous
+      // 0 low coverage
       // 1 support
       // 2 anti support
+      // 3 ambigous
       auto is_good_support =
         [&cov_down, &region_begin, &support_ratio]
           (long local_cov,
@@ -1285,7 +1286,7 @@ run_first_pass(bam1_t * hts_rec,
             return IS_ANY_HAP_SUPPORT;
           }
 
-          return 0;
+          return IS_ANY_ANTI_HAP_SUPPORT | IS_ANY_HAP_SUPPORT;
         };
 
       // check this bucket
@@ -2179,9 +2180,9 @@ realign_to_indels(std::vector<Tindel_events::iterator> const & realignment_indel
 
     double const count = correction * (indel_info.hq_count + indel_info.lq_count);
 
-    bool const is_good_count = (indel_info.hq_count >= 5 && count >= 8.0) ||
-                               (indel_info.span >= 9 && indel_info.hq_count >= 5 && count >= 7.0) ||
-                               (indel_info.span >= 18 && indel_info.hq_count >= 4 && count >= 6.0);
+    bool const is_good_count = (indel_info.hq_count >= 5 && count >= 5.5) ||
+                               (indel_info.span >= 5 && indel_info.hq_count >= 4 && count >= 5.0) ||
+                               (indel_info.span >= 15 && indel_info.hq_count >= 3 && count >= 4.5);
 
 #ifndef NDEBUG
     if (debug_event_type == indel.type &&

diff --git a/src/typer/event.cpp b/src/typer/event.cpp
@@ -321,11 +321,11 @@ EventSupport::is_good_indel(uint32_t eps) const
   assert(sequence_reversed <= depth);
   long const qual = 3 * get_log_qual(hq_count + lq_count, anti_count, eps);
 
-  if (qual < 60)
+  if (qual < 50)
     return false;
 
   double const qd = static_cast<double>(qual) / static_cast<double>(depth);
-  return qd >= 4.0;
+  return qd >= 3.5;
 }
 
 

diff --git a/src/typer/variant.cpp b/src/typer/variant.cpp
@@ -898,8 +898,8 @@ Variant::generate_infos()
 
         if (per_al.total_depth > 0 &&
             qd > 0.1 &&
-            per_al.maximum_alt_support >= 3 &&
-            per_al.maximum_alt_support_ratio >= 0.175)
+            per_al.maximum_alt_support >= 2 &&
+            per_al.maximum_alt_support_ratio >= 0.15)
         {
           double const alt_seq_depth = per_al.total_depth;
           double const _sb = 2.0 * ((static_cast<double>(sb_alt[s]) / alt_seq_depth) - 0.5);
@@ -1009,14 +1009,14 @@ Variant::generate_infos()
 #ifndef NDEBUG
       if (is_good_alt[a] == 0)
       {
-        BOOST_LOG_TRIVIAL(debug) << __HERE__ << " In variant=" << to_string(true) // skip calls
-                                 << " bad alt="
-                                 << std::string(seqs[a + 1].begin(), seqs[a + 1].end())
-                                 << " MaxAAS=" << per_al.maximum_alt_support
-                                 << " MaxAASR=" << per_al.maximum_alt_support_ratio
-                                 << " AAScore=" << aa_score[a]
-                                 << " ABHom=" << info_abhom
-                                 << " QDAlt=" << qd_alt[a];
+        BOOST_LOG_TRIVIAL(info) << __HERE__ << " In variant=" << to_string(true) // skip calls
+                                << " bad alt="
+                                << std::string(seqs[a + 1].begin(), seqs[a + 1].end())
+                                << " MaxAAS=" << per_al.maximum_alt_support
+                                << " MaxAASR=" << per_al.maximum_alt_support_ratio
+                                << " AAScore=" << aa_score[a]
+                                << " ABHom=" << info_abhom
+                                << " QDAlt=" << qd_alt[a];
       }
 #endif // NDEBUG
     }

diff --git a/test/data/reference/index_test.fa b/test/data/reference/index_test.fa
@@ -30,3 +30,6 @@ GGGGGGGGG
 >chr10
 GGGGGGGGGGGTGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
 GGGGGGGGG
+>chr11
+AAGGGCGGGGGTGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
+GGGGGGGGG
diff --git a/test/data/reference/index_test.fa.fai b/test/data/reference/index_test.fa.fai
@@ -8,3 +8,4 @@ chr7	280	878	70	71
 chr8	79	1168	70	71
 chr9	79	1255	70	71
 chr10	79	1343	70	71
+chr11	79	1431	70	71
diff --git a/test/data/reference/index_test.vcf b/test/data/reference/index_test.vcf
@@ -15,3 +15,7 @@ chr9	5	.	G	A	0	.	GT_ANTI_HAPLOTYPE=2;GT_ID=1
 chr9	10	.	G	GA	0	.	GT_ID=2
 chr10	5	.	G	A	0	.	GT_HAPLOTYPE=2;GT_ID=1
 chr10	10	.	G	GA	0	.	GT_ID=2
+chr11	2	.	AGGGC	A	0	.	GT_ANTI_HAPLOTYPE=2,3,4;GT_ID=1
+chr11	4	.	G	T	0	.	GT_ANTI_HAPLOTYPE=4;GT_HAPLOTYPE=3;GT_ID=2
+chr11	5	.	G	T	0	.	GT_ANTI_HAPLOTYPE=4;GT_ID=3
+chr11	6	.	C	T	0	.	GT_ID=4
diff --git a/test/data/reference/index_test.vcf.gz b/test/data/reference/index_test.vcf.gz
diff --git a/test/data/reference/index_test.vcf.gz.tbi b/test/data/reference/index_test.vcf.gz.tbi
diff --git a/test/graph/test_constructor.cpp b/test/graph/test_constructor.cpp
@@ -533,11 +533,30 @@ TEST_CASE("Construct test graph with anti events (chr10)")
   {
     REQUIRE(var_nodes[0].anti_events.size() == 1);
     REQUIRE(var_nodes[0].anti_events.count(2) == 1);
-    REQUIRE(var_nodes[1].anti_events.size() == 1);
-    REQUIRE(var_nodes[1].anti_events.count(-2) == 1);
+    REQUIRE(var_nodes[1].anti_events.size() == 0);
     REQUIRE(var_nodes[2].anti_events.size() == 0);
     REQUIRE(var_nodes[3].anti_events.size() == 0);
   }
 
   Options::instance()->add_all_variants = false;
 }
+
+
+TEST_CASE("Construct test graph with anti events (chr11)")
+{
+  using namespace gyper;
+  Options::instance()->add_all_variants = true;
+
+  BOOST_LOG_TRIVIAL(debug) << __HERE__ << " TEST CASE: Construct test graph (chr11).";
+  create_test_graph("/test/data/reference/index_test.fa", "/test/data/reference/index_test.vcf.gz", "chr11", true);
+
+  std::vector<gyper::RefNode> const & ref_nodes = graph.ref_nodes;
+  std::vector<gyper::VarNode> const & var_nodes = graph.var_nodes;
+
+  {
+    REQUIRE(ref_nodes.size() == 2);
+    REQUIRE(var_nodes.size() == 6);
+  }
+
+  Options::instance()->add_all_variants = false;
+}
diff --git a/test/index/test_index.cpp b/test/index/test_index.cpp
@@ -455,7 +455,7 @@ TEST_CASE("Test index chr10 with parity event")
   REQUIRE(labels[0].variant_id != labels[1].variant_id);
 
   labels = ph_index.get(gyper::to_uint64("AGGGGGGTGGGGGGGGGGGGGGGGGGGGGGGG", 0));
-  REQUIRE(labels.size() == 0);
+  REQUIRE(labels.size() == 2);
 
   labels = ph_index.get(gyper::to_uint64("AGGGGGAGTGGGGGGGGGGGGGGGGGGGGGGG", 0));