C++ Mapping::position方法代码示例

void flip_nodes(Alignment& a, set<int64_t> ids, const std::function<size_t(int64_t)>& node_length) {
    Path* path = a.mutable_path();
    for(size_t i = 0; i < path->mapping_size(); i++) {
        // Grab each mapping (includes its position)
        Mapping* mapping = path->mutable_mapping(i);
        if(ids.count(mapping->position().node_id())) {
            // We need to flip this mapping
            *mapping = reverse_mapping(*mapping, node_length);
        } 
    }
}

/// Find the region of the Mapping's node used by the Mapping, in forward strand space, as start to past_end.
static pair<size_t, size_t> mapping_to_range(const xg::XG* xg_index, const Mapping& mapping) {
    // How much of the node does it cover?
    auto mapping_length = mapping_from_length(mapping);
    
    // Work out where the start and past-end positions on the node's forward strand are.
    pair<size_t, size_t> node_range;
    if (mapping.position().is_reverse()) {
        // On the reverse strand we need the node length
        // TODO: getting it can be slow
        auto node_length = xg_index->node_length(mapping.position().node_id());
        
        node_range.first = node_length - mapping.position().offset() - mapping_length;
        node_range.second = node_length - mapping.position().offset();
    } else {
        // On the forward strand this is easy
        node_range.first = mapping.position().offset();
        node_range.second = node_range.first + mapping_length;
    }
    
    return node_range;
}

vector<int> Vectorizer::alignment_to_a_hot(Alignment a){
    int64_t entity_size = my_xg->node_count + my_xg->edge_count;
    vector<int> ret(entity_size, 0);
    Path path = a.path();
    for (int i = 0; i < path.mapping_size(); i++){
        Mapping mapping = path.mapping(i);
        if(! mapping.has_position()){
            continue;
        }
        Position pos = mapping.position();
        int64_t node_id = pos.node_id();
        int64_t key = my_xg->node_rank_as_entity(node_id);
        // Okay, solved the previous out of range errors:
        // We have to use an entity-space that is |nodes + edges + 1|
        // as nodes are indexed from 1, not from 0.
        //TODO: this means we may one day have to do the same bump up
        // by one for edges, as I assume they are also indexed starting at 1.
        //cerr << key << " - " << entity_size << endl;

        //Find edge by current / previous node ID
        // we can check the orientation, though it shouldn't **really** matter
        // whether we catch them in the forward or reverse direction.
        if (i > 0){
            Mapping prev_mapping = path.mapping(i - 1);
            Position prev_pos = prev_mapping.position();
            int64_t prev_node_id = prev_pos.node_id();
            if (my_xg->has_edge(prev_node_id, false, node_id, false)){
                int64_t edge_key = my_xg->edge_rank_as_entity(prev_node_id, false, node_id, false);
                vector<size_t> edge_paths = my_xg->paths_of_entity(edge_key);
                if (edge_paths.size() > 0){
                    ret[edge_key - 1] = 1;
                }
                else{
                    ret[edge_key - 1] = 2;
                }
            }
        }
        //Check if the node of interest is on a path
        vector<size_t> node_paths = my_xg->paths_of_node(node_id);
        if (node_paths.size() > 0){
            ret[key - 1] = 2;
        }
        else{
            ret[key - 1] = 1;
        }

    }

    return ret;

}

vector<double> Vectorizer::alignment_to_identity_hot(Alignment a){
    int64_t entity_size = my_xg->node_count + my_xg->edge_count;
    vector<double> ret(entity_size, 0.0);

    Path path = a.path();
    for (int i = 0; i < path.mapping_size(); i ++){
        Mapping mapping = path.mapping(i);
        if(! mapping.has_position()){
            continue;
        }
        Position pos = mapping.position();
        int64_t node_id = pos.node_id();
        int64_t key = my_xg->node_rank_as_entity(node_id);

        //Calculate % identity by walking the edits and counting matches.
        double pct_id = 0.0;
        double match_len = 0.0;
        double total_len = 0.0;

        for (int j = 0; j < mapping.edit_size(); j++){
            Edit e = mapping.edit(j);
            total_len += e.from_length();
            if (e.from_length() == e.to_length() && e.sequence() == ""){
                match_len += (double) e.to_length();
            }
            else if (e.from_length() == e.to_length() && e.sequence() != ""){
                // TODO if we map but don't match exactly, add half the average length to match_length
                //match_len += (double) (0.5 * ((double) e.to_length()));
            }
            else{
                
            }
            
        }
        pct_id = (match_len == 0.0 && total_len == 0.0) ? 0.0 : (match_len / total_len);
        ret[key - 1] = pct_id;

        if (i > 0){
            Mapping prev_mapping = path.mapping(i - 1);
            Position prev_pos = prev_mapping.position();
            int64_t prev_node_id = prev_pos.node_id();
            if (my_xg->has_edge(prev_node_id, false, node_id, false)){
                int64_t edge_key = my_xg->edge_rank_as_entity(prev_node_id, false, node_id, false);
                ret[edge_key - 1] = 1.0;
            }
        }
    }
    return ret;
}

    /**
     * Looks for alignments that change direction over their length.
     * This may happen because of:
     * 1. Mapping artifacts
     * 2. Cycles
     * 3. Highly repetitive regions
     * 4. Inversions (if you're lucky enough)
     *
     * Default behavior: if the Alignment reverses, return an empty Alignment.
     * inverse behavior: if the Alignment reverses, return the Alignment.
     */
    Alignment Filter::reversing_filter(Alignment& aln){

        Path path = aln.path();
        bool prev = false;

        for (int i = 1; i < path.mapping_size(); i++){
            Mapping mapping = path.mapping(i);
            Position pos = mapping.position();
            bool prev = path.mapping(i - 1).position().is_reverse();
            if (prev != pos.is_reverse()){
                return inverse ? aln : Alignment();
            }

        }
        return inverse ? Alignment() : aln;

    }

    /**
     *
     * Looks for alignments that transition from one path to another
     * over their length. This may occur for one of several reasons:
     * 1. The read covers a translocation
     * 2. The read looks a lot like two different (but highly-similar paths)
     * 3. The read is shattered (e.g. as in chromothripsis)
     *
     * Default behavior: if the Alignment is path divergent, return an empty Alignment, else return aln
     * Inverse behavior: if the Alignment is path divergent, return aln, else return an empty Alignment
     */
    Alignment Filter::path_divergence_filter(Alignment& aln){
        Path path = aln.path();
        for (int i = 1; i < path.mapping_size(); i++){
            Mapping mapping = path.mapping(i);
            Position pos = mapping.position();
            id_t current_node = pos.node_id();
            id_t prev_node = path.mapping(i - 1).position().node_id();
            bool paths_match = false;
            vector<size_t> paths_of_prev = my_xg_index->paths_of_node(prev_node);
            for (int i = 0; i < paths_of_prev.size(); i++){
                string p_name = my_xg_index->path_name(paths_of_prev[i]);
                if (my_xg_index->path_contains_node(p_name, current_node)){
                    paths_match = true;
                }
            }
            if (!paths_match){
                return inverse ? aln : Alignment();
            }

        }
        return inverse ? Alignment() : aln;
    }

bit_vector Vectorizer::alignment_to_onehot(Alignment a){
    // Make a vector as large as the | |nodes| + |edges| | space
    // TODO handle edges
    int64_t entity_size = my_xg->node_count + my_xg->edge_count;
    bit_vector ret(entity_size, 0);
    Path path = a.path();
    for (int i = 0; i < path.mapping_size(); i++){
        Mapping mapping = path.mapping(i);
        Position pos = mapping.position();
        int64_t node_id = pos.node_id();
        int64_t key = my_xg->node_rank_as_entity(node_id);
        // Okay, solved the previous out of range errors:
        // We have to use an entity-space that is |nodes + edges + 1|
        // as nodes are indexed from 1, not from 0.
        //TODO: this means we may one day have to do the same bump up
        // by one for edges, as I assume they are also indexed starting at 1.
        //cerr << key << " - " << entity_size << endl;

        //Find edge by current / previous node ID
        // we can check the orientation, though it shouldn't **really** matter
        // whether we catch them in the forward or reverse direction.
        if (i > 0){
            Mapping prev_mapping = path.mapping(i - 1);
            Position prev_pos = prev_mapping.position();
            int64_t prev_node_id = prev_pos.node_id();
            if (my_xg->has_edge(prev_node_id, false, node_id, false)){
                int64_t edge_key = my_xg->edge_rank_as_entity(prev_node_id, false, node_id, false);
                ret[edge_key - 1] = 1;
            }
        }
        //Find entity rank of edge

        ret[key - 1] = 1;
    }

    return ret;
}

void Pileups::compute_from_edit(NodePileup& pileup, int64_t& node_offset,
                                int64_t& read_offset,
                                const Node& node, const Alignment& alignment,
                                const Mapping& mapping, const Edit& edit) {
    string seq = edit.sequence();
    bool is_reverse = mapping.position().is_reverse();
    
    // ***** MATCH *****
    if (edit.from_length() == edit.to_length()) {
        assert (edit.from_length() > 0);
        make_match(seq, edit.from_length(), is_reverse);
        assert(seq.length() == edit.from_length());            
        int64_t delta = 1;
        for (int64_t i = 0; i < edit.from_length(); ++i) {
            BasePileup* base_pileup = get_create_base_pileup(pileup, node_offset);
            // reference_base if empty
            if (base_pileup->num_bases() == 0) {
                base_pileup->set_ref_base(node.sequence()[node_offset]);
            } else {
                assert(base_pileup->ref_base() == node.sequence()[node_offset]);
            }
            // add base to bases field (converting to ,. if match)
            char base = seq[i];
            if (!edit.sequence().empty() &&
                base_equal(seq[i], node.sequence()[node_offset], is_reverse)) {
                base = is_reverse ? ',' : '.';
            }
            *base_pileup->mutable_bases() += base;
            // add quality if there
            if (!alignment.quality().empty()) {
                *base_pileup->mutable_qualities() += alignment.quality()[read_offset];
            }
            // pileup size increases by 1
            base_pileup->set_num_bases(base_pileup->num_bases() + 1);
            // move right along read, and left/right depending on strand on reference
            node_offset += delta;
            ++read_offset;
        }
    }
    // ***** INSERT *****
    else if (edit.from_length() < edit.to_length()) {
        make_insert(seq, is_reverse);
        assert(edit.from_length() == 0);
        // we define insert (like sam) as insertion between current and next
        // position (on forward node coordinates). this means an insertion before
        // offset 0 is invalid! 
        int64_t insert_offset =  is_reverse ? node_offset : node_offset - 1;
        if (insert_offset >= 0) {        
            BasePileup* base_pileup = get_create_base_pileup(pileup, insert_offset);
            // reference_base if empty
            if (base_pileup->num_bases() == 0) {
                base_pileup->set_ref_base(node.sequence()[insert_offset]);
            } else {
                assert(base_pileup->ref_base() == node.sequence()[insert_offset]);
            }
            // add insertion string to bases field
            // todo: should we reverse complement this if mapping is reversed ??? 
            base_pileup->mutable_bases()->append(seq);
            if (!alignment.quality().empty()) {
                *base_pileup->mutable_qualities() += alignment.quality()[read_offset];
            }
            // pileup size increases by 1
            base_pileup->set_num_bases(base_pileup->num_bases() + 1);
        }
        else {
            // todo: need to either forget about these, or extend pileup format.
            // easy solution: change insert to come before position, and just add
            // optional pileup at n+1st base of node.  would like to figure out
            // how samtools does it first...
            /*
            stringstream ss;
            ss << "Warning: pileup does not support insertions before 0th base in node."
               << " Offending edit: " << pb2json(edit) << endl;
#pragma omp critical(cerr)
            cerr << ss.str();
            */
        }
        // move right along read (and stay put on reference)
        read_offset += edit.to_length();
    }
    // ***** DELETE *****
    else {
        assert(edit.to_length() == 0);
        assert(edit.sequence().empty());
        int64_t del_start = !is_reverse ? node_offset :
            node_offset - edit.from_length() + 1;
        seq = node.sequence().substr(del_start, edit.from_length());
        make_delete(seq, is_reverse);
        BasePileup* base_pileup = get_create_base_pileup(pileup, node_offset);
        // reference_base if empty
        if (base_pileup->num_bases() == 0) {
            base_pileup->set_ref_base(node.sequence()[node_offset]);
        } else {
            assert(base_pileup->ref_base() == node.sequence()[node_offset]);
        }
        // add deletion string to bases field
        // todo: should we reverse complement this if mapping is reversed ??? 
        base_pileup->mutable_bases()->append(seq);
        if (!alignment.quality().empty()) {
            *base_pileup->mutable_qualities() += alignment.quality()[read_offset];
//.........这里部分代码省略.........

    void Deconstructor::sb2vcf(string outfile){
        Header h;
        h.set_date();
        h.set_source("VG");
        h.set_reference("");
        h.set_version("VCF4.2");

        cout << h << endl;

        // for each superbubble:
        // Fill out a vcflib Variant
        // Check if it is masked by an input vcf
        // if not, print it to stdout



        map<id_t, vcflib::Variant> node_to_var;
        vcflib::VariantCallFile mask;
        if (!mask_file.empty()){
            //node_to_var = my_vg->get_node_to_variant(mask);
        }
        for (auto s : my_sbs){
            vcflib::Variant var;

            // Make subgraphs out of the superbubble:
            // Operating on a pair<id_t, id_t>, vector<id_t>
            // then enumerate k_paths through the SuperBubbles
            set<Node*> nodes;
            set<Edge*> edges;

            for (int i = 0; i < s.second.size(); i++){
                id_t n_id = s.second[i];
                //cerr << n_id << endl;
                Node* n_node = my_vg->get_node(n_id);
                vector<Edge*> e_end = my_vg->edges_from(n_node);
                nodes.insert(n_node);
                if (i < s.second.size() - 1){
                    edges.insert(e_end.begin(), e_end.end());
                }
            }

            vg::VG t_graph = vg::VG(nodes, edges);

            vector<Path> paths;

            std::function<void(NodeTraversal)> no_op = [](NodeTraversal n){};
            std::function<void(size_t, Path&)> extract_path = [&paths](size_t x_size, Path& path){
                paths.push_back(path);
            };

            t_graph.for_each_kpath(10000, false, 100, no_op, no_op, extract_path);

            std::function<std::vector<Path>(vector<Path>)> uniquify = [](vector<Path> v){
                map<string, Path> unqs;
                vector<Path> ret;
                for (auto x: v){
                    unqs[path_to_string(x)] = x;
                }

                for (auto y : unqs){
                    ret.push_back(y.second);
                }
                return ret;
            };

            paths = uniquify(paths);

            std::function<bool(Path)> all_ref = [&](Path p){
                for (int i = 0; i < p.mapping_size(); i++){
                    Mapping m = p.mapping(i);
                    Position pos = m.position();
                    vg::id_t pos_id = pos.node_id();
                    map<string, set<Mapping*> > path_to_mappings =  my_vg->paths.get_node_mapping(pos_id);

                    if (path_to_mappings.size() <= 0){
                        return false;
                    }
                }
                return true;
            };

            /*
             * This means we now have vectors for the superbubble
             * that have the paths through the nodes within it (including end nodes)
             * however, these paths are repeated several times.
             * We should find a way to prevent them being inserted once for each node.
             *
             * Next on the agenda: use the get_path_dist thing from vg call / vg stats
             * to get the distance to the head node.
             * Might need an XG index for this.
             *
             * Also need a way to deal with GAMs for this i.e. a way to 
             * count the number of times we see something come up in the gam
             */
            int first_len = (my_vg->get_node(1))->sequence().size();
            map<string, set<Mapping*> > p_to_mappings =  my_vg->paths.get_node_mapping(s.first.first);
            for (auto p_name : p_to_mappings){
                var.sequenceName = p_name.first;
            }
            var.position = my_xg->approx_path_distance(var.sequenceName, 1, s.first.first) + (s.first.first == 1 ? 0 : first_len);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        // node, so we can go through them while tampering with them.
        map<string, set<Mapping*> > mappings_by_path = graph.paths.get_node_mapping(graph.get_node(leaf->start().node_id()));
        
        // It's possible a path can enter the site through the end node and
        // never hit the start. So we're going to trim those back before we delete nodes and edges.
        map<string, set<Mapping*> > end_mappings_by_path = graph.paths.get_node_mapping(graph.get_node(leaf->end().node_id()));
        
        if (!drop_hairpin_paths) {
            // We shouldn't drop paths if they hairpin and can't be represented
            // in a simplified bubble. So we instead have to not simplify
            // bubbles that would have that problem.
            bool found_hairpin = false;
            
            for (auto& kv : mappings_by_path) {
                // For each path that hits the start node
                
                if (found_hairpin) {
                    // We only care if there are 1 or more hairpins, not how many
                    break;    
                }
                
                // Unpack the name
                auto& path_name = kv.first;
                
                for (Mapping* start_mapping : kv.second) {
                    // For each visit to the start node
                
                    if (found_hairpin) {
                        // We only care if there are 1 or more hairpins, not how many
                        break;    
                    }
                
                    // Determine what orientation we're going to scan in
                    bool backward = start_mapping->position().is_reverse();
                    
                    // Start at the start node
                    Mapping* here = start_mapping;
                    
                    while (here) {
                        // Until we hit the start/end of the path or the mapping we want
                        if (here->position().node_id() == leaf->end().node_id() &&
                            here->position().is_reverse() == (leaf->end().backward() != backward)) {
                            // We made it out.
                            // Stop scanning!
                            break;
                        }
                        
                        if (here->position().node_id() == leaf->start().node_id() &&
                            here->position().is_reverse() != (leaf->start().backward() != backward)) {
                            // We have encountered the start node with an incorrect orientation.
                            cerr << "warning:[vg simplify] Path " << path_name
                                << " doubles back through start of site "
                                << to_node_traversal(leaf->start(), graph) << " - "
                                << to_node_traversal(leaf->end(), graph) << "; skipping site!" << endl;
                                
                            found_hairpin = true;
                            break;
                        }
                        
                        // Scan left along ther path if we found the site start backwards, and right if we found it forwards.
                        here = backward ? graph.paths.traverse_left(here) : graph.paths.traverse_right(here);
                    }
                }
            }
            
            for (auto& kv : end_mappings_by_path) {

    Alignment Filter::depth_filter(Alignment& aln){
        if (use_avg && window_length != 0){

        }
        else if (use_avg != 0){

        }
        else{

        }

        Path path = aln.path();
        //TODO handle reversing mappings
        vector<int>* qual_window;
        if (window_length > 0){
            qual_window = new vector<int>();
        }

        for (int i = 0; i < path.mapping_size(); i++){
            Mapping mapping = path.mapping(i);
            Position start_pos = mapping.position();
            int64_t start_node = start_pos.node_id();
            int64_t start_offset = start_pos.offset();
            int64_t curr_offset_in_graph = 0;
            int64_t curr_offset_in_alignment = 0;
            stringstream pst;

            pst << start_node << "_" << curr_offset_in_graph;
            string p_hash = pst.str();
            for (int j = 0; j < mapping.edit_size(); j++){
                Edit ee = mapping.edit(j);
                if (ee.from_length() == ee.to_length() && ee.sequence() == ""){
                    if (!filter_matches){
                        continue;
                    }
                }
                stringstream est;
                est <<  ee.from_length() << "_" << ee.to_length() << "_" + ee.sequence();
                string e_hash = est.str();
#pragma omp critical(write)
                pos_to_edit_to_depth[p_hash][e_hash] += 1;
                /**
                 * If an edit fails the filter, either return a new empty alignment
                 * OR
                 * return a new alignment identical to the old one EXCEPT where
                 * the offending edit has been replaced by a match to the reference.
                 */
                if (pos_to_edit_to_depth[p_hash][e_hash] < min_depth){
                    if (!remove_failing_edits){
                        return inverse ? aln : Alignment();
                    }

                    else {
                        Alignment edited_aln = Alignment(aln);
                        edited_aln.mutable_path()->mutable_mapping(i)->mutable_edit(j)->set_sequence("");
                        edited_aln.mutable_path()->mutable_mapping(i)->mutable_edit(j)->set_from_length(ee.from_length());
                        edited_aln.mutable_path()->mutable_mapping(i)->mutable_edit(j)->set_to_length(ee.from_length());
                        return edited_aln;
                    }
                }
            }
            return inverse ? Alignment() : aln;
        }


    }