diff --git a/src/alignment.hpp b/src/alignment.hpp index cea9bc524d4..8778d4bc74c 100644 --- a/src/alignment.hpp +++ b/src/alignment.hpp @@ -482,12 +482,15 @@ Alignment target_alignment(const PathPositionHandleGraph* graph, const path_hand /// Returns indexes into an vector-like container of Alignments that correspond to supplementary alignments to the primary /// min_read_coverage Require the supplementaries and primary to cover this fraction of the read -/// max_separation Require the separation or overlap between supplementaries to be at most this many bases +/// max_separation Require the separation between supplementaries to be at most this many bases +/// max_overlap Require the overlap between supplementaries to be at most this many bases +/// max_uncovered_end Require the collection of supplementaries to leave at most this many bases uncovered on each end /// min_score_fraction Require each supplementary to have this fraction of the primary's score /// min_size Require each supplementary to align at least this many bases template -vector identify_supplementaries(const AlignmentVector& alignments, double min_read_coverage, size_t max_separation, - double min_score_fraction, size_t min_size, size_t primary_idx = 0) { +vector identify_supplementaries(const AlignmentVector& alignments, double min_read_coverage, size_t max_separation, + size_t max_overlap, size_t max_uncovered_end, double min_score_fraction, size_t min_size, + size_t primary_idx = 0) { assert(min_read_coverage >= 0.0 && min_read_coverage <= 1.0 && min_score_fraction >= 0.0 && min_score_fraction <= 1.0); @@ -502,12 +505,13 @@ vector identify_supplementaries(const AlignmentVector& alignments, doubl // do sparse DP over part of the read to determine which set of intervals achieves the highest read coverage, subject // to the constraints - auto do_dp = [&](vector>& intervals, size_t begin, size_t end) -> vector { + // the DP assumes that the primary interval is to the left and the end of the sequence is to the right + auto do_dp = [&](vector>& intervals, int64_t begin, int64_t end, bool* success) -> vector { std::sort(intervals.begin(), intervals.end()); // map from (end index -> (total coverage, final interval)) - map> dp; + map> dp; vector backpointer(intervals.size(), numeric_limits::max()); for (size_t i = 0; i < intervals.size(); ++i) { @@ -517,7 +521,7 @@ vector identify_supplementaries(const AlignmentVector& alignments, doubl // look for the best feasible previous DP entry // TODO: this could have better worst-case guarantees with a key-value RMQ auto max_it = dp.end(); - for (auto it = dp.lower_bound(get<0>(interval) - max_separation); it != dp.end() && it->first <= get<0>(interval) + max_separation; ++it) { + for (auto it = dp.lower_bound(get<0>(interval) - max_separation); it != dp.end() && it->first <= get<0>(interval) + max_overlap; ++it) { if (max_it == dp.end() || max_it->second.first - max(max_it->first - get<0>(interval), 0) < it->second.first) { max_it = it; } @@ -548,14 +552,26 @@ vector identify_supplementaries(const AlignmentVector& alignments, doubl // traceback the optimum vector traceback; - auto final_it = dp.lower_bound(max(begin, end - max_separation)); - if (final_it != dp.end()) { + auto final_it = dp.end(); + auto it = dp.lower_bound(max(begin, end - max_uncovered_end)); + while (it != dp.end()) { + if (final_it == dp.end() || it->second.first >= final_it->second.first) { + final_it = it; + } + ++it; + } + if (final_it != dp.end()) { // there is a feasible solution traceback.emplace_back(final_it->second.second); while (backpointer[traceback.back()] != numeric_limits::max()) { traceback.emplace_back(backpointer[traceback.back()]); } reverse(traceback.begin(), traceback.end()); + + *success = true; + } + else { + *success = false; } return traceback; @@ -566,7 +582,7 @@ vector identify_supplementaries(const AlignmentVector& alignments, doubl if (primary_interval.second - primary_interval.first >= min_size && alignments[primary_idx].score() >= min_score) { // records of (begin read pos, end read pos, idx of alignment) - vector> left_side, right_side; + vector> left_side, right_side; for (size_t i = 0; i < alignments.size(); ++i) { if (i == primary_idx) { @@ -575,23 +591,31 @@ vector identify_supplementaries(const AlignmentVector& alignments, doubl auto interval = aligned_interval(alignments[i]); // filter to alignments that meet the minimum thresholds if (alignments[i].score() >= min_score && interval.second - interval.first >= min_size) { - if (interval.second <= primary_interval.first + max_separation) { - left_side.emplace_back(interval.first, interval.second, i); + if (interval.second <= primary_interval.first + max_overlap) { + // negate the positions so that we can iterate left-to-right in the DP + left_side.emplace_back(-interval.second, -interval.first, i); } - else if (interval.first >= primary_interval.second - max_separation) { + else if (interval.first >= primary_interval.second - max_overlap) { right_side.emplace_back(interval.first, interval.second, i); } } } - // do DP on each side and combine the tracebacks + // do DP on each side and combine the tracebacks + + bool left_success = false, right_success = false; vector> full_traceback; - for (auto i : do_dp(left_side, 0, primary_interval.first)) { - full_traceback.push_back(left_side[i]); + // negate the interval of iteration so we can iterate left-to-right + for (auto i : do_dp(left_side, -primary_interval.first, 0, &left_success)) { + auto& interval = left_side[i]; + full_traceback.emplace_back(-get<1>(interval), -get<0>(interval), get<2>(interval)); } + // the negated interval is ordered in reverse, flip it back + std::reverse(full_traceback.begin(), full_traceback.end()); full_traceback.emplace_back(primary_interval.first, primary_interval.second, primary_idx); - for (auto i : do_dp(right_side, primary_interval.second, seq_size)) { - full_traceback.push_back(right_side[i]); + for (auto i : do_dp(right_side, primary_interval.second, seq_size, &right_success)) { + auto interval = right_side[i]; + full_traceback.emplace_back(get<0>(interval), get<1>(interval), get<2>(interval)); } // compute the total read coverage with sweep line algorithm @@ -599,6 +623,7 @@ vector identify_supplementaries(const AlignmentVector& alignments, doubl // TODO: is this even possible? maybe in some weird cases where the max separation is larger than the min size sort(full_traceback.begin(), full_traceback.end()); } + size_t total_cov = 0; pair curr_interval(0, 0); for (const auto& interval : full_traceback) { @@ -612,11 +637,10 @@ vector identify_supplementaries(const AlignmentVector& alignments, doubl } } total_cov += (curr_interval.second - curr_interval.first); - if (aligned_interval(alignments[get<2>(full_traceback.front())]).first <= max_separation && - aligned_interval(alignments[get<2>(full_traceback.back())]).second >= max(seq_size - max_separation, 0) && + if ((left_success || primary_interval.first <= max_uncovered_end) && + (right_success || seq_size - primary_interval.second <= max_uncovered_end) && total_cov >= min_total_cov) { // the supplementaries and primary jointly cover the entire read, the result of DP is feasible - for (auto& interval : full_traceback) { if (get<2>(interval) != primary_idx) { supplementaries.push_back(get<2>(interval)); diff --git a/src/minimizer_mapper.cpp b/src/minimizer_mapper.cpp index 5f0248d348c..56c663933f3 100644 --- a/src/minimizer_mapper.cpp +++ b/src/minimizer_mapper.cpp @@ -732,7 +732,7 @@ vector MinimizerMapper::map_from_extensions(Alignment& aln) { } // Could this cluster finish out a supplementary alignment? return (total_overlap <= 2 * max_supplementary_separation && - max(cluster_intervals.total_size() - total_overlap, 0) >= min_supplementary_read_coverage); + max(cluster_intervals.total_size() - total_overlap, 0) >= min_supplementary_filter_size_proportion * min_supplementary_size); }, cluster_coverage_threshold, min_extensions, max_extensions, rng, [&](size_t cluster_num, size_t item_count, bool escaped_threshold) -> bool { // Handle sufficiently good clusters in descending coverage order @@ -903,7 +903,7 @@ vector MinimizerMapper::map_from_extensions(Alignment& aln) { } // Could this cluster finish out a supplementary alignment? return (total_overlap <= 2 * max_supplementary_separation && - max(extension_intervals.total_size() - total_overlap, 0) >= min_supplementary_read_coverage); + max(extension_intervals.total_size() - total_overlap, 0) >= min_supplementary_filter_size_proportion * min_supplementary_size); }, extension_set_score_threshold, min_extension_sets, max_alignments, rng, [&](size_t extension_num, size_t item_count, bool escaped_threshold) -> bool { // This extension set is good enough. @@ -1000,6 +1000,12 @@ vector MinimizerMapper::map_from_extensions(Alignment& aln) { // Have a function to process the best alignments we obtained auto observe_alignment = [&](Alignment& aln) { + + if (find_supplementaries) { + auto interval = aligned_interval(aln); + current_read_coverage.add(interval.first, interval.second); + } + alignments.emplace_back(std::move(aln)); if (track_provenance) { @@ -3659,7 +3665,8 @@ MinimizerMapper::identify_supplementary_alignments(vector MinimizerMapper::identify_supplementary_alignments(vector{1, 2}); } SECTION("Respects the mininmum size constraint") { min_size = 6; - auto supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + auto supps = identify_supplementaries(alns, read_coverage, separation, separation, separation, score_fraction, min_size); REQUIRE(supps.size() == 2); min_size = 7; - supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + supps = identify_supplementaries(alns, read_coverage, separation, separation, separation, score_fraction, min_size); REQUIRE(supps.empty()); } SECTION("Respects the mininmum score constraint") { score_fraction = 5.0 / 8.0; - auto supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + auto supps = identify_supplementaries(alns, read_coverage, separation, separation, separation, score_fraction, min_size); REQUIRE(supps.size() == 2); score_fraction = 0.7; - supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + supps = identify_supplementaries(alns, read_coverage, separation, separation, separation, score_fraction, min_size); REQUIRE(supps.empty()); } SECTION("Respects read coverage constraint") { read_coverage = 1.0; - auto supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + auto supps = identify_supplementaries(alns, read_coverage, separation, separation, separation, score_fraction, min_size); REQUIRE(supps.size() == 2); auto m = alns[2].mutable_path()->mutable_mapping(0); @@ -664,18 +664,20 @@ TEST_CASE("Supplementary alignments can identified and processed", "[alignment][ e2->set_to_length(2); read_coverage = 0.95; - supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + supps = identify_supplementaries(alns, read_coverage, separation, separation, separation, score_fraction, min_size); REQUIRE(supps.size() == 2); read_coverage = 0.96; - supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + supps = identify_supplementaries(alns, read_coverage, separation, separation, separation, score_fraction, min_size); REQUIRE(supps.size() == 0); } SECTION("Respects separation constraint") { separation = 0; - auto supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + int overlap = 0; + int ends_uncovered = 0; + auto supps = identify_supplementaries(alns, read_coverage, separation, overlap, ends_uncovered, score_fraction, min_size); REQUIRE(supps.size() == 2); auto m = alns[2].mutable_path()->mutable_mapping(0); @@ -688,13 +690,16 @@ TEST_CASE("Supplementary alignments can identified and processed", "[alignment][ e2->set_from_length(2); e2->set_to_length(2); - supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + separation = 0; + supps = identify_supplementaries(alns, read_coverage, separation, overlap, ends_uncovered, score_fraction, min_size); REQUIRE(supps.size() == 0); separation = 1; - supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + supps = identify_supplementaries(alns, read_coverage, separation, overlap, ends_uncovered, score_fraction, min_size); REQUIRE(supps.size() == 2); + separation = 0; + // add a base of overlap e1->set_sequence("CAAGTTCTGCTTT"); e1->set_to_length(13); @@ -702,19 +707,54 @@ TEST_CASE("Supplementary alignments can identified and processed", "[alignment][ e2->set_to_length(4); m->mutable_position()->set_offset(3); - separation = 0; - supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + overlap = 0; + supps = identify_supplementaries(alns, read_coverage, separation, overlap, ends_uncovered, score_fraction, min_size); REQUIRE(supps.size() == 0); - separation = 1; - supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + overlap = 1; + supps = identify_supplementaries(alns, read_coverage, separation, overlap, ends_uncovered, score_fraction, min_size); REQUIRE(supps.size() == 2); + + overlap = 0; + + // return it to the original state e1->set_sequence("CAAGTTCTGCTTTC"); e1->set_to_length(14); e2->set_from_length(3); e2->set_to_length(3); + + // make a soft-clip at the end + { + auto m2 = alns[2].mutable_path()->mutable_mapping(1); + auto e3 = m2->mutable_edit(1); + e3->set_from_length(1); + e3->set_to_length(1); + auto e4 = m2->add_edit(); + e4->set_from_length(0); + e4->set_to_length(1); + e4->set_sequence("T"); + } + + ends_uncovered = 0; + supps = identify_supplementaries(alns, read_coverage, separation, overlap, ends_uncovered, score_fraction, min_size); + REQUIRE(supps.size() == 0); + + ends_uncovered = 1; + supps = identify_supplementaries(alns, read_coverage, separation, overlap, ends_uncovered, score_fraction, min_size); + REQUIRE(supps.size() == 2); + + ends_uncovered = 0; + + // convert back from soft-clip + { + auto m2 = alns[2].mutable_path()->mutable_mapping(1); + auto e3 = m2->mutable_edit(1); + e3->set_from_length(2); + e3->set_to_length(2); + m2->mutable_edit()->RemoveLast(); + } // split it into two alignments @@ -740,11 +780,11 @@ TEST_CASE("Supplementary alignments can identified and processed", "[alignment][ score_fraction = 0.1; separation = 1; - supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + supps = identify_supplementaries(alns, read_coverage, separation, overlap, ends_uncovered, score_fraction, min_size); REQUIRE(supps.size() == 3); separation = 0; - supps = identify_supplementaries(alns, read_coverage, separation, score_fraction, min_size); + supps = identify_supplementaries(alns, read_coverage, separation, overlap, ends_uncovered, score_fraction, min_size); REQUIRE(supps.size() == 0); } }