clip-vg.cpp

// 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 <limits>
#include <cmath>

#include "bdsg/packed_graph.hpp"
#include "bdsg/hash_graph.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
       << "    -a, --anchor PREFIX       If set, consider regions not aligned to a path with PREFIX unaligned (with -u)" << endl
       << "    -k, --flank N             Extend each clipped interval outward by up to N bp for as long as" << endl
       << "                              unaligned sequence stays dense (see -T), using the same test for" << endl
       << "                              unaligned that -u does.  Removes the fringe left where an aligner" << endl
       << "                              extended anchors into a repeat but not far enough to be clipped." << endl
       << "                              N is the main control: inside a repeat the gate seldom closes" << endl
       << "                              on its own, so the cap is what decides how much goes." << endl
       << "    -T, --flank-threshold F   Keep extending while more than this fraction of bases are" << endl
       << "                              unaligned (0-1).  Negative means calibrate against this graph," << endl
       << "                              which is strongly advisable: how much unaligned sequence a" << endl
       << "                              pangenome carries varies several-fold between graphs and is" << endl
       << "                              not predictable from any input property [-1]" << endl
       << "    -e, --ref-prefix STR      Forwardize (but don't clip) paths whose name begins with STR" << endl
       << "    -F, --forwardize-nonref   Also forwardize any node that no path visits forward.  Needs -e." << endl
       << "                              This mints new node ids, so it must not be used on a graph whose" << endl
       << "                              id space has to stay compatible with one produced earlier." << endl
       << "    -c, --allow-cycle         Do not fail with error when reference cycle detected" << 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
       << "    -L, --leave-aligned       When used in conjunction with -d, paths are prserved if they align to a non-dropped path" << endl 
       << "    -o, --out-bed FILE        Save all clipped intervals here" << 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, const string& anchor_prefix,
                                                                            const unordered_set<nid_t>& anchor_nodes);
static void build_anchor_nodes(const PathHandleGraph* graph, const string& anchor_prefix,
                               unordered_set<nid_t>& anchor_nodes_out);
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 allow_ref_cycles, bool progress);
static void forwardize_nonref_paths(MutablePathMutableHandleGraph* graph, const string& ref_prefix, bool progress);
static void flip_node(MutablePathMutableHandleGraph* graph, nid_t node_id);
static vector<unordered_set<nid_t>> weakly_connected_components(const HandleGraph* graph);
static void drop_paths(MutablePathMutableHandleGraph* graph, const string& drop_prefix, bool leave_aligned, bool progress);

// Is this step aligned?  With an anchor prefix that means the node sits on one of those paths;
// without one it means some other path visits the node.  find_unaligned uses this to cut intervals
// and extend_flanks uses it to decide how far to extend them, so they never disagree.
static bool step_is_aligned(const PathHandleGraph* graph, handle_t handle, path_handle_t path_handle,
                            const unordered_set<nid_t>& anchor_nodes, bool have_anchor_prefix) {
    if (anchor_nodes.count(graph->get_id(handle))) {
        return true;
    }
    if (have_anchor_prefix) {
        return false;
    }
    bool aligned = false;
    graph->for_each_step_on_handle(handle, [&](step_handle_t other) {
            if (graph->get_path_handle_of_step(other) != path_handle) {
                aligned = true;
            }
            return !aligned;
        });
    return aligned;
}

static void extend_flanks(const PathHandleGraph* graph,
                          unordered_map<string, vector<pair<int64_t, int64_t>>>& intervals,
                          int64_t max_flank, double threshold,
                          const unordered_set<nid_t>& anchor_nodes, bool have_anchor_prefix,
                          bool progress);

static unordered_map<string, vector<pair<int64_t, int64_t>>> get_path_intervals(const PathHandleGraph* graph);

static unordered_map<string, vector<pair<int64_t, int64_t>>> get_clipped_intervals(
    const unordered_map<string, vector<pair<int64_t, int64_t>>>& input_intervals,
    const unordered_map<string, vector<pair<int64_t, int64_t>>>& output_intervals);

// A path's name with any subrange stripped off, so that an input path and the subpaths clipped
// out of it share a key.  get_path_name() spells the subrange out as "name[start-end]", which would
// file every output subpath under a key no input path has.
static inline string strip_subpath_name(const string& path_name) {
    PathSense sense;
    string sample;
    string locus;
    size_t haplotype;
    size_t phase_block;
    subrange_t subrange;
    PathMetadata::parse_path_name(path_name, sense, sample, locus, haplotype, phase_block, subrange);
    return PathMetadata::create_path_name(sense, sample, locus, haplotype, phase_block,
                                          PathMetadata::NO_SUBRANGE);
}

// 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) {
    PathSense sense;
    string sample;
    string locus;
    size_t haplotype;
    size_t phase_block;
    subrange_t subrange;
    PathMetadata::parse_path_name(path_name, sense, sample, locus, haplotype, phase_block, subrange);
    subrange.first = subrange != PathMetadata::NO_SUBRANGE ? subrange.first : 0;
    subrange.first += offset;
    subrange.second = subrange.first + length;
    return PathMetadata::create_path_name(sense, sample, locus, haplotype, phase_block, subrange);
}

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

    string bed_path;
    // shared by find_unaligned and extend_flanks so both agree on what unaligned means
    unordered_set<nid_t> anchor_nodes;
    int64_t min_length = 0;
    int64_t max_unaligned = 0;
    int64_t flank = 0;
    double flank_threshold = -1.;
    bool flank_threshold_set = false;
    string anchor_prefix;
    string ref_prefix;
    bool forwardize_nonref = false;
    bool allow_ref_cycles = false;
    size_t input_count = 0;
    bool force_clip = false;
    bool orphan_filter = true;
    bool progress = false;
    vector<string> replace_list;
    string drop_prefix;
    bool leave_aligned_drop_paths = false;
    string out_bed_path;
    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'},
            {"anchor", required_argument, 0, 'a'},
            {"flank", required_argument, 0, 'k'},
            {"flank-threshold", required_argument, 0, 'T'},
            {"ref-prefix", required_argument, 0, 'e'},
            {"forwardize-nonref", no_argument, 0, 'F'},
            {"allow-cycle", no_argument, 0, 'c'},
            {"force-clip", no_argument, 0, 'f'},
            {"name-replace", required_argument, 0, 'r'},
            {"no-orphan-filter", no_argument, 0, 'n'},
            // the spelling that shipped, kept working so nobody's script breaks
            {"no-orphan_filter", no_argument, 0, 'n'},
            {"drop-prefix", required_argument, 0, 'd'},
            {"leave-aligned", no_argument, 0, 'L'},
            {"out-bed", required_argument, 0, 'o'},
            {"progress", no_argument, 0, 'p'},
            {0, 0, 0, 0}
        };

        int option_index = 0;

        c = getopt_long (argc, argv, "hpb:m:u:a:k:T:e:Fcfnr:d:Lo:",
                         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 'a':
            anchor_prefix = optarg;
            break;
        case 'k':
            flank = stol(optarg);
            break;
        case 'T':
            flank_threshold = stod(optarg);
            flank_threshold_set = true;
            break;
        case 'e':
            ref_prefix = optarg;
            break;
        case 'F':
            forwardize_nonref = true;
            break;
        case 'c':
            allow_ref_cycles = true;
            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 'L':
            leave_aligned_drop_paths = true;
            break;
        case 'o':
            out_bed_path = 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 (forwardize_nonref && ref_prefix.empty()) {
        cerr << "[clip-vg] error: -F/--forwardize-nonref requires -e/--ref-prefix" << endl;
        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;
    }
    if (!anchor_prefix.empty() && max_unaligned <= 0 && flank <= 0) {
        // -a used to matter only to -u.  It now also picks the test -k extends on, so it is
        // meaningful alongside -b, where the intervals come from a file but the gate still has to
        // decide what counts as unaligned.
        cerr << "[clip-vg] error: -a cannot be used without -u or -k" << endl;
        return 1;
    }
    
    if (leave_aligned_drop_paths && drop_prefix.empty()) {
        cerr << "[clip-vg] error: -L can only be used with -d" << endl;
        return 1;
    }

    // -0.0 is not less than zero, so it would slip past the "negative means calibrate" test in
    // extend_flanks() and land on threshold 0 -- the most aggressive setting there is, where a long
    // node scores 0 rather than negative and the walk can never die on clean sequence.  Python's
    // str(-0.0) is "-0.0", so a config carrying it reaches us in this form.  Normalise it here.
    if (std::signbit(flank_threshold)) {
        flank_threshold = -1.;
    }
    if (std::isnan(flank_threshold)) {
        cerr << "[clip-vg] error: -T/--flank-threshold is not a number.  Every node would score NaN,"
             << " no comparison against the running maximum would ever be true, and -k would do"
             << " nothing at all without saying so.  Use a negative value to calibrate." << endl;
        return 1;
    }
    if (flank_threshold >= 1.) {
        cerr << "[clip-vg] error: -T/--flank-threshold must be less than 1: at 1 or above every"
             << " node scores negative and no flank is ever trimmed.  Use a negative value to"
             << " calibrate against the graph." << endl;
        return 1;
    }

    // -k reshapes clipped intervals, so it needs a source of them.  -m makes exactly one interval
    // per path, spanning the whole path, which leaves nothing outside to extend into.  Warn rather
    // than fail: a wrapper may pass -k unconditionally and turn it off with -k 0, and failing a
    // multi-hour job over that would be worse.
    if (flank > 0) {
        if (input_count == 0) {
            cerr << "[clip-vg] warning: -k/--flank needs clipped intervals to act on, but none of"
                 << " -b, -m or -u was given" << endl;
        } else if (min_length != 0) {
            cerr << "[clip-vg] warning: -k/--flank has no effect with -m/--min-length, which clips"
                 << " whole paths" << endl;
        }
    }
    if (flank_threshold_set && flank <= 0) {
        cerr << "[clip-vg] warning: -T/--flank-threshold only applies with -k/--flank" << endl;
    }
    // 0 is in the documented range and is a coherent thing to ask for -- extend unconditionally,
    // which is the no-gate control -- but it does not read that way, so say what it does.  An
    // aligned base scores exactly 0 at this threshold rather than negative, so the running total
    // never falls and the walk cannot die on aligned sequence.
    if (flank > 0 && flank_threshold == 0.) {
        cerr << "[clip-vg] warning: -T 0 gates nothing: aligned sequence scores 0 rather than"
             << " negative, so every trim runs out to the -k cap.  Pass a positive threshold to"
             << " gate on unaligned density, or a negative one to calibrate." << endl;
    }

    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;
    }

    // Both the interval search and the flank gate ask the same question of a node, so build the
    // table once here.  -k can be used with -b, where find_unaligned never runs, and an empty
    // table would make every node look unaligned.
    build_anchor_nodes(graph.get(), anchor_prefix, anchor_nodes);

    // Catch a mistyped -e here, before anything has been clipped.  Unaligned regions are measured
    // against the reference, so a prefix that matches no path makes every base look unaligned and
    // the orphan filter then removes the entire graph -- silently, with exit 0.
    if (!ref_prefix.empty()) {
        bool found_ref_path = false;
        graph->for_each_path_handle([&](path_handle_t path_handle) {
                if (graph->get_path_name(path_handle).substr(0, ref_prefix.length()) == ref_prefix) {
                    found_ref_path = true;
                }
                return !found_ref_path;
            });
        if (!found_ref_path) {
            cerr << "[clip-vg] error: No path name begins with \"" << ref_prefix
                 << "\" given with -e/--ref-prefix" << endl;
            return 1;
        }
    }

    unordered_map<string, vector<pair<int64_t, int64_t>>> input_graph_intervals;
    if (!out_bed_path.empty()) {
        input_graph_intervals = get_path_intervals(graph.get());
        if (progress) {
            cerr << "[clip-vg]: Graph has " << input_graph_intervals.size() << " paths." << 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
        if (progress) {
            cerr << "[clip-vg]: Finding unaligned intervals >= " << max_unaligned
                 << " using anchor prefix " << anchor_prefix << " and ref prefix " << ref_prefix << endl;
        }
        bed_intervals = find_unaligned(graph.get(), max_unaligned, ref_prefix, anchor_prefix,
                                       anchor_nodes);
    }
    
    if (flank > 0 && !bed_intervals.empty()) {
        extend_flanks(graph.get(), bed_intervals, flank, flank_threshold, anchor_nodes,
                      !anchor_prefix.empty(), progress);

        // Mott's rule walks past a locally aligned stretch when unaligned sequence resumes beyond
        // it, so two intervals with a short gap between them can each extend across it and end up
        // overlapping.  chop_path requires sorted disjoint intervals, so collapse those.  Measured
        // 3 collapses on a 118-haplotype chrY and 93 on chr21, so this is not a corner case.
        size_t num_overlaps = 0;
        for (auto& bi : bed_intervals) {
            vector<pair<int64_t, int64_t>>& intervals = bi.second;
            if (intervals.size() < 2) {
                continue;
            }
            sort(intervals.begin(), intervals.end());
            vector<pair<int64_t, int64_t>> merged;
            merged.push_back(intervals[0]);
            for (size_t i = 1; i < intervals.size(); ++i) {
                // Only an actual overlap.  gap == 0 is book-ended, which load_bed has always
                // accepted and chop_path has always handled -- merging those would drop a
                // breakpoint and renumber the nodes a plain -b run has always emitted.
                if (intervals[i].first < merged.back().second) {
                    ++num_overlaps;
                    merged.back().second = max(merged.back().second, intervals[i].second);
                } else {
                    merged.push_back(intervals[i]);
                }
            }
            swap(intervals, merged);
        }
        if (progress && num_overlaps > 0) {
            cerr << "[clip-vg]: Collapsed " << num_overlaps
                 << " intervals that -k extended into each other" << endl;
        }
    }

    if (progress) {
        size_t num_intervals = 0;
        for (auto& bi : bed_intervals) {
            num_intervals += bi.second.size();
        }
        cerr << "[clip-vg]: Clipping " << num_intervals << " 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, allow_ref_cycles, progress);
        if (forwardize_nonref) {
            forwardize_nonref_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, leave_aligned_drop_paths, progress);
    }

    if (!out_bed_path.empty()) {
        unordered_map<string, vector<pair<int64_t, int64_t>>> output_graph_intervals = get_path_intervals(graph.get());
#ifdef debug
        for (const auto& xx : output_graph_intervals) {
            cerr << " got output intervals " << xx.first << " count = " << xx.second.size() << endl;
        }
#endif
        unordered_map<string, vector<pair<int64_t, int64_t>>> clipped_graph_intervals = get_clipped_intervals(input_graph_intervals, output_graph_intervals);
        ofstream out_bed_file(out_bed_path);
        size_t icount = 0;
        for (const auto& pi : clipped_graph_intervals) {
            for (const auto& i : pi.second) {
                out_bed_file << pi.first << "\t" << i.first << "\t" << i.second << "\n";
                ++icount;
            }
        }
        out_bed_file.flush();
        if (progress) {
            cerr << "[clip-vg]: Outputted " << icount << " clipped intervals to " << out_bed_path << endl;
        }
    }

    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;
}

// anchor-prefix means we consider a node unaligned if it doesn't align to a path with that prefix,
// so we need a table of the nodes on those paths.  Built here rather than inside find_unaligned
// because -k needs it too, and -k can be given with -b, where find_unaligned never runs.
void build_anchor_nodes(const PathHandleGraph* graph, const string& anchor_prefix,
                        unordered_set<nid_t>& anchor_nodes_out) {
    if (anchor_prefix.empty()) {
        return;
    }
    graph->for_each_path_handle([&](path_handle_t path_handle) {
            string path_name = graph->get_path_name(path_handle);
            if (path_name.compare(0, anchor_prefix.length(), anchor_prefix) == 0) {
                graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
                    anchor_nodes_out.insert(graph->get_id(graph->get_handle_of_step(step_handle)));
                });
            }
        });
}

unordered_map<string, vector<pair<int64_t, int64_t>>> find_unaligned(const PathHandleGraph* graph, int64_t max_unaligned,
                                                                     const string& ref_prefix, const string& anchor_prefix,
                                                                     const unordered_set<nid_t>& anchor_nodes) {
    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 = step_is_aligned(graph, handle, path_handle, anchor_nodes,
                                                       !anchor_prefix.empty());
                        // start an unaligned interval
                        if (start < 0 && aligned == false) {
                            start = offset;
                        }
                        // end an unaligned interval
                        if (aligned == true) {
                            if (start >= 0 && offset - 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 {
        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;
    // paths to destroy (faster to do in single api call)
    vector<path_handle_t> paths_to_destroy;
    
    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 << "adding path to destroy list" << graph->get_path_name(path_handle) << endl;
#endif
                paths_to_destroy.push_back(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;
        }
    }

    // delete all the paths
#ifdef debug
    cerr << "destroying " << paths_to_destroy.size() << " paths" << endl;
#endif
    graph->destroy_paths(paths_to_destroy);
    paths_to_destroy.clear();

    // 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);
                    });
                paths_to_destroy.push_back(subpath_handle);
                if (progress) {
                    cerr << "[clip-vg]: Removing orphaned subpath " << graph->get_path_name(subpath_handle) << endl;
                }
                ++removed_subpath_count;
            }
        }
        graph->destroy_paths(paths_to_destroy);

        // 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));
        });
    vector<path_handle_t> paths_to_destroy;
    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));
                });
            paths_to_destroy.push_back(path_handle);
        }
    }
    graph->destroy_paths(paths_to_destroy);
    
    if (progress) {
        cerr << "[clip-vg]: Replaced " << replacement_count << " substrings in " << path_count << " path names" << endl;
    }
}

// Replace a node by its reverse complement, so that every step that visited it backwards
// now visits it forwards (and every step that visited it forwards now visits it backwards).
// The incident edges and every step are rewired onto a new node, and the old one destroyed,
// so the node comes out with a fresh id.
// Note that this rewrites every step on the node, not just the ones going backwards, so it must
// not be called while a path containing the node is being traversed.
void flip_node(MutablePathMutableHandleGraph* graph, nid_t node_id) {
    handle_t handle = graph->flip(graph->get_handle(node_id));
    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);
    for (step_handle_t step : steps) {
        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});
    }
    assert(graph->steps_of_handle(handle).empty());
    dynamic_cast<DeletableHandleGraph*>(graph)->destroy_handle(handle);
}

// bits of the per-node census flags below
static const uint8_t CENSUS_SEEN = 1;
static const uint8_t CENSUS_REV = 2;
static const uint8_t CENSUS_MULTI = 4;

void forwardize_paths(MutablePathMutableHandleGraph* graph, const string& ref_prefix, bool allow_ref_cycles, bool progress) {

    // Census state, reused across paths.  A byte per node id rather than a hash entry: the map this
    // replaced cost ~76 bytes per reference-path node and was paid on every run, including the
    // overwhelmingly common one where nothing needs flipping at all.  Only nodes the path visits more
    // than once -- reference cycles, which are rare -- need their visits counted exactly, so those go
    // in a map that stays tiny.
    if (graph->get_node_count() == 0) {
        return;
    }
    const nid_t base_id = graph->min_node_id();
    vector<uint8_t> flags;
    vector<nid_t> order;                               // first-visit order, for determinism
    unordered_map<nid_t, pair<size_t, size_t>> multi;  // cycled node -> (visits, backwards visits)

    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;
                order.clear();
                multi.clear();
                // flipping mints ids above the old maximum, so the range can grow between paths
                size_t needed = graph->max_node_id() - base_id + 1;
                if (flags.size() < needed) {
                    flags.resize(needed, 0);
                }

                // Take a census of the path before touching anything.  Flipping from inside
                // for_each_step_in_path() is not safe: the traversal caches the path's last step up
                // front and stops only when it reaches it, but flipping rewrites every step on the
                // node.  A reference cycle whose other visit happens to be that last step therefore
                // destroys the loop's terminator and walks it off the end of the path.
                graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
                        handle_t handle = graph->get_handle_of_step(step_handle);
                        nid_t node_id = graph->get_id(handle);
                        uint8_t& f = flags[node_id - base_id];
                        bool is_rev = graph->get_is_reverse(handle);
                        if (f & CENSUS_SEEN) {
                            if (!(f & CENSUS_MULTI)) {
                                // second visit: start counting, seeding with the first
                                f |= CENSUS_MULTI;
                                multi[node_id] = make_pair((size_t)1, (size_t)((f & CENSUS_REV) ? 1 : 0));
                            }
                            pair<size_t, size_t>& visits = multi[node_id];
                            ++visits.first;
                            if (is_rev) {
                                ++visits.second;
                            }
                        } else {
                            f |= CENSUS_SEEN;
                            order.push_back(node_id);
                        }
                        if (is_rev) {
                            f |= CENSUS_REV;
                        }
                        ++total_steps;
                    });

                for (nid_t node_id : order) {
                    uint8_t f = flags[node_id - base_id];
                    flags[node_id - base_id] = 0;  // reset as we go, so the next path starts clean
                    if (!(f & CENSUS_REV)) {
                        // every visit is already forward
                        continue;
                    }
                    // visited once unless it is one of the rare cycled nodes
                    size_t ref_count = 1;
                    size_t rev_count = 1;
                    if (f & CENSUS_MULTI) {
                        const pair<size_t, size_t>& visits = multi[node_id];
                        ref_count = visits.first;
                        rev_count = visits.second;
                    }
                    if (ref_count > 1) {
                        if (!allow_ref_cycles) {
                            cerr << "[clip-vg] error: Cycle detected in reference path " << path_name << " at node " << node_id << endl;
                            exit(1);
                        }
                        if (rev_count != ref_count) {
                            // The reference walks this node both ways round, so no orientation
                            // satisfies every visit -- flipping would just move the reversal onto
                            // the other step.  Leave it be.
                            continue;
                        }
                        // Every visit is backwards, so one flip forwardizes all of them.
                    }
                    flip_node(graph, node_id);
                    fw_count += rev_count;
                }
                if (fw_count > 0 && progress) {
                    cerr << "[clip-vg]: Forwardized " << fw_count << " / " << total_steps << " steps in reference path " << path_name << endl;
                }
            }
        });

    if (!allow_ref_cycles) {
        // do a check just to be sure.  Not valid with -c, where a node the reference walks both
        // ways round is deliberately left reverse above.
        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);
                            }
                        });
                }
            });
    }
}

void forwardize_nonref_paths(MutablePathMutableHandleGraph* graph, const string& ref_prefix, bool progress) {

    if (graph->get_node_count() == 0) {
        return;
    }

    // Nodes on a reference path are off limits: forwardize_paths() has just made every
    // reference step forward, and flipping one of these would undo that.  A bit per node id
    // rather than a hash set -- this is one entry per reference base at whole-genome scale.
    const nid_t base_id = graph->min_node_id();
    vector<bool> ref_nodes(graph->max_node_id() - base_id + 1, false);
    size_t ref_path_count = 0;
    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) {
                ++ref_path_count;
                graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
                        ref_nodes[graph->get_id(graph->get_handle_of_step(step_handle)) - base_id] = true;
                    });
            }
        });

    // With no reference there is nothing for the rest of the graph to be forward with
    // respect to, and worse, every node would look eligible below -- a mistyped -e would
    // quietly rewrite the whole graph instead of doing nothing.  Say so and leave it alone.
    if (ref_path_count == 0) {
        cerr << "[clip-vg]: warning: no path name begins with \"" << ref_prefix
             << "\", so nothing was forwardized" << endl;
        return;
    }

    // A node that every path walks backwards can be flipped with nothing to weigh up:
    // afterwards every step on it is forward and no step has been made reverse.  It is
    // the whole-run case that matters here -- a contig aligned in reverse to the
    // reference gives a stretch of such nodes, and nothing downstream can flip them once
    // the graph is written, because doing so is a change of topology rather than of
    // annotation.  Nodes that some path walks forwards are left alone: flipping one of
    // those only moves the reversal onto a different path.
    //
    // Each node's decision depends only on its own steps, so it does not matter what
    // order they are visited in.  The flips themselves are applied in node id order so
    // that the ids handed out to the replacement nodes are deterministic too.
    vector<nid_t> to_flip;
    graph->for_each_handle([&](handle_t handle) {
            nid_t node_id = graph->get_id(handle);
            if (ref_nodes[node_id - base_id]) {
                return;
            }
            bool has_step = false;
            bool all_reverse = true;
            graph->for_each_step_on_handle(handle, [&](step_handle_t step_handle) {
                    has_step = true;
                    if (!graph->get_is_reverse(graph->get_handle_of_step(step_handle))) {
                        all_reverse = false;
                    }
                    return all_reverse;
                });
            if (has_step && all_reverse) {
                to_flip.push_back(node_id);
            }
        });
    std::sort(to_flip.begin(), to_flip.end());

    size_t flipped_bases = 0;
    for (nid_t node_id : to_flip) {
        flipped_bases += graph->get_length(graph->get_handle(node_id));
        flip_node(graph, node_id);
    }

    if (!to_flip.empty() && progress) {
        cerr << "[clip-vg]: Forwardized " << to_flip.size() << " nodes (" << flipped_bases
             << " bp) that no path visited forward" << endl;
    }

    // No verification sweep here on purpose.  It would be a second full pass over every handle and
    // every step -- 40% of this tool's runtime with -e -- and it cannot fire: to_flip was built from
    // a complete sweep with exactly the predicate it would re-test, flip_node() leaves every step on
    // the replacement node forward, and the ids it mints are never in ref_nodes.
}

// 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 leave_aligned, 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
            size_t offset = 0;
            vector<pair<int64_t, int64_t>> intervals;
            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;
                    size_t len = graph->get_length(handle);
                    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));
                        if (!intervals.empty() && offset == intervals.back().second) {
                            intervals.back().second += len;
                        } else {
                            intervals.push_back(make_pair(offset, offset + len));
                        }
                    }
                    offset += len;
                });
            if (leave_aligned && !intervals.empty()) {
                // chop the path (to keep fragments for any non-destroyed nodes)
                chop_path(graph, path_handle, intervals);
            }
            if (!leave_aligned || !intervals.empty()) {
                // 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;        
    }
}
// Extend each clipped interval outward along its path for as long as unaligned sequence stays
// dense, up to max_flank bp on each side.
//
// An aligner that extends an anchor into a repeat leaves a fringe: sequence that never reaches the
// anchor graph, in runs too short for -u to clip.  So we score each base +(1-threshold) when it
// is unaligned by the same test that cut the intervals, and -threshold otherwise, walk outward
// and stop at the furthest point where the running total was at its maximum.  That is Mott's
// trimming rule: the cut still lands past a locally aligned stretch if unaligned sequence resumes
// beyond it, but a sustained aligned stretch ends the trim.
//
// Two consequences of that rule, both intended, neither obvious.  Aligned sequence inside a walk's
// reach is removed -- crossing a local dip is the whole point, so an anchored stretch with fringe
// on the far side of it goes.  And an island between two clipped intervals is approached from both
// sides, so it survives only if it is wider than the two walks' reach combined, not one walk's.
// When they do meet, every base between them was claimed by one walk or the other, so collapsing
// the overlap downstream deletes nothing that was not already scored for deletion.
void extend_flanks(const PathHandleGraph* graph,
                   unordered_map<string, vector<pair<int64_t, int64_t>>>& intervals,
                   int64_t max_flank, double threshold,
                   const unordered_set<nid_t>& anchor_nodes, bool have_anchor_prefix,
                   bool progress) {
    int64_t total_added = 0;
    size_t num_extended = 0;

    // flatten a path into the length of the node at each step and the offset it starts at
    auto flatten = [&](path_handle_t path_handle, vector<int64_t>& lengths, vector<int64_t>& offsets,
                       vector<bool>& unaligned) {
        lengths.clear();
        offsets.clear();
        unaligned.clear();
        int64_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);
                int64_t len = graph->get_length(handle);
                offsets.push_back(offset);
                lengths.push_back(len);
                unaligned.push_back(!step_is_aligned(graph, handle, path_handle, anchor_nodes,
                                                     have_anchor_prefix));
                offset += len;
            });
        return offset;
    };

    if (threshold < 0) {
        // Calibrate against this graph rather than trusting a fixed number.  Over the 24 human
        // chromosomes of a 460-haplotype pangenome the near-field fraction ran 0.10 to 0.45 and
        // the threshold it implies 0.054 to 0.230 -- a 4.3x spread tracking neither haplotype
        // count nor chromosome size, with chr8 lowest, chr13 highest, and a 118-haplotype chrY
        // near the bottom.  A constant that gates one of those sensibly fails to gate another at
        // all.  Compare sequence just outside the clipped intervals against sequence far from any
        // of them, and put the threshold half way between.
        //
        // The two distances only have to bracket the fringe.  Its density decays over tens of kb,
        // so 5kb is inside it on every graph measured, and 200kb is outside: the far-field
        // fraction came back between 0.0036 and 0.0128 across those same 24 chromosomes, which is
        // background even alongside the largest arrays.
        const int64_t CALIB_NEAR = 5000;
        const int64_t CALIB_FAR = 200000;
        // Both samples need to be worth believing, but they are estimating opposite things, so
        // the floors are on different quantities.  far is establishing that a rate is *low*, which
        // takes length: at the 0.0036-0.0128 measured above, 100kb carries a few hundred unaligned
        // bases.  near is pinning down a rate that is large, so what it needs is those bases
        // themselves -- 500 of them, the same order 100kb of background yields, which at a 10-45%
        // fringe is a few kb of band rather than 100.  Demanding 100kb of near as well would
        // refuse any graph with a single clipped interval, where 10kb of band already carries a
        // perfectly good sample.  They arrive in whole nodes rather than independently, so these
        // are floors on plausibility, not confidence intervals.
        const int64_t CALIB_MIN_FAR = 100000;
        const int64_t CALIB_MIN_NEAR_UNALIGNED = 500;
        int64_t near_tot = 0, near_unaligned = 0, far_tot = 0, far_unaligned = 0;
        vector<int64_t> lengths, offsets;
        vector<bool> unaligned;

        for (auto& pi : intervals) {
            if (!graph->has_path(pi.first)) {
                continue;
            }
            flatten(graph->get_path_handle(pi.first), lengths, offsets, unaligned);
            if (lengths.empty()) {
                continue;
            }
            vector<int64_t> edges;
            for (auto& interval : pi.second) {
                edges.push_back(interval.first);
                edges.push_back(interval.second);
            }
            sort(edges.begin(), edges.end());

            for (size_t i = 0; i < lengths.size(); ++i) {
                int64_t mid = offsets[i] + lengths[i] / 2;
                // The intervals are disjoint and sorted, so edges alternate start,end,start,end...
                // and a position sits inside one exactly when an odd number of edges precede it.
                // Skip those -- they are going anyway -- without rescanning every interval.
                size_t after = upper_bound(edges.begin(), edges.end(), mid) - edges.begin();
                if (after % 2 == 1) {
                    continue;
                }
                int64_t dist = numeric_limits<int64_t>::max();
                if (after > 0) {
                    dist = min(dist, mid - edges[after - 1]);
                }
                if (after < edges.size()) {
                    dist = min(dist, edges[after] - mid);
                }
                if (dist <= CALIB_NEAR) {
                    near_tot += lengths[i];
                    if (unaligned[i]) near_unaligned += lengths[i];
                } else if (dist >= CALIB_FAR) {
                    far_tot += lengths[i];
                    if (unaligned[i]) far_unaligned += lengths[i];
                }
            }
        }

        double near_frac = near_tot ? (double)near_unaligned / near_tot : 0.;
        double far_frac = far_tot ? (double)far_unaligned / far_tot : 0.;
        // An unsampled background is not a clean one.  Left ungated, far_tot == 0 gives far_frac 0
        // -- the most permissive answer there is -- and the threshold below collapses to
        // near_frac / 2, which on uniformly unaligned sequence sits under the background rate by
        // construction: the walk never turns negative and every trim runs to its -k cap over
        // ordinary sequence.  Nothing reaches the far bucket unless some path runs more than
        // CALIB_FAR past a clipped interval, which no short contig does, and a handful of far bases
        // estimates a fraction no better than none of them do.  Demand a real sample or decline.
        // Both samples get a floor, and the same one.  near_frac sets the threshold directly --
        // with the background this far below it, threshold is within a few percent of near/2 --
        // so an unreliable near estimate is at least as damaging as an unreliable far one.
        if (far_tot < CALIB_MIN_FAR) {
            cerr << "[clip-vg]: Flank calibration sampled " << far_tot << "bp beyond " << CALIB_FAR
                 << "bp of a clipped interval and wants " << CALIB_MIN_FAR << "bp: too little"
                 << " background to estimate an unaligned rate from.  Skipping flank trimming --"
                 << " pass -T explicitly to trim without calibrating." << endl;
            return;
        }
        if (near_unaligned < CALIB_MIN_NEAR_UNALIGNED) {
            cerr << "[clip-vg]: Flank calibration found " << near_unaligned << "bp of unaligned"
                 << " sequence within " << CALIB_NEAR << "bp of a clipped interval and wants "
                 << CALIB_MIN_NEAR_UNALIGNED << "bp: too little fringe to measure.  Skipping flank"
                 << " trimming -- pass -T explicitly to trim without calibrating." << endl;
            return;
        }
        // A fringe has to stand out from the background, not merely exceed it.  The guard below
        // only bounds the sign of the difference, and a graph that is uniformly a few percent
        // unaligned would pass it and then calibrate to a threshold under its own ambient rate --
        // at which the walk never turns negative and every trim runs to the -k cap over ordinary
        // sequence.  Measured separations are 21x to 97x across 24 human chromosomes, so 2x
        // refuses the pathology without coming close to real data.
        if (near_frac <= 2. * far_frac) {
            cerr << "[clip-vg]: Flank calibration found no fringe to speak of (near " << near_frac
                 << " vs far " << far_frac << ", less than 2x); skipping flank trimming" << endl;
            return;
        }
        threshold = far_frac + 0.5 * (near_frac - far_frac);
        // Unconditional, like the two refusals above: this number decides how much sequence -k
        // deletes, it came from the graph rather than the command line, and a run that does not
        // pass -p still needs it in its log to be reproducible.
        cerr << "[clip-vg]: Flank calibration: " << near_frac << " of bases unaligned within "
             << CALIB_NEAR << "bp of a clipped interval vs " << far_frac << " beyond "
             << CALIB_FAR << "bp; using threshold " << threshold << endl;
    }

    for (auto& pi : intervals) {
        if (!graph->has_path(pi.first)) {
            continue;
        }
        path_handle_t path_handle = graph->get_path_handle(pi.first);

        vector<int64_t> lengths;
        vector<int64_t> offsets;
        vector<bool> unaligned;
        int64_t path_length = flatten(path_handle, lengths, offsets, unaligned);
        if (lengths.empty()) {
            continue;
        }

        // The score of one node, in Mott's formulation.  eff is how much of the node counts: the
        // part of it that lies outside the interval being extended and outside the neighbouring
        // intervals.  That is the whole node everywhere except at the two ends of a walk.
        auto node_score = [&](size_t idx, int64_t eff) -> double {
            return (unaligned[idx] ? (double)eff : 0.0) - threshold * (double)eff;
        };

        vector<pair<int64_t, int64_t>>& path_intervals = pi.second;
        int64_t prev_end = 0;
        for (size_t j = 0; j < path_intervals.size(); ++j) {
            pair<int64_t, int64_t>& interval = path_intervals[j];
            // Neither walk may score sequence that is already being clipped.  The intervals are
            // sorted and disjoint here, and the inside of a neighbour is unaligned by construction
            // -- that is why it was cut -- so a walk let into one would always find a new maximum
            // there and drag the boundary across the surviving island in between.  Both bounds are
            // the neighbours' original edges, so what an interval extends to does not depend on
            // the order the intervals are visited in.
            int64_t left_bound = prev_end;
            int64_t right_bound = j + 1 < path_intervals.size() ? path_intervals[j + 1].first : path_length;
            // Clamped and monotonic so a bad record below cannot hand the next interval a bound
            // outside the path.  Identical to interval.second for anything well formed, since the
            // intervals are sorted and disjoint.
            prev_end = max(prev_end, min(interval.second, path_length));

            // -b intervals come from a file, and load_bed only rejects strict overlaps: a record
            // can be empty, inverted, or off the end of the path.  Both walks would repair such a
            // record into a plausible interval that chop_path then honours, so a typo becomes
            // deleted sequence -- and for a wholly off-path record it converts an abort into
            // silent data loss.  Leave anything we cannot reason about exactly as it arrived, and
            // let the code downstream treat it as it always has.
            if (interval.first < 0 || interval.second > path_length ||
                interval.first >= interval.second) {
                continue;
            }

            // Walk left over [left_bound, interval.first), nearest node first.  A -b interval can
            // begin inside a node -- the ones find_unaligned makes never do -- so the first node
            // scored is whichever holds interval.first - 1, and only its part outside counts.
            size_t start_idx = lower_bound(offsets.begin(), offsets.end(), interval.first) - offsets.begin();
            double running = 0., best = 0.;
            int64_t best_pos = interval.first;
            for (int64_t i = (int64_t)start_idx - 1; i >= 0; --i) {
                int64_t from = max(offsets[i], left_bound);
                int64_t to = min(offsets[i] + lengths[i], interval.first);
                if (to <= from || interval.first - from > max_flank) {
                    break;
                }
                running += node_score(i, to - from);
                if (running > best) {
                    best = running;
                    best_pos = from;
                }
            }
            if (best_pos < interval.first) {
                total_added += interval.first - best_pos;
                ++num_extended;
                interval.first = best_pos;
            }

            // and right over (interval.second, right_bound], starting likewise on the node that
            // holds interval.second when the interval does not end on a boundary
            size_t end_idx = lower_bound(offsets.begin(), offsets.end(), interval.second) - offsets.begin();
            if (end_idx > 0 && offsets[end_idx - 1] + lengths[end_idx - 1] > interval.second) {
                --end_idx;
            }
            running = 0.; best = 0.;
            best_pos = interval.second;
            for (size_t i = end_idx; i < lengths.size(); ++i) {
                int64_t from = max(offsets[i], interval.second);
                int64_t to = min(offsets[i] + lengths[i], right_bound);
                if (to <= from || to - interval.second > max_flank) {
                    break;
                }
                running += node_score(i, to - from);
                if (running > best) {
                    best = running;
                    best_pos = to;
                }
            }
            if (best_pos > interval.second) {
                total_added += best_pos - interval.second;
                ++num_extended;
                interval.second = best_pos;
            }
        }
    }

    if (progress) {
        cerr << "[clip-vg]: Extended " << num_extended << " interval ends by " << total_added
             << " bp of unaligned flank" << endl;
    }
}

unordered_map<string, vector<pair<int64_t, int64_t>>> get_path_intervals(const PathHandleGraph* graph) {
    unordered_map<string, vector<pair<int64_t, int64_t>>> path_intervals;
    graph->for_each_path_handle([&](path_handle_t path_handle) {
            string path_name = strip_subpath_name(graph->get_path_name(path_handle));
            size_t path_len = 0;
            graph->for_each_step_in_path(path_handle, [&](step_handle_t step_handle) {
                    path_len += graph->get_length(graph->get_handle_of_step(step_handle));
                });
            subrange_t path_range = graph->get_subrange(path_handle);
            // .first is where the subpath starts; .second is where it ends, and using that put
            // every subpath's interval a whole path-length to the right of the truth
            int64_t start = path_range == PathMetadata::NO_SUBRANGE ? 0 : (int64_t)path_range.first;
            vector<pair<int64_t, int64_t>>& intervals = path_intervals[path_name];        
            intervals.push_back(make_pair(start, start + path_len));
        });

    for (auto& pi : path_intervals) {
        std::sort(pi.second.begin(), pi.second.end(), [](const pair<int64_t, int64_t>& i1, const pair<int64_t, int64_t>& i2) {
                return i1.first < i2.first || (i1.first == i2.first && i1.second < i2.second);
            });
    }

    return path_intervals;
}

static unordered_map<string, vector<pair<int64_t, int64_t>>> get_clipped_intervals(
    const unordered_map<string, vector<pair<int64_t, int64_t>>>& input_intervals,
    const unordered_map<string, vector<pair<int64_t, int64_t>>>& output_intervals) {

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

    for (const auto& input_pi : input_intervals) {
        const string& path_name = input_pi.first;
        const auto& in_intervals = input_pi.second;
        if (!output_intervals.count(path_name)) {
            // path doesn't appear in output -> everything was clipped
            clipped_intervals[path_name].insert(clipped_intervals[path_name].end(), in_intervals.begin(), in_intervals.end());
#ifdef debug
            cerr << "clippin everything for " << path_name << endl;
#endif
        } else {
#ifdef debug
            cerr << "doin frag clip for " << path_name << endl;
#endif
            const auto& out_intervals = output_intervals.at(path_name);
            // note: clipping here is fairly simple because output intervals are a subset
            // and nothing overlaps
            int64_t j = 0;
            int64_t k = 0;
            for (int64_t i = 0; i < (int64_t)in_intervals.size(); ++i) {
                // j: first out_interval that's not completely left of in_interval;
                for (; j < (int64_t)out_intervals.size() && out_intervals[j].second < in_intervals[i].first; ++j);
                // k: first out_interval that's completely right of in_interval
                for (k = j; k < (int64_t)out_intervals.size() && out_intervals[k].first < in_intervals[i].second; ++k);
                // [j, k) all out_intervals that are sub intervals of in_interval[i]
                vector<pair<int64_t, int64_t>>&  gap_intervals = clipped_intervals[path_name];
                if (k == j) {
                    // nothing overlaps -- everything clipped
                    gap_intervals.push_back(in_intervals[i]);
                } else {
                    // left of first interval                    
                    if (out_intervals[j].first > in_intervals[i].first) {
                        gap_intervals.push_back(make_pair(in_intervals[i].first, out_intervals[j].first));
                    }
                    // gaps between the out_intervals that fall inside this in_interval.  Bounded
                    // by [j, k) rather than the whole vector: stopping at size - 2 dropped the last
                    // gap, and scanning all of them reported every gap once per in_interval.
                    for (int64_t l = j; l + 1 < k; ++l) {
                        if (out_intervals[l+1].first > out_intervals[l].second) {
                            gap_intervals.push_back(make_pair(out_intervals[l].second, out_intervals[l+1].first));
                        }
                    }
                    // right of the last interval inside this one
                    if (in_intervals[i].second > out_intervals[k-1].second) {
                        gap_intervals.push_back(make_pair(out_intervals[k-1].second, in_intervals[i].second));
                    }
                }
            }
        }
    }

    return clipped_intervals;
}