Merge remote-tracking branch 'origin/master' into deconstruct

Author: glennhickey

src/algorithms/chain_items.cpp

@@ -148,15 +148,30 @@ class Anchor {
 
     /// Compose a read start position, graph start position, and match length into an Anchor.
     /// Can also bring along a distance hint and a seed number.
-    inline Anchor(size_t read_start, const pos_t& graph_start, size_t length, size_t margin_before, size_t margin_after, int score, size_t seed_number = std::numeric_limits<size_t>::max(), ZipCode* hint = nullptr, size_t hint_start = 0, bool skippable = false) : start(read_start), size(length), margin_before(margin_before), margin_after(margin_after), start_pos(graph_start), end_pos(advance(graph_start, length)), points(score), start_seed(seed_number), end_seed(seed_number), start_zip(hint), end_zip(hint), start_offset(hint_start), end_offset(length - hint_start), seed_length(margin_before + length + margin_after), skippable(skippable) {
+        inline Anchor(size_t read_start, const pos_t& graph_start, size_t length, size_t margin_before, size_t margin_after,
+        int score, size_t seed_number = std::numeric_limits<size_t>::max(), 
+        ZipCode* hint = nullptr, size_t hint_start = 0, bool skippable = false) : 
+        start(read_start), size(length), margin_before(margin_before), margin_after(margin_after), 
+        start_pos(graph_start), end_pos(advance(graph_start, length)), points(score), 
+        start_seed(seed_number), end_seed(seed_number), start_zip(hint), end_zip(hint), 
+        start_offset(hint_start), end_offset(length - hint_start), 
+        seed_length(margin_before + length + margin_after), skippable(skippable) {
         // Nothing to do!
     }
 
     /// Compose two Anchors into an Anchor that represents coming in through
     /// the first one and going out through the second, like a tunnel. Useful
     /// for representing chains as chainable items.
-    inline Anchor(const Anchor& first, const Anchor& last, size_t extra_margin_before, size_t extra_margin_after, int score) : start(first.read_start()), size(last.read_end() - first.read_start()), margin_before(first.margin_before + extra_margin_before), margin_after(last.margin_after + extra_margin_after), start_pos(first.graph_start()), end_pos(last.graph_end()), points(score), start_seed(first.seed_start()), end_seed(last.seed_end()), start_zip(first.start_hint()), end_zip(last.end_hint()), start_offset(first.start_offset), end_offset(last.end_offset), seed_length((first.base_seed_length() + last.base_seed_length()) / 2), skippable(first.is_skippable() || last.is_skippable()) {
-        // Nothing to do!
+    inline Anchor(const Anchor& first, const Anchor& last, 
+                  size_t extra_margin_before, size_t extra_margin_after, int score) : 
+                  start(first.read_start()), size(last.read_end() - first.read_start()), 
+                  margin_before(first.margin_before + extra_margin_before), 
+                  margin_after(last.margin_after + extra_margin_after), start_pos(first.graph_start()), 
+                  end_pos(last.graph_end()), points(score), start_seed(first.seed_start()), end_seed(last.seed_end()), 
+                  start_zip(first.start_hint()), end_zip(last.end_hint()), 
+                  start_offset(first.start_offset), end_offset(last.end_offset), 
+                  seed_length((first.base_seed_length() + last.base_seed_length()) / 2), 
+                  skippable(first.is_skippable() || last.is_skippable()) {
     }
 
     // Act like data
@@ -319,6 +334,37 @@ transition_iterator zip_tree_transition_iterator(const std::vector<SnarlDistance
                                                  size_t max_graph_lookback_bases,
                                                  size_t max_read_lookback_bases);
 
+/**
+ * Walk through the ziptree from left to right to find all seeds,
+ * then generate all possible transitions from each given seed.
+ * 
+ * Calls ZipCodeTree.find_distances() as the core of the algorithm.
+ * Used as a helper by zip_tree_transition_iterator().
+ * 
+ * Transitions are (source anchor, destination anchor, graph distance).
+ */
+std::vector<std::tuple<size_t, size_t, size_t>> generate_zip_tree_transitions(
+    const std::vector<SnarlDistanceIndexClusterer::Seed>& seeds,
+    const ZipCodeTree& zip_code_tree,
+    size_t max_graph_lookback_bases,
+    const std::unordered_map<size_t, size_t>& seed_to_starting, 
+    const std::unordered_map<size_t, size_t>& seed_to_ending);
+
+/**
+ * Calculate read distances for each of the zip tree's transitions.
+ * Also filters out transitions that can't be used,
+ * e.g. not reachable in the read.
+ * 
+ * Used as a helper by zip_tree_transition_iterator().
+ * 
+ * Transitions are taken as (source anchor, destination anchor, graph distance)
+ * and output as (source anchor, destination anchor, graph distance, read distance).
+ */
+std::vector<std::tuple<size_t, size_t, size_t, size_t>> calculate_transition_read_distances(
+    const std::vector<std::tuple<size_t, size_t, size_t>>& all_transitions,
+    const VectorView<Anchor>& to_chain,
+    size_t max_read_lookback_bases);
+
 /**
  * Fill in the given DP table for the explored chain scores ending with each
  * item. Returns the best observed score overall from that table, with

src/algorithms/chain_items.hpp

@@ -614,6 +614,8 @@ string alignment_to_sam_internal(const Alignment& alignment,
                                  bool paired,
                                  const int32_t tlen_max) {
 
+    
+
     // Determine flags, using orientation, next/prev fragments, and pairing status.
     int32_t flags = determine_flag(alignment, refseq, refpos, refrev, mateseq, matepos, materev, tlen, paired, tlen_max);
 
@@ -691,8 +693,12 @@ int32_t determine_flag(const Alignment& alignment,
     // Determine flags, using orientation, next/prev fragments, and pairing status.
     int32_t flags = sam_flag(alignment, refrev, paired);
 
-    // We've observed some reads with the unmapped flag set and also a CIGAR string set, which shouldn't happen.
-    // We will check for this. The CIGAR string will only be set in the output if the alignment has a path.
+    // We've observed some reads with the unmapped flag set and also a CIGAR
+    // string set, which is allowed by the SAM spec but which we are not
+    // supposed to generate.
+    //
+    // We will check for this. The CIGAR string will only be set in the output
+    // if the alignment has a path.
     assert((bool)(flags & BAM_FUNMAP) != (alignment.has_path() && alignment.path().mapping_size()));
 
     if (!((bool)(flags & BAM_FUNMAP)) && paired && !refseq.empty() && refseq == mateseq) {
@@ -725,6 +731,10 @@ int32_t determine_flag(const Alignment& alignment,
         flags |= BAM_FMREVERSE;
     }
 
+    if (is_supplementary(alignment)) {
+        flags |= BAM_FSUPPLEMENTARY;
+    }
+    
     return flags;
 }
 
@@ -805,7 +815,7 @@ bam1_t* alignment_to_bam_internal(bam_hdr_t* header,
                                   const bool refrev,
                                   const vector<pair<int, char>>& cigar,
                                   const string& mateseq,
-                                  const int32_t matepos,
+                                  int32_t matepos,
                                   bool materev,
                                   const int32_t tlen,
                                   bool paired,
@@ -1362,6 +1372,131 @@ vector<pair<int, char>> cigar_against_path(const Alignment& alignment, bool on_r
     return cigar;
 }
 
+vector<pair<int, char>> spliced_cigar_against_path(const Alignment& aln, const PathPositionHandleGraph& graph, const string& path_name, 
+                                                   int64_t pos, bool rev, int64_t min_splice_length) {
+    // the return value
+    vector<pair<int, char>> cigar;
+    
+    if (aln.has_path() && aln.path().mapping_size() > 0) {
+        // the read is aligned to the path
+        
+        path_handle_t path_handle = graph.get_path_handle(path_name);
+        step_handle_t step = graph.get_step_at_position(path_handle, pos);
+        
+        // to indicate whether we've found the edit that corresponds to the BAM position
+        bool found_pos = false;
+        
+        const Path& path = aln.path();
+        for (size_t i = 0; i < path.mapping_size(); ++i) {
+            
+            // we traverse backwards on a reverse strand mapping
+            const Mapping& mapping = path.mapping(rev ? path.mapping_size() - 1 - i : i);
+            
+            for (size_t j = 0; j < mapping.edit_size(); ++j) {
+                                
+                // we traverse backwards on a reverse strand mapping
+                const Edit& edit = mapping.edit(rev ? mapping.edit_size() - 1 - j : j);
+                
+                if (!found_pos) {
+                    // we may still be searching through an initial softclip to find
+                    // the edit that corresponds to the BAM position
+                    if (edit.to_length() > 0 && edit.from_length() == 0) {
+                        append_cigar_operation(edit.to_length(), 'S', cigar);
+                        // skip the main block where we assign cigar operations
+                        continue;
+                    }
+                    else {
+                        found_pos = true;
+                    }
+                }
+                
+                // identify the cigar operation
+                char cigar_code;
+                int length;
+                if (edit.from_length() == edit.to_length()) {
+                    cigar_code = 'M';
+                    length = edit.from_length();
+                }
+                else if (edit.from_length() > 0 && edit.to_length() == 0) {
+                    cigar_code = 'D';
+                    length = edit.from_length();
+                }
+                else if (edit.to_length() > 0 && edit.from_length() == 0) {
+                    cigar_code = 'I';
+                    length = edit.to_length();
+                }
+                else {
+                    throw std::runtime_error("Spliced CIGAR construction can only convert simple edits");
+                }
+                
+                append_cigar_operation(length, cigar_code, cigar);
+            } // close loop over edits
+            
+            if (found_pos && i + 1 < path.mapping_size()) {
+                // we're anchored on the path by the annotated position, and we're transitioning between
+                // two mappings, so we should check for a deletion/splice edge
+                
+                step_handle_t next_step = graph.get_next_step(step);
+                
+                handle_t next_handle = graph.get_handle_of_step(next_step);
+                const Position& next_pos = path.mapping(rev ? path.mapping_size() - 2 - i : i + 1).position();
+                if (graph.get_id(next_handle) != next_pos.node_id()
+                    || (graph.get_is_reverse(next_handle) != next_pos.is_reverse()) != rev) {
+                    
+                    // the next mapping in the alignment is not the next mapping on the path, so we must have
+                    // taken a deletion
+                    
+                    // find the closest step that is further along the path than the current one
+                    // and matches the next position (usually there will only be one)
+                    size_t curr_offset = graph.get_position_of_step(step);
+                    size_t nearest_offset = numeric_limits<size_t>::max();
+                    graph.for_each_step_on_handle(graph.get_handle(next_pos.node_id()),
+                                                  [&](const step_handle_t& candidate) {
+                        
+                        if (graph.get_path_handle_of_step(candidate) == path_handle) {
+                            size_t candidate_offset = graph.get_position_of_step(candidate);
+                            if (candidate_offset < nearest_offset && candidate_offset > curr_offset) {
+                                nearest_offset = candidate_offset;
+                                next_step = candidate;
+                            }
+                        }
+                    });
+                    
+                    if (nearest_offset == numeric_limits<size_t>::max()) {
+                        throw std::runtime_error("Spliced BAM conversion could not find path steps that match alignment");
+                    }
+                    
+                    // the gap between the current step and the next one along the path
+                    size_t deletion_length = (nearest_offset - curr_offset -
+                                              graph.get_length(graph.get_handle_of_step(step)));
+                    
+                    // add to the cigar
+                    if (deletion_length >= min_splice_length) {
+                        // long enough to be a splice
+                        append_cigar_operation(deletion_length, 'N', cigar);
+                    }
+                    else if (deletion_length) {
+                        // create or extend a deletion
+                        append_cigar_operation(deletion_length, 'D', cigar);
+                    }
+                }
+                
+                // iterate along the path
+                step = next_step;
+            }
+        } // close loop over mappings
+        
+        if (cigar.back().second == 'I') {
+            // the final insertion is actually a softclip
+            cigar.back().second = 'S';
+        }
+    }
+    
+    simplify_cigar(cigar);
+    
+    return cigar;
+}
+
 void simplify_cigar(vector<pair<int, char>>& cigar) {
 
     size_t removed = 0;
@@ -1664,7 +1799,7 @@ int32_t sam_flag(const Alignment& alignment, bool on_reverse_strand, bool paired
         // TODO: catch reads with pair partners when they shouldn't be paired?
     }
 
-    if (!alignment.has_path() || alignment.path().mapping_size() == 0) {
+    if (!is_mapped(alignment)) {
         // unmapped
         flag |= BAM_FUNMAP;
     } 
@@ -1823,18 +1958,28 @@ Alignment bam_to_alignment(const bam1_t *b,
 
     }
     alignment.set_read_paired((b->core.flag & BAM_FPAIRED) != 0);
+
+    // Consult the One True Place in SAM for knowing if a read is aligned or not.
+    bool is_aligned = !(b->core.flag & BAM_FUNMAP);
 
-    if (graph != nullptr && bh != nullptr && b->core.tid >= 0) {
+    if (is_aligned) {
+        // Set the little-used GAM field that exists to represent this.
+        alignment.set_read_mapped(is_aligned);
+
+        // Use the MAPQ
         alignment.set_mapping_quality(b->core.qual);
-        alignment.set_read_mapped(true);
-        // Look for the path handle this is against.
-        auto found = tid_path_handle.find(b->core.tid);
-        if (found == tid_path_handle.end()) {
-            cerr << "[vg::alignment.cpp] error: alignment references path not present in graph: "
-                 << bh->target_name[b->core.tid] << endl;
-            exit(1);
+        
+        if (graph != nullptr && bh != nullptr && b->core.tid >= 0) {
+            // We can actually build the path for this read.
+            // Look for the path handle this is against.
+            auto found = tid_path_handle.find(b->core.tid);
+            if (found == tid_path_handle.end()) {
+                cerr << "[vg::alignment.cpp] error: alignment references path not present in graph: "
+                     << bh->target_name[b->core.tid] << endl;
+                exit(1);
+            }
+            mapping_against_path(alignment, b, found->second, graph, b->core.flag & BAM_FREVERSE);
         }
-        mapping_against_path(alignment, b, found->second, graph, b->core.flag & BAM_FREVERSE);
     }
 
     // TODO: htslib doesn't wrap this flag for some reason.
@@ -1856,7 +2001,8 @@ Alignment bam_to_alignment(const bam1_t *b,
             }
             ++removed;
         }
-        else if (tag_name == "AS") {
+        else if (tag_name == "AS" && is_aligned) {
+            // Set the score, for aligned reads.
             alignment.set_score(parse<int64_t>(tag.substr(5, string::npos)));
             ++removed;
         }
@@ -2228,9 +2374,8 @@ int softclip_trim(Alignment& alignment) {
 }
 
 bool is_perfect(const Alignment& alignment) {
-    // Non-aligned paths can't be perfect (code from subcommand/stats_main.cpp)
-    bool has_alignment = alignment.score() > 0 || alignment.path().mapping_size() > 0;
-    if (!has_alignment) return false;
+    // Non-aligned paths can't be perfect
+    if (!is_mapped(alignment)) return false;
 
     // Check that the path is perfect
     for (size_t i = 0; i < alignment.path().mapping_size(); ++i) {
@@ -2242,6 +2387,79 @@ bool is_perfect(const Alignment& alignment) {
     return true;
 }
 
+bool is_supplementary(const Alignment& alignment) {
+    if (!has_annotation(alignment, "supplementary")) {
+        return false;
+    }
+    else {
+        return get_annotation<bool>(alignment, "supplementary");
+    }
+}
+
+pair<int64_t, int64_t> aligned_interval(const Alignment& aln) {
+    return pair<int64_t, int64_t>(softclip_start(aln), aln.sequence().size() - softclip_end(aln));
+}
+
+string mate_info(const string& path, int32_t pos, bool rev_strand, bool is_read1) {
+    subrange_t subrange;
+    string path_name = Paths::strip_subrange(path, &subrange);
+    if (subrange != PathMetadata::NO_SUBRANGE) {
+        pos += subrange.first;
+    }
+    string info;
+    info.append(path_name);
+    info.push_back('\t');
+    info.append(to_string(pos));
+    info.push_back('\t');
+    info.push_back(rev_strand ? '1' : '0');
+    info.push_back(is_read1 ? '1' : '0');
+    return info;
+}
+
+tuple<string, int32_t, bool, bool> parse_mate_info(const string& info) {
+    tuple<string, int64_t, bool, bool> parsed;
+    size_t i = 0;
+    while (info[i] != '\t') {
+        ++i;
+    }
+    get<0>(parsed) = info.substr(0, i);
+    ++i;
+    size_t j = i;
+    while (info[j] != '\t') {
+        ++j;
+    }
+    get<1>(parsed) = stoi(info.substr(i, j - i));
+    get<2>(parsed) = info[j + 1] == '1';
+    get<3>(parsed) = info[j + 2] == '1';
+    return parsed;
+}
+
+bool is_mapped(const Alignment& alignment) {
+    if (alignment.read_mapped()) {
+        // The flag for being mapped is set.
+        return true;
+    }
+
+    // Otherwise it's not set, but we don't consistently set it in vg. (SAM
+    // uses an *un*mapped flag so you don't forget to set it on mapped reads
+    // ever.) So second-guess the flag and up on round-tripping SAM with
+    // alignment fields set for unmapped reads.
+    
+    if (alignment.path().mapping_size() > 0) {
+        // If an alignment path is filled in, it is mapped.
+        return true;
+    }
+
+    if (alignment.score() > 0) {
+        // If there's no alignment actually there but anonzero score, we
+        // represent a mapped read, we're just not sure where to.
+        return true;
+    }
+
+    // If there's no evidence of being mapped, we're unmapped.
+    return false;
+}
+
 int query_overlap(const Alignment& aln1, const Alignment& aln2) {
     if (!alignment_to_length(aln1) || !alignment_to_length(aln2)
         || !aln1.path().mapping_size() || !aln2.path().mapping_size()