<|start_filename|>filter-paf-deletions.cpp<|end_filename|>
// Filter big deletions from PAF.  These can sometimes arise from minigraph split alignments.  They are rare, but can really
// mess up topology of graph if even one gets into cactus.

// 1) Estimate anchors along reference path for every node in input graph
// 2) Scan every query in PAF in order, and look at target blocks
// 3) Use table from 1) in order to estimate the distances between target blocks
// 4) If two conesecutive target blocks span too big of a distance, delete the smaller block

//#define debug

#include <cstdlib>
#include <iostream>
#include <cassert>
#include <fstream>
#include <deque>
#include <unordered_map>
#include <unistd.h>
#include <getopt.h>
#include <omp.h>

#include "bdsg/packed_graph.hpp"
#include "bdsg/hash_graph.hpp"
#include "bdsg/odgi.hpp"

#include "IntervalTree.h"

using namespace std;
using namespace handlegraph;
using namespace bdsg;

struct Anchor {
    path_handle_t path_handle;
    int64_t max_offset;
    int64_t min_offset;
};

static unique_ptr<MutablePathMutableHandleGraph> load_graph(istream& graph_stream);
static pair<unordered_map<string, nid_t>, unordered_map<nid_t, string>> load_trans(const string& trans_path);
static vector<string> &split_delims(const string &s, const string& delims, vector<string> &elems);
static unordered_map<nid_t, Anchor> index_graph(const PathHandleGraph* graph,
                                                const string& ref_prefix);
static IntervalTree<int64_t, int64_t> index_deletions(const PathHandleGraph* graph, const unordered_map<nid_t, Anchor>& index);
static int64_t filter_paf(const PathHandleGraph* graph, ifstream& paf_file, const unordered_map<nid_t, Anchor>& ref_index,
                          const IntervalTree<int64_t, int64_t>& ref_deletions,
                          const unordered_map<string, nid_t>& mg_to_vg, int64_t max_deletion,
                          double overlap_threshold, vector<bool>& filtered_lines, bool progress, bool verbose);

void help(char** argv) {
    cerr << "usage: " << argv[0] << " [options] <graph.vg> <trans> <aln.paf> <threshold>\n" << endl
         << "Use distances from graph to filter out implied deletions from PAF (cigars not considered, only blocks)" << endl
         << "  <graph.vg> : minigraph as obtained from vg convert -g graph.gfa" << endl
         << "  <trans> : node translation from vg convert -g -T" << endl
         << "  <aln.paf> : paf alignment from cactus-graphmap" << endl
         << "  <threshold> : only remove deletions greater than this" << endl
         << endl
         << "options: " << endl
         << "    -r, --ref-prefix STR      Only consider paths whose names start with STR" << endl
         << "    -p, --progress            Print progress" << endl
         << "    -v, --verbose             Print deletions" << endl
         << "    -t, --threads N           number of threads to use (used only for computing snarls) [default: all available]" << endl
       << endl;
}    

int main(int argc, char** argv) {

    string ref_prefix;
    bool progress = false;
    bool verbose = false;
    // only filter deletions that don't overlap an existing deletion by at least this much
    // (doesn't seem to a factor -- most big deletions not in minigraph)
    double overlap_threshold = 0.5;
    int c;
    optind = 1; 
    while (true) {

        static const struct option long_options[] = {
            {"ref-prefix", required_argument, 0, 'r'},
            {"help", no_argument, 0, 'h'},
            {"progress", no_argument, 0, 'p'},
            {"verbose", no_argument, 0, 'v'},
            {"threads", required_argument, 0, 't'},
            {0, 0, 0, 0}
        };

        int option_index = 0;

        c = getopt_long (argc, argv, "hpvt:",
                         long_options, &option_index);

        // Detect the end of the options.
        if (c == -1)
            break;

        switch (c)
        {
        case 'r':
            ref_prefix = optarg;
            break;
        case 'v':
            verbose = true;
            break;
        case 'p':
            progress = true;
            break;
        case 't':
        {
            int num_threads = stoi(optarg);
            if (num_threads <= 0) {
                cerr << "[filter-paf-deletions] error: Thread count (-t) set to " << num_threads << ", must set to a positive integer." << endl;
                exit(1);
            }
            omp_set_num_threads(num_threads);
            break;
        }                        
        case 'h':
        case '?':
            /* getopt_long already printed an error message. */
            help(argv);
            exit(1);
            break;
        default:
            abort ();
        }
    }

    if (argc <= 4) {
        cerr << "[filter-paf-deletions] error: too few arguments\n" << endl;
        help(argv);
        return 1;
    }

    // Parse the positional argument
    if (optind >= argc) {
        cerr << "[filter-paf-deletions] error: too few arguments\n" << endl;
        help(argv);
        return 1;
    }

    if (optind != argc - 4) {
        cerr << "[filter-paf-deletions] error: too many arguments\n" << endl;
        help(argv);
        return 1;
    }
    
    string graph_path = argv[optind++];
    string trans_path = argv[optind++];
    string paf_path = argv[optind++];
    int64_t max_deletion = stol(argv[optind++]);

    // load the graph
    ifstream graph_stream(graph_path);
    if (!graph_stream) {
        cerr << "[filter-paf-deletions] error: Unable to open input graph " << graph_path << endl;
        return 1;
    }    
    unique_ptr<PathHandleGraph> graph = load_graph(graph_stream);
    graph_stream.close();
    if (progress) {
        cerr << "[filter-paf-deletions]: Loaded graph" << endl;
    }

    // load the minigraph <-> vg id translation table (because our PAF is expressed in terms of the minigraph
    // ids but we lose them when converting to vg.)
    unordered_map<string, nid_t> mg_to_vg;
    unordered_map<nid_t, string> vg_to_mg;
    std::tie(mg_to_vg, vg_to_mg) = load_trans(trans_path);

    if (progress) {
        cerr << "[filter-paf-deletions]: Loaded " << mg_to_vg.size() << " translations." << endl;
    }

    // open the paf
    ifstream paf_file(paf_path);
    if (!paf_file) {
        cerr << "[filter-paf-deletions] error: Unable to open PAF" << endl;
        return 1;
    }

    // index the minigraph
    // this maps each node in the graph to a (maximal) reference interval
    unordered_map<nid_t, Anchor> ref_index = index_graph(graph.get(), ref_prefix);
    if (progress) {
        cerr << "[filter-paf-deletions]: Created reference path index" << endl;
    }

    IntervalTree<int64_t, int64_t> ref_deletions = index_deletions(graph.get(), ref_index);
    if (progress) {
        cerr << "[filter-paf-deletions]: Created reference deletion index" << endl;
    }

#ifdef debug
    for (auto fam : ref_index) {
        cerr << fam.first << " -> " << graph->get_path_name(fam.second.path_handle) << " " << fam.second.min_offset << " " << fam.second.max_offset << endl;
    }
#endif

    // we have everything needed to filter the paf
    vector<bool> filtered_lines;
    int64_t filtered_line_total = 0;
    int64_t filtered_line_it = 0;
    int64_t iteration = 0;
    do {
        filtered_line_it = filter_paf(graph.get(), paf_file, ref_index, ref_deletions, mg_to_vg, max_deletion,
                                      overlap_threshold, filtered_lines, progress, verbose);
        filtered_line_total += filtered_line_it;
        if (progress) {
            cerr << "[filter-paf-deletions]: Iteration " << iteration << ": Found " << filtered_line_it << " lines to filter" << endl;
        }
        ++iteration;
    } while (filtered_line_it > 0);

    if (progress) {
        cerr << "[filter-paf-deletions]: Filtering out " << filtered_line_total << " paf lines" << endl;
    }

    // output the unfiltered lines
    paf_file.clear();
    paf_file.seekg(0, ios::beg) ;
    string buffer;
    for (int64_t line_no = 0; line_no < filtered_lines.size(); ++line_no) {
        const auto& ret = getline(paf_file, buffer);
        assert(ret);
        if (filtered_lines[line_no] == false) {
            cout << buffer << "\n";
        }
    }
    cout << flush;

    return 0;
}

static string strip_prefix(const string& name) {
    if (name.compare(0, 3, "id=") == 0) {
        size_t p = name.find('|', 3);
        assert(p != string::npos);
        return name.substr(p + 1);
    }
    return name;
}

int64_t filter_paf(const PathHandleGraph* graph, ifstream& paf_file, const unordered_map<nid_t, Anchor>& ref_index,
                   const IntervalTree<int64_t, int64_t>& ref_deletions,
                   const unordered_map<string, nid_t>& mg_to_vg, int64_t max_deletion,
                   double overlap_threshold, vector<bool>& filtered_lines, bool progress, bool verbose) {

    paf_file.clear();
    paf_file.seekg(0, ios::beg) ;
    string buffer;

    unordered_set<string> query_set;
    string prev_query_name;
    int64_t prev_query_start;
    int64_t prev_query_end;
    nid_t prev_target_id = 0;
    int64_t prev_ref_start;
    int64_t prev_ref_end;

    // store matches in current block (run of mappings from same query that don't span a max_deletion)
    vector<int64_t> line_to_block; // map paf line number to a block (this vector contains offsets to below vector)
    vector<int64_t> blocks;
    unordered_set<int64_t> cut_blocks; // offsets in above that mark block before max_deletion

    size_t line_no = 0;
    size_t unfiltered_line_no = 0;
    while (getline(paf_file, buffer)) {
        if (unfiltered_line_no >= filtered_lines.size()) {
            // first pass
            filtered_lines.push_back(false);
        } else if (filtered_lines[unfiltered_line_no]) {
            ++unfiltered_line_no;
            continue;
        }
        vector<string> toks;
        split_delims(buffer, "\t\n", toks);

        // pretty silly not to have a paf parser in this repo
        const string& query_name = toks[0];
        int64_t query_start = stol(toks[2]);
        int64_t query_end = stol(toks[3]);
        string target_name = strip_prefix(toks[5]);
        int64_t target_length = stol(toks[6]);
        int64_t target_start = stol(toks[7]);
        int64_t target_end = stol(toks[8]);
        int64_t matches = stol(toks[9]);

        // map the target interval onto the reference path
        if (!mg_to_vg.count(target_name)) {
            cerr << "[filter-paf-deletions] error: target name from PAF not found in translation: " << target_name << endl;
            exit(1);
        }
        nid_t target_id = mg_to_vg.at(target_name);
        const Anchor& anchor = ref_index.at(target_id);
        int64_t ref_start = anchor.min_offset + target_start;
        int64_t ref_end = anchor.max_offset - (target_length - target_end);

#ifdef debug
        cerr << "anchor " << target_id << " = " << anchor.min_offset << " , " << anchor.max_offset << endl;
        cerr << "ref start = " << anchor.min_offset  << " + " << target_start << "   ref_end = " << anchor.max_offset << " - ("
             << target_length << " - " << target_end << ");" << endl;
#endif
        bool continues_block = query_name == prev_query_name;
        if (continues_block) {

            if (prev_query_start > query_start) {
                cerr << "Input paf not in order.  Sort with \"sort -k 1,1 -k 3,3n\" first.  Offending line:\n" << buffer << endl;
                exit(1);
            }
            int64_t query_delta = abs(query_start - prev_query_end);
            int64_t ref_delta = abs(ref_start - prev_ref_end);
            int64_t delta = abs(ref_delta - query_delta);
#ifdef debug
            cerr << "query delta " << query_delta << " ref delta = " << ref_delta << " delta = " << delta << endl;
#endif

            int64_t ref_overlap_size = 0;
            if (delta > 1) {
                vector<Interval<int64_t, int64_t>> overlaps = ref_deletions.findOverlapping(min(prev_query_end, query_start), max(prev_query_end, query_start));
                for (const auto& overlap : overlaps) {
                    int64_t intersection_start = max(prev_query_end, overlap.start);
                    int64_t intersection_stop = min(query_start, overlap.stop);
                    ref_overlap_size = max(ref_overlap_size, intersection_stop - intersection_start);
                }
            }

            if (delta > max_deletion && ((double)ref_overlap_size / delta < overlap_threshold)) {
                if (verbose) {
                    cerr << "[filter-paf-deletions]: detected deletion of size " << delta << " with overlap " << ref_overlap_size <<
                        " on following paf line:"
                         << " \n  " << buffer << endl;
                }
                continues_block = false;
                cut_blocks.insert(blocks.size() - 1);
            }
        }

        // update the current block may adding the matches (or assigng a new block if necessary)
        if (!continues_block) {
            blocks.push_back(0);
        }
        assert(!blocks.empty());
        blocks.back() += matches;
        line_to_block.push_back(blocks.size() - 1);

        prev_query_name = query_name;
        prev_query_start = query_start;
        prev_query_end = query_end;
        prev_target_id = target_id;
        prev_ref_start = ref_start;
        prev_ref_end = ref_end;
        
        ++line_no;
        ++unfiltered_line_no;
    }

    // second pass to do filtering now that we have block sizes
    paf_file.clear();
    paf_file.seekg(0, ios::beg) ;

    size_t filtered_line_count = 0;
        
    line_no = 0;
    unfiltered_line_no = 0;
    while (getline(paf_file, buffer)) {
        if (filtered_lines[unfiltered_line_no]) {
            ++unfiltered_line_no;
            continue;
        }

        bool filter_out = false;        
        int64_t cur_block = line_to_block[line_no];
        if (cur_block > 0) {
            int64_t prev_block = cur_block - 1;
            if (cut_blocks.count(prev_block) && blocks[prev_block] >= blocks[cur_block]) {
                filter_out = true;
            }
        }
        if (!filter_out && cur_block < blocks.size() - 1) {
            int64_t next_block = cur_block + 1;
            if (cut_blocks.count(cur_block) && blocks[next_block] > blocks[cur_block]) {
                filter_out =true;
            }
        }

        if (filter_out) {
            filtered_lines[unfiltered_line_no] = true;
            ++filtered_line_count;
        }
        ++unfiltered_line_no;
        ++line_no;
    }

    return filtered_line_count;
}

unordered_map<nid_t, Anchor> index_graph(const PathHandleGraph* graph,
                                         const string& ref_prefix) {

    // start by making a path position index
    // minigraph assumption: no more than one path per handle!
    unordered_map<handle_t, int64_t> position_index;
    graph->for_each_path_handle([&](path_handle_t path_handle) {            
            if (graph->get_path_name(path_handle).compare(0, ref_prefix.length(), ref_prefix) == 0) {
                size_t offset = 0;
                graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
                        handle_t handle = graph->get_handle_of_step(step_handle);
                        size_t len = graph->get_length(handle);
                        assert(len > 0);
                        assert(!position_index.count(handle));
                        position_index[handle] = offset;
                        position_index[graph->flip(handle)] = offset + len - 1;
                        offset += len;
                    });
            }
        });

    if (position_index.empty()) {
        cerr << "[filter-paf-deletions] error: no reference path found" << endl;
        exit(0);
    }
    
    vector<unordered_map<nid_t, Anchor>> thread_results(get_thread_count());
    
    // really slow brute-force relies on minigraph not having too many nodes
    graph->for_each_handle([&](handle_t handle) {
            unordered_map<nid_t, Anchor>& ref_index = thread_results[omp_get_thread_num()];

            // find all reference nodes that are connected via BFS
            unordered_set<handle_t> context;
            unordered_set<handle_t> ref_handles;
            vector<handle_t> cur_handles = {handle};
            while (!cur_handles.empty()) {
                vector<handle_t> next_handles;
                for (auto& h : cur_handles) {
                    if (!context.count(h)) {
                        context.insert(h);
                        if (position_index.count(h)) {
                            // dead-end on reference
                            ref_handles.insert(h);
                        } else {
                            graph->follow_edges(h, false, [&](handle_t n) {
                                    next_handles.push_back(n);
                                });
                            graph->follow_edges(h, true, [&](handle_t p) {
                                    next_handles.push_back(p);
                                });
                        }
                    }
                }
                cur_handles = std::move(next_handles);
            }

            // update the index with reference offsets
            unordered_set<path_handle_t> ref_path_set;
            int64_t min_ref_offset = numeric_limits<int64_t>::max();
            int64_t max_ref_offset = -1;
            for (handle_t ref_handle : ref_handles) {
                vector<step_handle_t> steps = graph->steps_of_handle(ref_handle);
                assert(steps.size() == 1);
                path_handle_t ref_path_handle = graph->get_path_handle_of_step(steps.back());
                ref_path_set.insert(ref_path_handle);
                // assumption: only one reference path in component
                // (fair for minigraph, but may need to do better than prefix for path selection)
                assert(ref_path_set.size() == 1);
                int64_t ref_offset = position_index.at(ref_handle);
                int64_t ref_offset_rev = position_index.at(graph->flip(ref_handle));
                min_ref_offset = std::min(min_ref_offset, min(ref_offset, ref_offset_rev));
                max_ref_offset = std::max(max_ref_offset, max(ref_offset, ref_offset_rev));
                assert(max_ref_offset >= min_ref_offset);
            }
            if (!ref_path_set.empty()) {
                assert(ref_path_set.size() == 1);
                Anchor& anchor = ref_index[graph->get_id(handle)];
                anchor.path_handle = *ref_path_set.begin();
                anchor.min_offset = min_ref_offset;
                anchor.max_offset = max_ref_offset;
                assert(anchor.max_offset >= anchor.min_offset);
            }
        }, true);

    // merge up the indexes
    for (size_t i = 1; i < thread_results.size(); ++i) {
        for (const auto& id_anchor : thread_results[i]) {
            thread_results[0][id_anchor.first] = id_anchor.second;
        }
        thread_results[i].clear();
    }

    return thread_results[0];
}
    
IntervalTree<int64_t, int64_t> index_deletions(const PathHandleGraph* graph, const unordered_map<nid_t, Anchor>& index) {

    vector<vector<Interval<int64_t, int64_t>>> thread_deletions(get_thread_count());

    // get approximate deletion intervals using the index
    graph->for_each_edge([&](edge_t edge) {
            const Anchor& a1 = index.at(graph->get_id(edge.first));
            const Anchor& a2 = index.at(graph->get_id(edge.second));
            Interval<int64_t, int64_t> interval(0, 0, 0);
            if (a1.min_offset < a2.min_offset) {
                interval.start = a1.max_offset;
                interval.stop = a2.min_offset;
            } else {
                interval.start = a2.max_offset;
                interval.stop = a1.min_offset;                
            }
            interval.value = interval.stop - interval.start;
            if (interval.value > 1) {
                thread_deletions[omp_get_thread_num()].push_back(interval);
            }            
        }, true);

    for (size_t i = 1; i < thread_deletions.size(); ++i) {
        for (const auto& interval : thread_deletions[i]) {
            thread_deletions[0].push_back(interval);
        }
        thread_deletions[i].clear();
    }
    IntervalTree<int64_t, int64_t> tree(thread_deletions[0]);
    return tree;
}


pair<unordered_map<string, nid_t>, unordered_map<nid_t, string>> load_trans(const string& trans_path) {
    ifstream trans_file(trans_path);
    if (!trans_file) {
        cerr << "[filter-paf-deletions] error: Unable to load trans file" << endl;
        exit(1);
    }

    unordered_map<string, nid_t> mg_to_vg;
    unordered_map<nid_t, string> vg_to_mg;

    string buffer;
    while (getline(trans_file, buffer)) {
        vector<string> toks;
        split_delims(buffer, "\t\n", toks);
        assert(toks.size() == 3 && toks[0] == "T");
        mg_to_vg[toks[1]] = stol(toks[2]);
        vg_to_mg[stol(toks[2])] = toks[1];
    }

    return make_pair(mg_to_vg, vg_to_mg);    
}


vector<string> &split_delims(const string &s, const string& delims, vector<string> &elems) {
    size_t start = string::npos;
    for (size_t i = 0; i < s.size(); ++i) {
        if (delims.find(s[i]) != string::npos) {
            if (start != string::npos && i > start) {
                elems.push_back(s.substr(start, i - start));
            }
            start = string::npos;
        } else if (start == string::npos) {
            start = i;
        }
    }
    if (start != string::npos && start < s.size()) {
        elems.push_back(s.substr(start, s.size() - start));
    }
    return elems;
}

unique_ptr<MutablePathMutableHandleGraph> load_graph(istream& graph_stream) {

    char magic_bytes[4];
    graph_stream.read(magic_bytes, 4);
    uint32_t magic_number = ntohl(*((uint32_t*) magic_bytes));
    graph_stream.clear();
    graph_stream.seekg(0, ios::beg);

    MutablePathMutableHandleGraph* graph;
    if (magic_number == PackedGraph().get_magic_number()) {
        graph = new PackedGraph();
    } else if (magic_number == HashGraph().get_magic_number()) {
        graph = new HashGraph();
    } else if (magic_number == ODGI().get_magic_number()) {
        graph = new ODGI();
    } else {
        cerr << "Unable to parse input graph with magic number " << magic_number << endl;
        exit(1);
    }
    dynamic_cast<SerializableHandleGraph*>(graph)->deserialize(graph_stream);

    return unique_ptr<MutablePathMutableHandleGraph>(graph);
}



<|start_filename|>clip-vg.cpp<|end_filename|>
// Chop regions (from BED File) out of paths in vg graphs, creating subpath names and cutting out nodes or parts of nodes
// Assumes that:
// - regions don't overlap (error otherwise)
// - regions are only ever part of at most one path each (assert false otherwise)

//#define debug

#include <cstdlib>
#include <iostream>
#include <cassert>
#include <fstream>
#include <deque>
#include <unordered_map>
#include <unistd.h>
#include <getopt.h>

#include "bdsg/packed_graph.hpp"
#include "bdsg/hash_graph.hpp"
#include "bdsg/odgi.hpp"

using namespace std;
using namespace handlegraph;
using namespace bdsg;

void help(char** argv) {
  cerr << "usage: " << argv[0] << " [options] <graph>" << endl
       << "Chop out path intervals from a vg graph" << endl
       << endl
       << "options: " << endl
       << "    -b, --bed FILE            Intervals to clip in BED format" << endl
       << "    -m, --min-length N        Only clip paths of length < N" << endl
       << "    -u, --max-unaligned N     Clip out unaligned regions of length > N" << endl
       << "    -e, --ref-prefix STR      Forwardize (but don't clip) paths whose name begins with STR" << endl
       << "    -f, --force-clip          Don't abort with error if clipped node overlapped by multiple paths" << endl
       << "    -r, --name-replace S1>S2  Replace (first occurrence of) S1 with S2 in all path names" << endl
       << "    -n, --no-orphan-filter    Don't filter out new subpaths that don't align to anything" << endl
       << "    -d, --drop-path PREFIX    Remove all paths with given PREFIX, and all nodes that are on no other paths (done after other filters)" << endl
       << "    -p, --progress            Print progress" << endl
       << endl;
}    

static unordered_map<string, vector<pair<int64_t, int64_t>>> load_bed(istream& bed_stream, const string& ref_prefix);
static unordered_map<string, vector<pair<int64_t, int64_t>>> find_unaligned(const PathHandleGraph* graph, int64_t max_unaligned,
                                                                            const string& ref_prefix);
static unique_ptr<MutablePathMutableHandleGraph> load_graph(istream& graph_stream);
static vector<string> &split_delims(const string &s, const string& delims, vector<string> &elems);
static void chop_path_intervals(MutablePathMutableHandleGraph* graph,
                                const unordered_map<string, vector<pair<int64_t, int64_t>>>& bed_intervals,
                                bool force_clip, bool orphan_filter,
                                const string& ref_prefix,
                                bool progress);
static pair<unordered_set<handle_t>, vector<path_handle_t>> chop_path(MutablePathMutableHandleGraph* graph,
                                                                      path_handle_t path_handle,
                                                                      const vector<pair<int64_t, int64_t>>& intervals);
static void replace_path_name_substrings(MutablePathMutableHandleGraph* graph, const vector<string>& to_replace,
                                         bool progress);
static void forwardize_paths(MutablePathMutableHandleGraph* graph, const string& ref_prefix, bool progress);
static vector<unordered_set<nid_t>> weakly_connected_components(const HandleGraph* graph);
static void drop_paths(MutablePathMutableHandleGraph* graph, const string& drop_prefix, bool progress);


// from path.cpp in vg
// Check if using subpath naming scheme.  If it is return true,
// the root path name, and the offset (false otherwise)
// formats accepted:
// name[start]  (0-based)
// name[start-end] (0-based, open-ended like BED, end defaults to 0 in above case)
static tuple<bool, string, size_t, size_t>parse_subpath_name(const string& path_name) {
    size_t tag_offset = path_name.rfind("[");
    if (tag_offset == string::npos || tag_offset + 2 >= path_name.length() || path_name.back() != ']') {
        return make_tuple(false, "", 0, 0);
    } else {
        string offset_str = path_name.substr(tag_offset + 1, path_name.length() - tag_offset - 2);
        size_t sep_offset = offset_str.find("-");
        string end_offset_str;
        if (sep_offset != string::npos && !offset_str.empty() && sep_offset < offset_str.length() - 1) {
            end_offset_str = offset_str.substr(sep_offset + 1, offset_str.length() - sep_offset - 1);
            offset_str = offset_str.substr(0, sep_offset);
        }
        size_t offset_val;
        size_t end_offset_val = 0;
        try {
           offset_val = std::stol(offset_str);
           if (!end_offset_str.empty()) {
               end_offset_val = std::stol(end_offset_str);
           }
        } catch(...) {
            return make_tuple(false, "", 0, 0);
        }
        return make_tuple(true, path_name.substr(0, tag_offset), offset_val, end_offset_val);
    }
}

// Create a subpath name (todo: make same function in vg consistent (it only includes start))
static inline string make_subpath_name(const string& path_name, size_t offset, size_t length) {
    tuple<bool, string, size_t, size_t> parsed_name = parse_subpath_name(path_name);
    size_t cur_offset = 0;
    const string* cur_name = &path_name;
    if (get<0>(parsed_name)) {
        cur_offset = get<2>(parsed_name);
        cur_name = &get<1>(parsed_name);
    }
    return *cur_name + "[" + std::to_string(cur_offset + offset) + "-" + std::to_string(cur_offset + offset + length) + "]";
}

int main(int argc, char** argv) {

    string bed_path;
    int64_t min_length = 0;
    int64_t max_unaligned = 0;
    string ref_prefix;
    size_t input_count = 0;
    bool force_clip = false;
    bool orphan_filter = true;
    bool progress = false;
    vector<string> replace_list;
    string drop_prefix;
    int c;
    optind = 1; 
    while (true) {

        static const struct option long_options[] = {
            {"help", no_argument, 0, 'h'},
            {"bed", required_argument, 0, 'b'},
            {"min-length", required_argument, 0, 'm'},
            {"max-unaligned", required_argument, 0, 'u'},
            {"ref-prefix", required_argument, 0, 'e'},
            {"force-clip", no_argument, 0, 'f'},
            {"name-replace", required_argument, 0, 'r'},
            {"no-orphan_filter", no_argument, 0, 'n'},
            {"drop-prefix", required_argument, 0, 'd'},
            {"progress", no_argument, 0, 'p'},
            {0, 0, 0, 0}
        };

        int option_index = 0;

        c = getopt_long (argc, argv, "hpb:m:u:e:fnr:d:",
                         long_options, &option_index);

        // Detect the end of the options.
        if (c == -1)
            break;

        switch (c)
        {
        case 'b':
            bed_path = optarg;
            ++input_count;
            break;
        case 'm':
            min_length = stol(optarg);
            ++input_count;
            break;
        case 'u':
            max_unaligned = stol(optarg);
            ++input_count;
            break;
        case 'e':
            ref_prefix = optarg;
            break;
        case 'f':
            force_clip = true;
            break;
        case 'n':
            orphan_filter = false;
            break;
        case 'r':
            replace_list.push_back(optarg);
            break;
        case 'd':
            drop_prefix = optarg;
            break;
        case 'p':
            progress = true;
            break;
        case 'h':
        case '?':
            /* getopt_long already printed an error message. */
            help(argv);
            exit(1);
            break;
        default:
            abort ();
        }
    }

    if (argc <= 1) {
        help(argv);
        return 1;
    }

    // Parse the positional argument
    if (optind >= argc) {
        cerr << "[clip-vg] error: too few arguments" << endl;
        help(argv);
        return 1;
    }

    if (optind != argc - 1) {
        cerr << "[clip-vg] error: too many arguments" << endl;
        help(argv);
        return 1;
    }

    if (input_count > 1) {
        cerr << "[clip-vg] error: at most one of -b, -m or -u can be used at a time" << endl;
        return 1;
    }

    if (input_count == 0 && replace_list.empty() && ref_prefix.empty()) {
        cerr <<  "[clip-vg] error: at east one of -b, -m, -u, -e or -r must be specified" << endl;
        return 1;
    }

    string graph_path = argv[optind++];
    ifstream graph_stream(graph_path);
    if (!graph_stream) {
        cerr << "[clip-vg] error: Unable to open input graph " << graph_path << endl;
        return 1;
    }    
    unique_ptr<MutablePathMutableHandleGraph> graph = load_graph(graph_stream);
    graph_stream.close();
    if (progress) {
        cerr << "[clip-vg]: Loaded graph" << endl;
    }

    unordered_map<string, vector<pair<int64_t, int64_t>>> bed_intervals;

    if (!bed_path.empty()) {
        ifstream bed_stream(bed_path);
        if (!bed_stream) {
            cerr << "[clip-vg] error: Unable to open input BED file " << bed_path << endl;
            return 1;
        }
        bed_intervals = load_bed(bed_stream, ref_prefix);
    } else if (min_length != 0) {
        // apply min length to all paths to get intervals
        graph->for_each_path_handle([&](path_handle_t path_handle) {
                string path_name = graph->get_path_name(path_handle);
                if (ref_prefix.empty() || path_name.substr(0, ref_prefix.length()) != ref_prefix) {
                    int64_t path_length = 0;
                    graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
                            path_length += graph->get_length(graph->get_handle_of_step(step_handle));
                            return path_length < min_length;
                        });
                    if (path_length < min_length) {
                        bed_intervals[path_name].push_back(make_pair(0, path_length));
                    }
                }
            });
    } else if (max_unaligned != 0) {
        // apply max unaligned length to all paths
        bed_intervals = find_unaligned(graph.get(), max_unaligned, ref_prefix);
    }
    
    if (progress) {
        size_t num_intervals = 0;
        for (auto& bi : bed_intervals) {
            num_intervals += bi.second.size();
        }
        cerr << "[clip-vg]: Loaded " << num_intervals << " BED intervals over " << bed_intervals.size() << " sequences" << endl;
    }

    if (!bed_intervals.empty()) {
        chop_path_intervals(graph.get(), bed_intervals, force_clip, orphan_filter, ref_prefix, progress);
    }

    if (!ref_prefix.empty()) {
        forwardize_paths(graph.get(), ref_prefix, progress);  
    }
    
    if (!replace_list.empty()) {
        replace_path_name_substrings(graph.get(), replace_list, progress);
    }

    if (!drop_prefix.empty()) {
        drop_paths(graph.get(), drop_prefix, progress);
    }

    dynamic_cast<SerializableHandleGraph*>(graph.get())->serialize(cout);

    return 0;
}

unordered_map<string, vector<pair<int64_t, int64_t>>> load_bed(istream& bed_stream, const string& ref_prefix) {
    // load bed
    unordered_map<string, vector<pair<int64_t, int64_t>>> intervals;
    string buffer;
    while (getline(bed_stream, buffer)) {
        vector<string> toks;
        split_delims(buffer, "\t\n", toks);
        if (toks.size() >= 3) {
            string& name = toks[0];
            if (ref_prefix.empty() || name.substr(0, ref_prefix.length()) != ref_prefix) {
                int64_t start = stol(toks[1]);
                int64_t end = stol(toks[2]);
                intervals[name].push_back(make_pair(start, end));
            }
        }
    }
    // verify bed
    for (auto& seq_intervals : intervals) {
        sort(seq_intervals.second.begin(), seq_intervals.second.end(),
             [](const pair<int64_t, int64_t>& b1, const pair<int64_t, int64_t>& b2) {
                 return b1.first < b2.first || (b1.first == b2.first && b1.second < b2.second);
             });
        for (size_t i = 1; i < seq_intervals.second.size(); ++i) {
            if (seq_intervals.second[i].first < seq_intervals.second[i-1].second) {
                cerr << "Overlapping bed intervals found:\n"
                     << " " << seq_intervals.first << "\t"
                     << seq_intervals.second[i-1].first << "\t"
                     << seq_intervals.second[i-1].second << endl
                     << " " << seq_intervals.first << "\t"
                     << seq_intervals.second[i].first << "\t"
                     << seq_intervals.second[i].second << endl
                     << "These are not supported.  Please clean up (ex with bedools merge) first" << endl;
                exit(1);
            }
        }
    }
    return intervals;
}

unordered_map<string, vector<pair<int64_t, int64_t>>> find_unaligned(const PathHandleGraph* graph, int64_t max_unaligned,
                                                                     const string& ref_prefix) {
    unordered_map<string, vector<pair<int64_t, int64_t>>> intervals;

    graph->for_each_path_handle([&](path_handle_t path_handle) {
            string path_name = graph->get_path_name(path_handle);
            if (ref_prefix.empty() || path_name.substr(0, ref_prefix.length()) != ref_prefix) {
                int64_t offset = 0;
                int64_t start = -1;
                graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
                        handle_t handle = graph->get_handle_of_step(step_handle);
                        int64_t len = (int64_t)graph->get_length(handle);
                        bool aligned = false;
                        graph->for_each_step_on_handle(handle, [&](step_handle_t step_handle_2) {
                                if (graph->get_path_handle_of_step(step_handle_2) != path_handle) {
                                    aligned = true;
                                }
                                return !aligned;
                            });
                        // start an unaligned interval
                        if (start < 0 && aligned == false) {
                            start = offset;
                        }
                        // end an unaligned interval
                        if (aligned == true) {
                            if (start >= 0 && offset + len - start > max_unaligned) {
                                intervals[path_name].push_back(make_pair(start, offset));
                            }
                            start = -1;
                        }
                        offset += len;
                    });
                if (start >= 0 && offset - start > max_unaligned) {
                    intervals[path_name].push_back(make_pair(start, offset));
                }
            }
        });
    return intervals;
}


unique_ptr<MutablePathMutableHandleGraph> load_graph(istream& graph_stream) {

    char magic_bytes[4];
    graph_stream.read(magic_bytes, 4);
    uint32_t magic_number = ntohl(*((uint32_t*) magic_bytes));
    graph_stream.clear();
    graph_stream.seekg(0, ios::beg);

    MutablePathMutableHandleGraph* graph;
    if (magic_number == PackedGraph().get_magic_number()) {
        graph = new PackedGraph();
    } else if (magic_number == HashGraph().get_magic_number()) {
        graph = new HashGraph();
    } else if (magic_number == ODGI().get_magic_number()) {
        graph = new ODGI();
    } else {
        cerr << "Unable to parse input graph with magic number " << magic_number << endl;
        exit(1);
    }
    dynamic_cast<SerializableHandleGraph*>(graph)->deserialize(graph_stream);

    return unique_ptr<MutablePathMutableHandleGraph>(graph);
}

vector<string> &split_delims(const string &s, const string& delims, vector<string> &elems) {
    size_t start = string::npos;
    for (size_t i = 0; i < s.size(); ++i) {
        if (delims.find(s[i]) != string::npos) {
            if (start != string::npos && i > start) {
                elems.push_back(s.substr(start, i - start));
            }
            start = string::npos;
        } else if (start == string::npos) {
            start = i;
        }
    }
    if (start != string::npos && start < s.size()) {
        elems.push_back(s.substr(start, s.size() - start));
    }
    return elems;
}

void chop_path_intervals(MutablePathMutableHandleGraph* graph,
                         const unordered_map<string, vector<pair<int64_t, int64_t>>>& bed_intervals,
                         bool force_clip, bool orphan_filter, const string& ref_prefix,
                         bool progress) {

    // keep some stats to print
    size_t chopped_paths = 0;
    size_t chopped_nodes = 0;
    size_t chopped_bases = 0;
    
    // careful not to iterate and chop, as we could hit new subpaths made
    vector<path_handle_t> path_handles;    
    graph->for_each_path_handle([&](path_handle_t path_handle) {
            path_handles.push_back(path_handle);
        });

    // when force_clip is true, store handles here to given them second chance at destruction
    // after all paths are deleted
    unordered_set<nid_t> to_destroy;
    // newly created subpaths
    vector<path_handle_t> subpaths;
    
    for (auto path_handle : path_handles) {
        string path_name = graph->get_path_name(path_handle);
        auto it = bed_intervals.find(path_name);
        bool was_chopped = false;
        if (it != bed_intervals.end()) {
            if (progress) {
                cerr << "[clip-vg]: Clipping " << it->second.size() << " intervals from path " << path_name << endl;
            }
            auto chopped_handles_subpaths = chop_path(graph, path_handle, it->second);
            auto& chopped_handles = chopped_handles_subpaths.first;
            subpaths.insert(subpaths.end(), chopped_handles_subpaths.second.begin(), chopped_handles_subpaths.second.end());
            if (!chopped_handles.empty()) {
#ifdef debug
                cerr << "destroying path " << graph->get_path_name(path_handle) << endl;
#endif
                graph->destroy_path(path_handle);
                for (handle_t handle : chopped_handles) {
                    if (force_clip) {
                        to_destroy.insert(graph->get_id(handle));
                    } else {
                        vector<step_handle_t> steps = graph->steps_of_handle(handle);
                        bool aligned = false;
                        for (size_t i = 0; i < steps.size() && !aligned; ++i) {
                            string other_path_name = graph->get_path_name(graph->get_path_handle_of_step(steps[i]));
                            if (path_name.substr(0, path_name.rfind("[")) !=
                                other_path_name.substr(0, other_path_name.rfind("["))) {
                                aligned = true;
                            }
                        }
                        if (!aligned) {
                            chopped_bases += graph->get_length(handle);
                            was_chopped = true;
                            ++chopped_nodes;
                            to_destroy.insert(graph->get_id(handle));
                        } else {
                            cerr << "[clip-vg]: Unable to clip node " << graph->get_id(handle) << ":" << graph->get_is_reverse(handle)
                                 << " in path " << path_name << " because it is found in the following other paths:\n";
                            for (step_handle_t step : graph->steps_of_handle(handle)) {
                                cerr <<"\t" << graph->get_path_name(graph->get_path_handle_of_step(step)) << endl;
                            }
                            cerr << " Use the -f option to not abort in this case" << endl;
                            exit(1);
                        }
                    }
                }
            }
        }
        if (was_chopped) {
            ++chopped_paths;
        }
    }

    // trim out fragments between clipped regions that would otherwise be left disconnected from the graph
    size_t removed_subpath_count = 0;
    size_t removed_subpath_base_count = 0;
    size_t removed_component_count = 0;
    size_t removed_component_base_count = 0;
    if (orphan_filter) {
        for (path_handle_t subpath_handle : subpaths) {
            bool connected = false;
            graph->for_each_step_in_path(subpath_handle, [&](step_handle_t step_handle) {
                    connected = graph->steps_of_handle(graph->get_handle_of_step(step_handle)).size() > 1;
                    return !connected;
                });
            if (!connected) {
                graph->for_each_step_in_path(subpath_handle, [&](step_handle_t step_handle) {
                        handle_t handle = graph->get_handle_of_step(step_handle);
                        to_destroy.insert(graph->get_id(handle));
                        removed_subpath_base_count += graph->get_length(handle);
                    });
                graph->destroy_path(subpath_handle);
                if (progress) {
                    cerr << "[clip-vg]: Removing orphaned subpath " << graph->get_path_name(subpath_handle) << endl;
                }
                ++removed_subpath_count;
            }
        }

        // use the reference path prefix (if given) to clip out components that aren't anchored to it
        // (this would take care of above filter, but we leave that one as it's not dependent on path name)
        if (!ref_prefix.empty()) {
            vector<unordered_set<nid_t>> components = weakly_connected_components(graph);
            for (auto& component : components) {
                bool ref_anchored = false;
                for (auto ni = component.begin(); !ref_anchored && ni != component.end(); ++ni) {
                    vector<step_handle_t> steps = graph->steps_of_handle(graph->get_handle(*ni));
                    for (size_t si = 0; !ref_anchored && si < steps.size(); ++si) {
                        string step_path_name = graph->get_path_name(graph->get_path_handle_of_step(steps[si]));
                        if (step_path_name.substr(0, ref_prefix.length()) == ref_prefix) {
                            ref_anchored = true;
                        }
                    }
                }
                if (!ref_anchored) {
                    ++removed_component_count;
                    for (auto node_id : component) {
                        handle_t node_handle = graph->get_handle(node_id);
                        removed_component_base_count += graph->get_length(node_handle);
                        // destroy here instead of adding to to_destroy, becuase we don't care
                        // if there are paths or not (so don't require -f)
                        dynamic_cast<DeletableHandleGraph*>(graph)->destroy_handle(node_handle);
                        if (to_destroy.count(node_id)) {
                            to_destroy.erase(node_id);
                        }
                    }
                }
            }
        }
    }
    
    for (nid_t nid : to_destroy) {
        assert(graph->has_node(nid));
        handle_t handle = graph->get_handle(nid);
        if (graph->steps_of_handle(handle).empty()) {
            chopped_bases += graph->get_length(handle);
            ++chopped_nodes;
            dynamic_cast<DeletableHandleGraph*>(graph)->destroy_handle(handle);
#ifdef debug
            cerr << "force destroying handle " << graph->get_id(handle) << ":" << graph->get_is_reverse(handle) << endl;
#endif
        }
    }
    
    if (progress) {
        cerr << "[clip-vg]: Clipped "
             << chopped_bases << " bases from "
             << chopped_nodes << " nodes";
        if (!force_clip) {
            cerr << " in " << chopped_paths << " paths";
        }
        cerr << endl;        
        if (removed_subpath_count > 0) {
            cerr << "[clip-vg]: " << removed_subpath_count << " orphaned subpaths were removed with total "
                 << removed_subpath_base_count << " bases" << endl;
        }
        if (removed_component_count > 0) {
            cerr << "[clip-vg]: " << removed_component_count << " orphaned connected components were removed with total "
                 << removed_component_base_count << " bases" << endl;
        }
    }
}

pair<unordered_set<handle_t>, vector<path_handle_t>> chop_path(MutablePathMutableHandleGraph* graph,
                                                               path_handle_t path_handle,
                                                               const vector<pair<int64_t, int64_t>>& intervals) {

    // get the breakpoints
    set<int64_t> breakpoints;
    for (const pair<int64_t, int64_t>& interval : intervals) {
        breakpoints.insert(interval.first);
        breakpoints.insert(interval.second); // we're cutting before offset, so the open coordinate is what we want
    }

    // to be safe, don't cut and iterate at the same time, so load up steps here
    vector<handle_t> steps;
    graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
            steps.push_back(graph->get_handle_of_step(step_handle));
        });
    
    // cut the nodes to ensure breakpoints at node boundaries
    int64_t offset = 0;
    for (auto handle : steps) {
        int64_t len = graph->get_length(handle);
        // find breakpoints in node
        vector<size_t> cut_points;
        for (auto i = breakpoints.lower_bound(offset); i != breakpoints.end() && *i - offset < len; ++i) {
            int64_t cut_point = *i - offset;
            // libbdsg is buggy and can't accept cutpoints on ends on reverse strand
            if (cut_point > 0 && cut_point < len) {
                cut_points.push_back(cut_point);
            }
        }
        // chop the node
        if (!cut_points.empty()) {
#ifdef debug
            cerr << "dividing node_id=" << graph->get_id(handle) << ":" << graph->get_is_reverse(handle) << " seq=" << graph->get_sequence(handle)
                 << " for path " << graph->get_path_name(path_handle) << " at cut points:";
            for (auto cp : cut_points) {
                cerr << " " << cp;
            }
            cerr << endl;
#endif        
            size_t total_pieces_length = 0;
            vector<handle_t> pieces = graph->divide_handle(handle, cut_points) ;
            for (size_t i = 0; i < pieces.size(); ++i) {
                handle_t& piece = pieces[i];
                size_t piece_length = graph->get_length(piece);
                if (i == 0) {
                    assert(piece_length == cut_points[0]);
                } else if (i < pieces.size() - 1) {
                    assert(piece_length == cut_points[i] - cut_points[i-1]);
                }
                total_pieces_length += piece_length;
#ifdef debug
                cerr << " piece " << graph->get_id(piece) << ":" << graph->get_is_reverse(piece) << " " << graph->get_sequence(piece)
                     << " tlen=" << total_pieces_length << "/" << len << endl;
#endif
            }
            // bugs in divide-handle turning out to be a real issue.  add this sanity check to catch them early.
            assert(total_pieces_length == (size_t)len);
        }
        offset += len;
    }
    
    steps.clear();
    int64_t original_path_length = offset;
    unordered_set<handle_t> chopped_handles;
    offset = 0;
    step_handle_t current_step = graph->path_begin(path_handle);
#ifdef debug
    cerr << "init step to " << graph->get_id(graph->get_handle_of_step(current_step)) << ":" << graph->get_is_reverse(graph->get_handle_of_step(current_step))
         << " seq=" <<graph->get_sequence(graph->get_handle_of_step(current_step)) << endl;
#endif
    vector<path_handle_t> subpaths;
    
    // cut out a subpath and make a new path out of it
    function<void(int64_t)> cut_to = [&](int64_t end_offset) {
#ifdef debug
        cerr << "\ncut_to " << end_offset << " where current offset is " << offset << endl;
#endif
        vector<handle_t> steps;
        int64_t start_offset = offset;
        int64_t path_length = 0;
        while (offset < end_offset && current_step != graph->path_end(path_handle)) {
            handle_t handle = graph->get_handle_of_step(current_step);
            steps.push_back(handle);
            offset += graph->get_length(handle);
            current_step = graph->get_next_step(current_step);
            path_length += graph->get_length(handle);
        }
#ifdef debug
        cerr << "start offset=" << start_offset << " path length=" << path_length << " end offset=" << end_offset << endl;
#endif
        assert(start_offset + path_length == end_offset);

        if (path_length > 0) {
            path_handle_t subpath_handle = graph->create_path_handle(make_subpath_name(graph->get_path_name(path_handle), start_offset, path_length));
            for (auto step : steps) {
#ifdef debug
                cerr << " pushing subpath step " << graph->get_id(step) << ":" << graph->get_is_reverse(step)
                     << " len=" << graph->get_length(step) <<  " to " << graph->get_path_name(subpath_handle) << endl;
#endif
                graph->append_step(subpath_handle, step);
            }
            subpaths.push_back(subpath_handle);
        }
    };

    
    for (size_t i = 0; i < intervals.size(); ++i) {
        if (intervals[i].first > offset) {
            // cut everythign left of the interval
            cut_to(intervals[i].first);
        }
        // scan past the interval
        while (offset < intervals[i].second && current_step != graph->path_end(path_handle)) {
            handle_t handle = graph->get_handle_of_step(current_step);
            offset += graph->get_length(handle);
            current_step = graph->get_next_step(current_step);
#ifdef debug
            cerr << "adding to delete set: " << graph->get_id(handle) << endl;
#endif
            chopped_handles.insert(handle);
        }
    }

    // cut the last bit
    if (offset < original_path_length) {
        cut_to(original_path_length);
    }
    
    return make_pair(chopped_handles, subpaths);    
}

void replace_path_name_substrings(MutablePathMutableHandleGraph* graph, const vector<string>& to_replace,
                                  bool progress) {
    // parse the strings
    vector<pair<string, string>> replace;
    for (const string& repstring : to_replace) {
        size_t sep = repstring.find('>');
        if (sep == string::npos || sep == 0 || sep == repstring.length() - 1) {
            cerr << "[clip-vg]: Unable to find separator '>' in " << repstring << ". Replacement must be"
                 << " specified with \"s1>s2\"" << endl;
            exit(1);
        }
        replace.push_back(make_pair(repstring.substr(0, sep), repstring.substr(sep + 1)));
        if (replace.back().first == replace.back().second) {
            replace.pop_back();
        }
    }

    size_t replacement_count = 0;
    size_t path_count = 0;
    // take care to not modify path handles while iterating path handles, just in case
    vector<string> path_names;
    graph->for_each_path_handle([&](path_handle_t path_handle) {
            path_names.push_back(graph->get_path_name(path_handle));
        });
    for (string& path_name : path_names) {
        path_handle_t path_handle = graph->get_path_handle(path_name);
        bool changed = false;
        for (auto& rep : replace) {
            size_t p = path_name.find(rep.first);
            if (p != string::npos) {
                path_name.replace(p, rep.first.length(), rep.second);
                ++replacement_count;
                changed = true;
            }
        }
        if (changed) {
            ++path_count;
            if (graph->has_path(path_name)) {
                cerr << "[clip-vg] error: cannot change name of path from " << graph->get_path_name(path_handle) << " to "
                     << path_name << " because the latter already exists in the graph" << endl;
                exit(1);
            }
            path_handle_t new_path_handle = graph->create_path_handle(path_name, graph->get_is_circular(path_handle));
            graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
                    graph->append_step(new_path_handle, graph->get_handle_of_step(step_handle));
                });
            graph->destroy_path(path_handle);
        }
    }
    
    if (progress) {
        cerr << "[clip-vg]: Replaced " << replacement_count << " substrings in " << path_count << " path names" << endl;
    }
}

void forwardize_paths(MutablePathMutableHandleGraph* graph, const string& ref_prefix, bool progress) {
    
    graph->for_each_path_handle([&](path_handle_t path_handle) {
            string path_name = graph->get_path_name(path_handle);
            if (path_name.substr(0, ref_prefix.length()) == ref_prefix) {
                size_t fw_count = 0;
                size_t total_steps = 0;
                graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
                        handle_t handle = graph->get_handle_of_step(step_handle);
                        if (graph->get_is_reverse(handle)) {
                            handle_t flipped_handle = graph->create_handle(graph->get_sequence(handle));
			    graph->follow_edges(handle, true, [&](handle_t prev_handle) {
                                    if (graph->get_id(prev_handle) != graph->get_id(handle)) {
                                        graph->create_edge(prev_handle, flipped_handle);
                                    }
			      });
			    graph->follow_edges(handle, false, [&](handle_t next_handle) {
                                    if (graph->get_id(handle) != graph->get_id(next_handle)) {
                                        graph->create_edge(flipped_handle, next_handle);
                                    }
			      });
                            // self-loop cases we punted on above:
                            if (graph->has_edge(handle, handle)) {
                                graph->create_edge(flipped_handle, flipped_handle);
                            }
                            if (graph->has_edge(handle, graph->flip(handle))) {
                                graph->create_edge(flipped_handle, graph->flip(flipped_handle));                                
                            }
                            if (graph->has_edge(graph->flip(handle), handle)) {
                                graph->create_edge(graph->flip(flipped_handle), flipped_handle);
                            }
                            vector<step_handle_t> steps = graph->steps_of_handle(handle);
                            size_t ref_count = 0;
                            for (step_handle_t step : steps) {
                                if (graph->get_path_handle_of_step(step) == path_handle) {
                                    ++ref_count;
                                }
                                step_handle_t next_step = graph->get_next_step(step);
                                handle_t new_handle = graph->get_is_reverse(graph->get_handle_of_step(step)) ? flipped_handle :
                                    graph->flip(flipped_handle);
                                graph->rewrite_segment(step, next_step, {new_handle});
                            }
                            if (ref_count > 1) {
                                cerr << "[clip-vg] error: Cycle detected in reference path " << path_name << " at node " << graph->get_id(handle) << endl;
                                exit(1);
                            }
                            ++fw_count;
                            assert(graph->steps_of_handle(handle).empty());
                            dynamic_cast<DeletableHandleGraph*>(graph)->destroy_handle(handle);
                        }
                        ++total_steps;
                    });
                if (fw_count > 0 && progress) {
                    cerr << "[clip-vg]: Forwardized " << fw_count << " / " << total_steps << " steps in reference path " << path_name << endl;
                }
            }
        });

    // do a check just to be sure
    graph->for_each_path_handle([&](path_handle_t path_handle) {
            string path_name = graph->get_path_name(path_handle);
            if (path_name.substr(0, ref_prefix.length()) == ref_prefix) {
                graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
                        handle_t handle = graph->get_handle_of_step(step_handle);
                        if (graph->get_is_reverse(handle)) {
                            cerr << "[clip-vg] error: Failed to fowardize node " << graph->get_id(handle) << " in path " << path_name << endl;
                            exit(1);
                        }
                    });
            }
        });
}

// this is pasted from libhandlegraph
// todo: update libhandlegraph to version that contains algorithms!!!
vector<unordered_set<nid_t>> weakly_connected_components(const HandleGraph* graph) {
    vector<unordered_set<nid_t>> to_return;
    
    // This only holds locally forward handles
    unordered_set<handle_t> traversed;
    
    graph->for_each_handle([&](const handle_t& handle) {
        
        // Only think about it in the forward orientation
        auto forward = graph->forward(handle);
        
        if (traversed.count(forward)) {
            // Already have this node, so don't start a search from it.
            return;
        }
        
        // The stack only holds locally forward handles
        vector<handle_t> stack{forward};
        to_return.emplace_back();
        while (!stack.empty()) {
            handle_t here = stack.back();
            stack.pop_back();
            
            traversed.insert(here);
            to_return.back().insert(graph->get_id(here));
            
            // We have a function to handle all connected handles
            auto handle_other = [&](const handle_t& other) {
                // Again, make it forward
                auto other_forward = graph->forward(other);
                
                if (!traversed.count(other_forward)) {
                    stack.push_back(other_forward);
                }
            };
            
            // Look at edges in both directions
            graph->follow_edges(here, false, handle_other);
            graph->follow_edges(here, true, handle_other);
            
        }
    });
    return to_return;
}

// this was written to filter out minigraph-only nodes from the graph (and all the minigraph paths)
// this used to get done by hal2vg, but it's useful to keep them around so that -u option will work
// better.  ie this way, a path that's private to a sample but still in the minigraph will be kept
// because it aligns to two paths whereas if minigraph paths weren't in, it'd be deleted
void drop_paths(MutablePathMutableHandleGraph* graph, const string& drop_prefix, bool progress) {

    unordered_set<nid_t> to_destroy;
    size_t removed_path_count = 0;
    size_t removed_base_count = 0;

    // careful not to iterate and chop, as we could hit new subpaths made
    vector<path_handle_t> path_handles;    
    graph->for_each_path_handle([&](path_handle_t path_handle) {
            path_handles.push_back(path_handle);
        });

    for (path_handle_t& path_handle : path_handles) {
        string path_name = graph->get_path_name(path_handle);
        if (path_name.compare(0, drop_prefix.length(), drop_prefix) == 0) {
            // we've found a path with the given prefix: now destroy all handles that don't touch
            // any path *without* the prefix
            graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
                    handle_t handle = graph->get_handle_of_step(step_handle);
                    bool has_other_path = false;
                    graph->for_each_step_on_handle(handle, [&](step_handle_t step_handle_2) {
                            path_handle_t other_path_handle = graph->get_path_handle_of_step(step_handle_2);
                            if (other_path_handle != path_handle &&
                                graph->get_path_name(other_path_handle).compare(0, drop_prefix.length(), drop_prefix) != 0) {
                                has_other_path = true;
                                return false;
                            }
                            return true;
                        });
                    if (!has_other_path) {
                        to_destroy.insert(graph->get_id(handle));
                    }
                });
            // detroy the path
            graph->destroy_path(path_handle);
            removed_path_count++;
        }
        }

    // destroy the nodes
    size_t removed_node_count = to_destroy.size();
    for (nid_t node_id : to_destroy) {
        handle_t node_handle = graph->get_handle(node_id);
        removed_base_count += graph->get_length(node_handle);
        dynamic_cast<DeletableHandleGraph*>(graph)->destroy_handle(node_handle);        
    }

    if (progress) {
        cerr << "[clip-vg]: Drop prefix removed " << removed_base_count << " bases from " << removed_node_count << " nodes in " << removed_path_count << " paths" << endl;        
    }
}
