Bug-fixes and performance improvements for short-read Giraffe's supplementary identification algorithm

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<class AlignmentVector>
-vector<size_t> 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<size_t> 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<size_t> 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<tuple<size_t, size_t, size_t>>& intervals, size_t begin, size_t end) -> vector<size_t> {
+        // 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<tuple<int64_t, int64_t, size_t>>& intervals, int64_t begin, int64_t end, bool* success) -> vector<size_t> {
 
             std::sort(intervals.begin(), intervals.end());  
 
             // map from (end index -> (total coverage, final interval))
-            map<size_t, pair<size_t, size_t>> dp;
+            map<int64_t, pair<size_t, size_t>> dp;
             vector<size_t> backpointer(intervals.size(), numeric_limits<size_t>::max());
 
             for (size_t i = 0; i < intervals.size(); ++i) {
@@ -517,7 +521,7 @@ vector<size_t> 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<int64_t>(max_it->first - get<0>(interval), 0) < it->second.first) {
                         max_it = it;
                     } 
@@ -548,14 +552,26 @@ vector<size_t> identify_supplementaries(const AlignmentVector& alignments, doubl
 
             // traceback the optimum
             vector<size_t> traceback;
-            auto final_it = dp.lower_bound(max<int64_t>(begin, end - max_separation));
-            if (final_it != dp.end()) {
+            auto final_it = dp.end();
+            auto it = dp.lower_bound(max<int64_t>(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<size_t>::max()) {
                     traceback.emplace_back(backpointer[traceback.back()]);
                 }
                 reverse(traceback.begin(), traceback.end());
+
+                *success = true;
+            }
+            else {
+                *success = false;
             }
 
             return traceback;
@@ -566,7 +582,7 @@ vector<size_t> 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<tuple<size_t, size_t, size_t>> left_side, right_side;
+            vector<tuple<int64_t, int64_t, size_t>> left_side, right_side;
 
             for (size_t i = 0; i < alignments.size(); ++i) {
                 if (i == primary_idx) {
@@ -575,30 +591,39 @@ vector<size_t> 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<tuple<size_t, size_t, size_t>> 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
             if (!is_sorted(full_traceback.begin(), full_traceback.end())) {
                 // 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<size_t, size_t> curr_interval(0, 0);
             for (const auto& interval : full_traceback) {
@@ -612,11 +637,10 @@ vector<size_t> 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<int64_t>(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));

src/minimizer_mapper.cpp

@@ -732,7 +732,7 @@ vector<Alignment> MinimizerMapper::map_from_extensions(Alignment& aln) {
             }
             // Could this cluster finish out a supplementary alignment?
             return (total_overlap <= 2 * max_supplementary_separation && 
-                    max<int64_t>(cluster_intervals.total_size() - total_overlap, 0) >= min_supplementary_read_coverage);
+                    max<int64_t>(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<Alignment> MinimizerMapper::map_from_extensions(Alignment& aln) {
             }
             // Could this cluster finish out a supplementary alignment?
             return (total_overlap <= 2 * max_supplementary_separation && 
-                    max<int64_t>(extension_intervals.total_size() - total_overlap, 0) >= min_supplementary_read_coverage);
+                    max<int64_t>(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<Alignment> 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<std::array<vector<Alig
             auto& read_suppl_candidates = candidates[r];
 
             auto supplementary_idxs = identify_supplementaries(read_suppl_candidates, min_supplementary_read_coverage, 
-                                                               max_supplementary_separation, min_supplementary_score_fraction, 
+                                                               max_supplementary_separation, max_supplementary_overlap, 
+                                                               max_supplementary_uncovered_end, min_supplementary_score_fraction, 
                                                                min_supplementary_size, 0);
 
             if (!supplementary_idxs.empty()) {
@@ -3678,12 +3685,12 @@ MinimizerMapper::identify_supplementary_alignments(vector<std::array<vector<Alig
                 // Translate back from vector indexes to alignment indexes
                 auto& read_supplementaries = supplementaries[r];
                 for (auto i : supplementary_idxs) {
-                    if (i <= paired_suppl_source.size()) {
+                    if (i <= paired_suppl_source[r].size()) {
                         const auto& source = paired_alignments[paired_suppl_source[r][i - 1]][r];
                         read_supplementaries.push_back({source.fragment, source.alignment});
                     }
                     else {
-                        const auto& source = unpaired_alignments[unpaired_suppl_source[r][i - paired_suppl_source.size() - 1]];
+                        const auto& source = unpaired_alignments[unpaired_suppl_source[r][i - paired_suppl_source[r].size() - 1]];
                         read_supplementaries.push_back({source.fragment, source.alignment});
                     }
                 }
@@ -3802,7 +3809,8 @@ vector<Alignment> MinimizerMapper::identify_supplementary_alignments(vector<Alig
         }
 
         auto supplementary_idxs = identify_supplementaries(alignments, min_supplementary_read_coverage, 
-                                                           max_supplementary_separation, min_supplementary_score_fraction, 
+                                                           max_supplementary_separation, max_supplementary_overlap, 
+                                                           max_supplementary_uncovered_end, min_supplementary_score_fraction, 
                                                            min_supplementary_size, primary_idx);
 
         if (!supplementary_idxs.empty()) {

src/minimizer_mapper.hpp

@@ -493,18 +493,31 @@ class MinimizerMapper : public AlignerClient {
     static constexpr size_t default_min_supplementary_size = 40;
     size_t min_supplementary_size = default_min_supplementary_size;
 
-    /// The maximum amount of read separation or overlap between supplementary alignment(s) and the primary alignment
-    static constexpr size_t default_max_supplementary_separation = 10;
+    /// Allow clusters to pass through filters if its independent read coverage is this proportion of the min supplementary size
+    static constexpr double default_min_supplementary_filter_size_proportion = 0.75;
+    double min_supplementary_filter_size_proportion = default_min_supplementary_filter_size_proportion;
+
+    /// The maximum amount of read separation between supplementary alignment(s) and the primary alignment
+    static constexpr size_t default_max_supplementary_separation = 40;
     size_t max_supplementary_separation = default_max_supplementary_separation;
 
+    /// The maximum amount of overlap between supplementary alignment(s) and the primary alignment
+    static constexpr size_t default_max_supplementary_overlap = 10;
+    size_t max_supplementary_overlap = default_max_supplementary_overlap;
+
+    /// The maximum amount of each end of the read that can be left unaligned by the collection of primary and supplementary
+    /// alignment(s)
+    static constexpr size_t default_max_supplementary_uncovered_end = 10;
+    size_t max_supplementary_uncovered_end = default_max_supplementary_uncovered_end;
+
     /// The minimum score of a supplementary as a fraction of the primary alignment score
-    static constexpr double default_min_supplementary_score_fraction = 0.4;
-    size_t min_supplementary_score_fraction = default_min_supplementary_score_fraction;
+    static constexpr double default_min_supplementary_score_fraction = 0.0;
+    double min_supplementary_score_fraction = default_min_supplementary_score_fraction;
 
     /// The minimum fraction of the read that the primary and the supplementaries must jointly align in order for
     /// supplementary alignments to be reported from disjoint graph regions
-    static constexpr double default_min_supplementary_read_coverage = 0.9;
-    size_t min_supplementary_read_coverage = default_min_supplementary_read_coverage;
+    static constexpr double default_min_supplementary_read_coverage = 0.8;
+    double min_supplementary_read_coverage = default_min_supplementary_read_coverage;
 
     /// Apply this sample name
     string sample_name;

src/surjecting_alignment_emitter.cpp

@@ -71,7 +71,7 @@ void SurjectingAlignmentEmitter::surject_paired_alignments_in_place(vector<Align
                                mate_info(primary_pos2.name(), primary_pos2.offset(), primary_pos2.is_reverse(), false));
             }
             for (size_t j = 0; j < surjected2.size(); ++j) {
-                if (j == primary_idx1) {
+                if (j == primary_idx2) {
                     continue;
                 }
                 supplementary_alns.emplace_back(std::move(surjected2[j]));