diff --git a/CMakeLists.txt b/CMakeLists.txt index 147428e..d2d7b0b 100644 --- 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 index 1e40e92..74db138 100644 --- 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 index 26d217b..cb64219 100644 --- a/src/graph/constructor.cpp +++ b/src/graph/constructor.cpp @@ -1692,7 +1692,7 @@ add_var_record(std::vector & 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 index 9240fc0..ea92db1 100644 --- 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, ®ion_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 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 index a6b7587..4f97783 100644 --- 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(qual) / static_cast(depth); - return qd >= 4.0; + return qd >= 3.5; } diff --git a/src/typer/variant.cpp b/src/typer/variant.cpp index 22a3635..e6bad29 100644 --- 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(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 index 3c9484e..9708cae 100644 --- 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 index 0ce21b4..aa4d3e3 100644 --- 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 index c653c93..f151040 100644 --- 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 index fa56d17..d910d6f 100644 Binary files a/test/data/reference/index_test.vcf.gz and b/test/data/reference/index_test.vcf.gz differ diff --git a/test/data/reference/index_test.vcf.gz.tbi b/test/data/reference/index_test.vcf.gz.tbi index 7c95085..a3da4ed 100644 Binary files a/test/data/reference/index_test.vcf.gz.tbi and b/test/data/reference/index_test.vcf.gz.tbi differ diff --git a/test/graph/test_constructor.cpp b/test/graph/test_constructor.cpp index 3ee67f7..3234681 100644 --- 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 const & ref_nodes = graph.ref_nodes; + std::vector 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 index ce87aed..8e2f1da 100644 --- 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));