diff --git a/.github/workflows/build_and_test_on_push.yml b/.github/workflows/build_and_test_on_push.yml index 582a1b658..37e626a12 100644 --- a/.github/workflows/build_and_test_on_push.yml +++ b/.github/workflows/build_and_test_on_push.yml @@ -1,9 +1,9 @@ -on: [ push ] +on: [ push, pull_request ] name: build and test jobs: - build_and_test: + linux: runs-on: ubuntu-20.04 steps: - uses: actions/checkout@v2 @@ -27,3 +27,20 @@ jobs: run: ASAN_OPTIONS=detect_leaks=1:symbolize=1 LSAN_OPTIONS=verbosity=0:log_threads=1 bin/odgi test - name: Run remaining tests run: ctest --test-dir build -E odgi-test --verbose + + macos: + runs-on: macos-latest + steps: + - uses: actions/checkout@v2 + - name: Install required packages + run: brew update && brew install cmake llvm jemalloc + - name: Init and update submodules + run: git submodule update --init --recursive + - name: Build odgi + run: | + CC=$(brew --prefix llvm)/bin/clang CXX=$(brew --prefix llvm)/bin/clang++ LDFLAGS=-L/opt/homebrew/lib cmake . -Bbuild + CC=$(brew --prefix llvm)/bin/clang CXX=$(brew --prefix llvm)/bin/clang++ LDFLAGS=-L/opt/homebrew/lib cmake --build build -- -j + - name: Run odgi program tests + run: ASAN_OPTIONS=detect_leaks=1:symbolize=1 LSAN_OPTIONS=verbosity=0:log_threads=1 bin/odgi test + - name: Run remaining tests + run: ctest --test-dir build -E odgi-test --verbose diff --git a/.readthedocs.yaml b/.readthedocs.yaml new file mode 100644 index 000000000..31dbf0d70 --- /dev/null +++ b/.readthedocs.yaml @@ -0,0 +1,35 @@ +# Read the Docs configuration file for Sphinx projects +# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details + +# Required +version: 2 + +# Set the OS, Python version and other tools you might need +build: + os: ubuntu-22.04 + tools: + python: "3.12" + # You can also specify other tool versions: + # nodejs: "20" + # rust: "1.70" + # golang: "1.20" + +# Build documentation in the "docs/" directory with Sphinx +sphinx: + configuration: docs/conf.py + # You can configure Sphinx to use a different builder, for instance use the dirhtml builder for simpler URLs + # builder: "dirhtml" + # Fail on all warnings to avoid broken references + # fail_on_warning: true + +# Optionally build your docs in additional formats such as PDF and ePub +# formats: +# - pdf +# - epub + +# Optional but recommended, declare the Python requirements required +# to build your documentation +# See https://docs.readthedocs.io/en/stable/guides/reproducible-builds.html +python: + install: + - requirements: docs/requirements.txt diff --git a/README.md b/README.md index 945f6bd49..d8789239b 100644 --- a/README.md +++ b/README.md @@ -51,6 +51,14 @@ Static builds are unlikely to be supported on OSX, and require appropriate stati For more information on optimisations, debugging and GNU Guix builds, see [INSTALL.md](./INSTALL.md) and [CMakeLists.txt](./CMakeLists.txt). +### Nix build + +If you have `nix`, build and installation in your profile are as simple as: + +``` +nix-build && nix-env -i ./result +``` + #### Notes for distribution If you need to avoid machine-specific optimizations, use the `CMAKE_BUILD_TYPE=Generic` build type: diff --git a/default.nix b/default.nix new file mode 100644 index 000000000..756bd8a73 --- /dev/null +++ b/default.nix @@ -0,0 +1,5 @@ +{ pkgs ? import {} }: + +pkgs.callPackage ./odgi.nix { + inherit (pkgs) stdenv fetchFromGitHub cmake jemalloc pkg-config python3 gcc zlib htslib gsl; +} diff --git a/docs/conf.py b/docs/conf.py index 9950659bc..11f169916 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -20,13 +20,13 @@ # -- Project information ----------------------------------------------------- project = u'odgi' -copyright = '2020-2023, *Guarracino A., *Heumos S., Nahnsen S., Prins P., Garrison E. Revision v0.8.2-1fa78aa' +copyright = '2020-2024, *Guarracino A., *Heumos S., Nahnsen S., Prins P., Garrison E. Revision v0.8.4-a19163ea' author = u'*Andrea Guarracino, *Simon Heumos, Sven Nahnsen, Pjotr Prins, Erik Garrison' # The short X.Y version -version = 'v0.8.2' +version = 'v0.8.4' # The full version, including alpha/beta/rc tags -release = '1fa78aa' +release = 'a19163ea' # -- General configuration --------------------------------------------------- diff --git a/docs/img/DRB1-3123.du.png b/docs/img/DRB1-3123.du.png index 5810380a9..665eb8101 100644 Binary files a/docs/img/DRB1-3123.du.png and b/docs/img/DRB1-3123.du.png differ diff --git a/docs/img/DRB1-3123.png b/docs/img/DRB1-3123.png index 58bd8d6c2..48163fa7a 100644 Binary files a/docs/img/DRB1-3123.png and b/docs/img/DRB1-3123.png differ diff --git a/docs/img/DRB1-3123.z.png b/docs/img/DRB1-3123.z.png index d90d157b9..7702cb6d7 100644 Binary files a/docs/img/DRB1-3123.z.png and b/docs/img/DRB1-3123.z.png differ diff --git a/docs/img/DRB1-3123_sorted.U1000.png b/docs/img/DRB1-3123_sorted.U1000.png new file mode 100644 index 000000000..73efc3645 Binary files /dev/null and b/docs/img/DRB1-3123_sorted.U1000.png differ diff --git a/docs/img/DRB1-3123_sorted.j10000.png b/docs/img/DRB1-3123_sorted.j10000.png new file mode 100644 index 000000000..d3b2e39e8 Binary files /dev/null and b/docs/img/DRB1-3123_sorted.j10000.png differ diff --git a/docs/img/DRB1-3123_sorted.x2.png b/docs/img/DRB1-3123_sorted.x2.png new file mode 100644 index 000000000..d6be38793 Binary files /dev/null and b/docs/img/DRB1-3123_sorted.x2.png differ diff --git a/docs/img/DRB1-3123_sorting_layouting.png b/docs/img/DRB1-3123_sorting_layouting.png index f001bd45f..45c2192c6 100644 Binary files a/docs/img/DRB1-3123_sorting_layouting.png and b/docs/img/DRB1-3123_sorting_layouting.png differ diff --git a/docs/img/LPA.b.png b/docs/img/LPA.b.png index 1c6ce8bc3..a4ef11ea3 100644 Binary files a/docs/img/LPA.b.png and b/docs/img/LPA.b.png differ diff --git a/docs/img/LPA.bm.VNTRs.png b/docs/img/LPA.bm.VNTRs.png index 7549c7a4b..c2faa8932 100644 Binary files a/docs/img/LPA.bm.VNTRs.png and b/docs/img/LPA.bm.VNTRs.png differ diff --git a/docs/img/LPA.bm.png b/docs/img/LPA.bm.png index 8ad5cabb2..88ca91864 100644 Binary files a/docs/img/LPA.bm.png and b/docs/img/LPA.bm.png differ diff --git a/docs/index.rst b/docs/index.rst index 4615dbbae..1e7adf14a 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -79,7 +79,7 @@ Core Functionalities :target: rst/tutorials/extract_selected_loci.html .. |sorting_layouting| image:: img/DRB1-3123_sorting_layouting.png - :target: rst/tutorials/sorting_layouting.html + :target: rst/tutorials/sort_layout.html .. |navigating_and_annotating_graphs| image:: img/nav_welcome.png :target: rst/tutorials/navigating_and_annotating_graphs.html diff --git a/docs/requirements.txt b/docs/requirements.txt index 0da073719..47a606942 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,3 +1,4 @@ breathe m2r2 asciidoc +sphinx-rtd-theme diff --git a/docs/rst/commands/odgi_crush.rst b/docs/rst/commands/odgi_crush.rst index ef61a54cf..d70a9b8a4 100644 --- a/docs/rst/commands/odgi_crush.rst +++ b/docs/rst/commands/odgi_crush.rst @@ -15,7 +15,7 @@ SYNOPSIS DESCRIPTION =========== -Divide nodes into smaller pieces preserving node topology and order. +Replaces runs of Ns with single Ns (for example, ANNNT becomes ANT). OPTIONS ======= diff --git a/docs/rst/commands/odgi_sort.rst b/docs/rst/commands/odgi_sort.rst index 274664a90..5f6c8a17c 100644 --- a/docs/rst/commands/odgi_sort.rst +++ b/docs/rst/commands/odgi_sort.rst @@ -49,9 +49,8 @@ order: force-directed graph drawing algorithm minimizes the graph’s energy function or stress level. It applies stochastic gradient descent (SGD) to move a single pair of nodes at a time. The path index is - used to pick the terms to move stochastically. If ran with 1 thread - only, the resulting order of the graph is deterministic. The seed is - adjustable. + used to pick the terms to move stochastically. For more details about + the algorithm, please take a look at https://www.biorxiv.org/content/10.1101/2023.09.22.558964v2. Sorting the paths in a graph my refine the sorting process. For the users’ convenience, it is possible to specify a whole pipeline of sorts diff --git a/docs/rst/multiqc.rst b/docs/rst/multiqc.rst index 4ac7e04e4..ed17cdc1a 100644 --- a/docs/rst/multiqc.rst +++ b/docs/rst/multiqc.rst @@ -38,7 +38,7 @@ To see the full statistics in YAML format of the graph, execute: .. code-block:: bash - odgi stats -i DRB1-3123.gfa.og -m + odgi stats -i DRB1-3123.gfa.og -m -sgdl This prints the following YAML to stdout: @@ -89,7 +89,7 @@ Let's save the statistics this time: .. code-block:: bash - odgi stats -i DRB1-3123.gfa.og -m > DRB1-3123.gfa.og.stats.yaml + odgi stats -i DRB1-3123.gfa.og -m -sgdl > DRB1-3123.gfa.og.stats.yaml .. note:: @@ -167,7 +167,7 @@ Assuming, we have several graphs, of which we want to compare the statistics fro .. code-bock:: bash odgi build -g LPA.gfa -o LPA.gfa.og - odgi stats -i LPA.gfa.og -y > LPA.gfa.og.stats.yaml + odgi stats -i LPA.gfa.og -m -sgdl > LPA.gfa.og.stats.yaml odgi viz -i LPA.gfa.og -o LPA.gfa.og.viz_mqc.png odgi layout -i LPA.gfa.og -o LPA.gfa.og.lay odgi draw -i LPA.gfa.og -c LPA.gfa.og.lay -p LPA.gfa.og.lay.draw_mqc.png -w 10 -C diff --git a/docs/rst/quick_start.rst b/docs/rst/quick_start.rst index c61a171c3..c399191b9 100644 --- a/docs/rst/quick_start.rst +++ b/docs/rst/quick_start.rst @@ -18,12 +18,12 @@ version 1 (`GFAv1 `_) Build graph from GFA ---------------------------- -Assuming that your current working directory is the root of the ``odgi`` project, to construct an ``odgi`` file from a -``GFA`` file, execute: +To construct an ``odgi`` file from a ``GFA`` file, execute: .. code-block:: bash - odgi build -g test/DRB1-3123.gfa -o DRB1-3123.og + wget https://raw.githubusercontent.com/pangenome/odgi/master/test/DRB1-3123.gfa + odgi build -g DRB1-3123.gfa -o DRB1-3123.og The command creates a file called ``DRB1-3123.og``, which contains the input graph in ``odgi`` format. diff --git a/docs/rst/tutorials/exploratory_analysis.rst b/docs/rst/tutorials/exploratory_analysis.rst index 7db08d80b..35cd38b1d 100644 --- a/docs/rst/tutorials/exploratory_analysis.rst +++ b/docs/rst/tutorials/exploratory_analysis.rst @@ -76,7 +76,7 @@ Color with respect to the node position This is a linearized visualization, but the pangenome graphs are not linear when the embedded genomes present structural variation. However, a graph can be optimized for being better visualized in 1-Dimension by sorting its nodes properly -(see the :ref:`sorting-layouting` tutorial for more information). +(see the :ref:`sort-layout` tutorial for more information). To color the bars with respect to the node position in each path, execute: diff --git a/docs/rst/tutorials/sort_layout.rst b/docs/rst/tutorials/sort_layout.rst index 8cbf5dd61..517a3fd4f 100644 --- a/docs/rst/tutorials/sort_layout.rst +++ b/docs/rst/tutorials/sort_layout.rst @@ -1,4 +1,4 @@ -.. _sorting-layouting: +.. _sort-layout: ############### Sort and Layout @@ -16,6 +16,8 @@ a 1D and 2D layout to simplify these complex regions. This tutorial shows how to sort and visualize a graph in 1D. It explains how to generate a 2D layout of a graph, and how to take a look at the calculated layout using static and interactive tools. +For more details about the applied algorithm, please take a look at https://www.biorxiv.org/content/10.1101/2023.09.22.558964v2. + .. Pangenome graphs embed linear pangenomic sequences as paths in .. the graph, but to our knowledge, no algorithm takes into account this biological information in the sorting. Moreover, .. existing 2D layout methods struggle to deal with large graphs. ``odgi`` implements a new layout algorithm to simplify a pangenome @@ -39,12 +41,12 @@ to take a look at the calculated layout using static and interactive tools. Build the unsorted DRB1-3123 graph ---------------------------------- -Assuming that your current working directory is the root of the ``odgi`` project, to construct an ``odgi`` graph from the -``DRB1-3123`` dataset in ``GFA`` format, execute: +To construct an ``odgi`` graph from the ``DRB1-3123`` dataset in ``GFA`` format, execute: .. code-block:: bash - odgi build -g test/DRB1-3123_unsorted.gfa -o DRB1-3123_unsorted.og + wget https://raw.githubusercontent.com/pangenome/odgi/master/test/DRB1-3123_unsorted.gfa + odgi build -g DRB1-3123_unsorted.gfa -o DRB1-3123_unsorted.og The command creates a file called ``DRB1-3123_unsorted.og``, which contains the input graph in ``odgi`` format. This graph contains 12 ALT sequences of the `HLA-DRB1 gene `_ from the GRCh38 reference genome. @@ -129,6 +131,22 @@ nodes. .. note:: The PG-SGD is not deterministic, because of its `Hogwild! `_ approach. + For more details about the applied algorithm, please take a look at https://www.biorxiv.org/content/10.1101/2023.09.22.558964v2. + +.. note:: + The 1D PG-SGD implementation comes with a huge amount of tunable parameters. Based on our experience applying it to hundreds of graphs, the current + defaults usually work well for most graphs. However, if you feel the sorting did not work well enough, there are 2 key parameters one can tune: + + | **-G, --path-sgd-min-term-updates-paths**\ =\ *N*: The minimum number of terms to be + updated before a new path-guided + linear 1D SGD iteration with adjusted + learning rate eta starts, expressed as + a multiple of total path steps (default: 1.0). + | **-x, --path-sgd-iter-max**\ =\ *N*: The maximum number of iterations for path-guided linear 1D SGD model (default: 100). + + Increasing both can lead to a better sorted graph. For example, one can start optimizing with setting **-x, --path-sgd-iter-max**\ =\ *200*. + For more parameter details please take + a look at :ref:`odgi sort`. .. To reproduce the visualization below, the sorted graph can be found under ``test/DRB1-3123_sorted.og``. @@ -169,6 +187,73 @@ This prints to stdout: Compared to before, these metrics show that the goodness of the sorting of the graph improved significantly. +-------------------------------------------- +Playing around with the 1D PG-SGD parameters +-------------------------------------------- + +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +What happens if the maximum number of iterations is very low? +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: bash + + odgi sort -i DRB1-3123_unsorted.og --threads 2 -P -Y -x 2 -o DRB1-3123_sorted.x2.og + odgi viz -i DRB1-3123_sorted.x2.og -o DRB1-3123_sorted.x2.png + +.. image:: /img/DRB1-3123_sorted.x2.png + +The graph appears very complex and not quite human readable. That's because in total there were two times the number +of total path steps node position updates instead of one hundred times the number of total path steps, which is the current default. +For very complex graphs, one may have to increase this number even further. + +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +What happens if the minimum number of term updates is very high? +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: bash + + odgi sort -i DRB1-3123_unsorted.og --threads 2 -P -Y -U 1000 -o DRB1-3123_sorted.U1000.og + odgi viz -i DRB1-3123_sorted.U1000.og -o DRB1-3123_sorted.U1000.png + +.. image:: /img/DRB1-3123_sorted.U1000.png + +The graph lost it's complexity and is now linear. Compared to the 1D visualization using the default parameters, it is hard +to spot any differences. So let's take a look at the metrics: + +.. code-block:: bash + + odgi stats -i DRB1-3123_sorted.U1000.og -s -d -l -g + +This prints to stdout: + +.. code-block:: bash + + #mean_links_length + path in_node_space in_nucleotide_space num_links_considered num_gap_links_not_penalized + all_paths 1.00361 8.30677 21870 15195 + #sum_of_path_node_distances + path in_node_space in_nucleotide_space nodes nucleotides num_penalties num_penalties_different_orientation + all_paths 3.23238 3.73489 21882 163416 3750 1 + +We actually were able to improve the metrics compared to using default parameters. However, the runtime increased from under 1 second to ~30 seconds. +So one needs to be careful with such a parameter. Compared to the gains in linearity, such an additional time usage would be a huge +waste with very large graphs. + +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ +What happens if the threshold of the maximum distance of two nodes is very high? +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. code-block:: bash + + odgi sort -i DRB1-3123_unsorted.og --threads 2 -P -Y -j 10000 -o DRB1-3123_sorted.j10000.og + odgi viz -i DRB1-3123_sorted.j10000.og -o DRB1-3123_sorted.j10000.png + +.. image:: /img/DRB1-3123_sorted.j10000.png + +The graph appears very complex and not quite human readable. That's because the iterations are terminated as soon as the +expected distance of two nodes, the nucleotide distance given by two randomly chosen path steps, is as close as 10000. +Naturally, this happens very soon. + ========================================================= 1D reference-guided grooming and reference-guided sorting ========================================================= @@ -267,6 +352,8 @@ We can clearly observe, that the path positions of the two reference now define 2D layout ========= +The 2D PG-SGD layout algorithm is described in https://www.biorxiv.org/content/10.1101/2023.09.22.558964v2. + ----------------------------------------- 2D layout of the unsorted DRB1-3123 graph ----------------------------------------- @@ -277,6 +364,23 @@ We want to have a 2D layout of our DRB1-3123 graph: odgi layout -i DRB1-3123_unsorted.og -o DRB1-3123_unsorted.og.lay -P --threads 2 +.. note:: + The 2D PG-SGD implementation comes with a huge amount of tunable parameters. Based on our experience applying it to hundreds of graphs, the current + defaults usually work well for most graphs. However, if you feel the resulting 2D layout is not of a good enough quality, there are 2 key parameters one can tune: + + | **-G, --path-sgd-min-term-updates-paths**\ =\ *N*: Minimum number of terms N to be + updated before a new path-guided 2D + SGD iteration with adjusted learning + rate eta starts, expressed as a + multiple of total path length + (default: 10). + | **-x, --path-sgd-iter-max**\ =\ *N*: The maximum number of iterations N for + the path-guided 2D SGD model (default: + 30). + + Increasing both can lead to a better graph layout. For example, one can start optimizing with setting **-x, --path-sgd-iter-max**\ =\ *100*. + For more parameter details please take a look at :ref:`odgi layout`. + -------------------------------------------- Drawing the 2D layout of the DRB1-3123 graph -------------------------------------------- diff --git a/odgi.nix b/odgi.nix new file mode 100644 index 000000000..e13b16e59 --- /dev/null +++ b/odgi.nix @@ -0,0 +1,40 @@ +{ lib, stdenv, fetchFromGitHub, cmake, jemalloc, python3, pkg-config, gcc, git, zlib, htslib, gsl }: + +stdenv.mkDerivation rec { + pname = "odgi"; + version = "0.8.3"; + + src = fetchFromGitHub { + owner = "pangenome"; + repo = "odgi"; + rev = "86e62bac42d737808a83bb00a0e7a60069494dcf"; + sha256 = "sha256-IKHQyP3E01LZvui6ykUiWPr2zuD2iqq55ARG+O2KCxM="; + fetchSubmodules = true; + }; + + nativeBuildInputs = [ cmake pkg-config ]; + + buildInputs = [ + jemalloc + gcc + zlib + htslib + gsl + python3 + ]; + + postPatch = '' + mkdir -p include + echo "#define ODGI_GIT_VERSION \"${version}\"" > include/odgi_git_version.hpp + ''; + + makeFlags = [ "CC=${gcc}/bin/gcc" ]; + + meta = with lib; { + description = "odgi optimized dynamic sequence graph implementation"; + homepage = "https://github.com/pangenome/odgi"; + license = licenses.mit; + platforms = platforms.linux; + maintainers = [ maintainers.yourNameHere ]; # Replace with your name + }; +} diff --git a/scripts/heaps_fit.R b/scripts/heaps_fit.R index f77b8a8a6..6e2338acd 100644 --- a/scripts/heaps_fit.R +++ b/scripts/heaps_fit.R @@ -32,5 +32,5 @@ print(z * (f(n) - f(n-1))) print(z * (f(2) - f(1))) #print(f(n) - f(n-1)) pdf(NULL) -ggplot(x, aes(x=nth.genome, y=base.pairs/1e9)) + geom_point(alpha=I(1/10)) + stat_function(fun=function(x) (fit$par[1] * x^fit$par[2] + fit$par[3]) * m) + scale_y_continuous("observed pangenome size (Gbp)") + scale_x_continuous(paste("Nth included genome (", max(x$permutation)+1 ," permutations) with gamma=", round(fit$par[2], digits=3), sep = "")) +ggplot(x, aes(x=nth.genome, y=base.pairs/1e9)) + geom_point(alpha=I(1/10)) + stat_function(fun=function(x) (fit$par[1] * x^fit$par[2] + fit$par[3]) * m) + scale_y_continuous("observed pangenome size (Gbp)") + scale_x_continuous(paste("Nth included genome (", max(x$permutation)+1 ," permutations) with gamma=", round(fit$par[2], digits=3), sep = "")) + expand_limits(x = 0, y = 0) ggsave(args[2], height=5, width=9) diff --git a/src/algorithms/atomic_image.cpp b/src/algorithms/atomic_image.cpp index 7f0a836b8..870e2759e 100644 --- a/src/algorithms/atomic_image.cpp +++ b/src/algorithms/atomic_image.cpp @@ -84,6 +84,25 @@ color_t mix(const color_t& a, const color_t& b, const double& f) { return out; } +std::string to_hex(const color_t& c) { + std::stringstream ss; + ss << "#"; + ss << std::hex << std::setfill('0') << std::setw(2) << (int)c.c.r; + ss << std::hex << std::setfill('0') << std::setw(2) << (int)c.c.g; + ss << std::hex << std::setfill('0') << std::setw(2) << (int)c.c.b; + return ss.str(); +} + +std::string to_rgba(const color_t& c) { + std::stringstream ss; + ss << "rgba("; + ss << (int)c.c.r << ","; + ss << (int)c.c.g << ","; + ss << (int)c.c.b << ","; + ss << (int)c.c.a << ")"; + return ss.str(); +} + // helpers double u_ipart(double x) { return std::floor(x); } diff --git a/src/algorithms/atomic_image.hpp b/src/algorithms/atomic_image.hpp index a9b380397..05472902f 100644 --- a/src/algorithms/atomic_image.hpp +++ b/src/algorithms/atomic_image.hpp @@ -10,6 +10,7 @@ #include #include #include +#include #include "picosha2.h" namespace odgi { @@ -92,6 +93,9 @@ color_t darkest(const color_t& a, const color_t& b); color_t layer(const color_t& a, const color_t& b, const double& f); color_t mix(const color_t& a, const color_t& b, const double& f); +std::string to_hex(const color_t& c); +std::string to_rgba(const color_t& c); + const color_t COLOR_BLACK = { 0xff000000 }; const color_t COLOR_LIGHTGRAY = { 0xffD3D3D3 }; const color_t COLOR_WHITE = { 0xffffffff }; diff --git a/src/algorithms/draw.cpp b/src/algorithms/draw.cpp index 7d86d3b75..7d606c8f2 100644 --- a/src/algorithms/draw.cpp +++ b/src/algorithms/draw.cpp @@ -104,7 +104,9 @@ void draw_svg(std::ostream &out, const std::vector &Y, const PathHandleGraph &graph, const double& scale, - const double& border) { + const double& border, + const double& line_width, + std::vector& node_id_to_color) { std::vector> weak_components; coord_range_2d_t rendered_range; @@ -124,9 +126,10 @@ void draw_svg(std::ostream &out, << "viewBox=\"" << viewbox_x1 << " " << viewbox_y1 << " " << width << " " << height << "\"" << " xmlns=\"http://www.w3.org/2000/svg\">" - << "" + // interferes with the styling of the lines + //<< "" << std::endl; auto range_itr = component_ranges.begin(); @@ -134,8 +137,10 @@ void draw_svg(std::ostream &out, auto& range = *range_itr++; auto& x_off = range.x_offset; auto& y_off = range.y_offset; + //const algorithms::color_t node_color = !node_id_to_color.empty() ? node_id_to_color[graph.get_id(handle)] : COLOR_BLACK; for (auto& handle : component) { uint64_t a = 2 * number_bool_packing::unpack_number(handle); + algorithms::color_t color = node_id_to_color.empty() ? COLOR_BLACK : node_id_to_color[graph.get_id(handle)]; out << "" << std::endl; diff --git a/src/algorithms/draw.hpp b/src/algorithms/draw.hpp index 8a1f1be50..6315e90f9 100644 --- a/src/algorithms/draw.hpp +++ b/src/algorithms/draw.hpp @@ -69,7 +69,9 @@ void draw_svg(std::ostream &out, const std::vector &Y, const PathHandleGraph &graph, const double& scale, - const double& border); + const double& border, + const double& line_width, + std::vector& node_id_to_color); std::vector rasterize(const std::vector &X, const std::vector &Y, diff --git a/src/algorithms/subgraph/extract.cpp b/src/algorithms/subgraph/extract.cpp index c570068f8..2818f6aef 100644 --- a/src/algorithms/subgraph/extract.cpp +++ b/src/algorithms/subgraph/extract.cpp @@ -53,9 +53,10 @@ namespace odgi { } path_handle_t create_subpath(graph_t &subgraph, const string &subpath_name, const bool is_circular) { - if (subgraph.has_path(subpath_name)) { - subgraph.destroy_path(subgraph.get_path_handle(subpath_name)); - } + // The function assumes that every path is new and unique + // if (subgraph.has_path(subpath_name)) { + // subgraph.destroy_path(subgraph.get_path_handle(subpath_name)); + // } return subgraph.create_path_handle(subpath_name, is_circular); }; diff --git a/src/subcommand/crush_main.cpp b/src/subcommand/crush_main.cpp index cd61a31a2..5f3149f59 100644 --- a/src/subcommand/crush_main.cpp +++ b/src/subcommand/crush_main.cpp @@ -19,7 +19,7 @@ int main_crush(int argc, char **argv) { argv[0] = (char *) prog_name.c_str(); --argc; - args::ArgumentParser parser("Divide nodes into smaller pieces preserving node topology and order."); + args::ArgumentParser parser("Replaces runs of Ns with single Ns (for example, ANNNT becomes ANT)."); args::Group mandatory_opts(parser, "[ MANDATORY ARGUMENTS ]"); args::ValueFlag og_in_file(mandatory_opts, "FILE", "Load the succinct variation graph in ODGI format from this *FILE*. The file name usually ends with *.og*. It also accepts GFAv1, but the on-the-fly conversion to the ODGI format requires additional time!", {'i', "idx"}); args::ValueFlag og_out_file(mandatory_opts, "FILE", "Write the N-crushed succinct variation graph in ODGI format to *FILE*. A file ending of *.og* is recommended.", diff --git a/src/subcommand/depth_main.cpp b/src/subcommand/depth_main.cpp index dbe33e62b..558175a7e 100644 --- a/src/subcommand/depth_main.cpp +++ b/src/subcommand/depth_main.cpp @@ -80,6 +80,10 @@ namespace odgi { "merging regions not separated by more than LEN bp." " When TIPS=1, retain only tips.", {'W', "windows-out"}); + args::Flag window_unique_depth(depth_opts, "window-unique-depth", + "For --window-in and --window-out, count UNIQUE depth, not total node depth", + {'U', "window-unique-depth"}); + args::Group threading_opts(parser, "[ Threading ] "); args::ValueFlag _num_threads(threading_opts, "N", "Number of threads to use in parallel operations.", {'t', "threads"}); @@ -437,7 +441,11 @@ namespace odgi { graph.for_each_handle( [&](const handle_t& h) { auto id = graph.get_id(h); - depths[id - shift] = get_graph_node_depth(graph, id, paths_to_consider).first; + if (window_unique_depth) { + depths[id - shift] = get_graph_node_depth(graph, id, paths_to_consider).second; // Unique depth + } else { + depths[id - shift] = get_graph_node_depth(graph, id, paths_to_consider).first; // Total depth + } }, true); auto in_bounds = diff --git a/src/subcommand/draw_main.cpp b/src/subcommand/draw_main.cpp index 12edf2e60..620e69d21 100644 --- a/src/subcommand/draw_main.cpp +++ b/src/subcommand/draw_main.cpp @@ -41,7 +41,7 @@ int main_draw(int argc, char **argv) { args::ValueFlag png_height(visualizations_opts, "FILE", "Height of PNG rendering (default: 1000).", {'H', "png-height"}); args::ValueFlag png_border(visualizations_opts, "FILE", "Size of PNG border in bp (default: 10).", {'E', "png-border"}); args::Flag color_paths(visualizations_opts, "color-paths", "Color paths (in PNG output).", {'C', "color-paths"}); - args::ValueFlag render_scale(visualizations_opts, "N", "Image scaling (default 1.0).", {'R', "scale"}); + args::ValueFlag render_scale(visualizations_opts, "N", "Image scaling (default 0.001).", {'R', "scale"}); args::ValueFlag render_border(visualizations_opts, "N", "Image border (in approximate bp) (default 100.0).", {'B', "border"}); args::ValueFlag png_line_width(visualizations_opts, "N", "Line width (in approximate bp) (default 0.0).", {'w', "line-width"}); //args::ValueFlag png_line_overlay(parser, "N", "line width (in approximate bp) (default 10.0)", {'O', "line-overlay"}); @@ -174,7 +174,7 @@ int main_draw(int argc, char **argv) { const double _png_line_width = png_line_width ? args::get(png_line_width) : 0; const bool _color_paths = args::get(color_paths); const double _png_path_line_spacing = png_path_line_spacing ? args::get(png_path_line_spacing) : 0.0; - const double svg_scale = !render_scale ? 1.0 : args::get(render_scale); + const double svg_scale = !render_scale ? 0.01 : args::get(render_scale); size_t max_node_depth = 0; graph.for_each_handle( [&](const handle_t& h) { @@ -216,7 +216,7 @@ int main_draw(int argc, char **argv) { // todo could be done with callbacks std::vector X = layout.get_X(); std::vector Y = layout.get_Y(); - algorithms::draw_svg(f, X, Y, graph, svg_scale, border_bp); + algorithms::draw_svg(f, X, Y, graph, svg_scale, border_bp, _png_line_width, node_id_to_color); f.close(); } diff --git a/src/subcommand/extract_main.cpp b/src/subcommand/extract_main.cpp index e1737eb46..ce0ecc022 100644 --- a/src/subcommand/extract_main.cpp +++ b/src/subcommand/extract_main.cpp @@ -31,16 +31,16 @@ namespace odgi { args::Group graph_files_io_opts(parser, "[ Graph Files IO ]"); args::ValueFlag og_out_file(graph_files_io_opts, "FILE", "Store all subgraphs in this FILE. The file name usually ends with *.og*.", {'o', "out"}); - args::ValueFlag _max_dist_subpaths(mandatory_opts, "N", + args::Group extract_opts(parser, "[ Extract Options ]"); + args::ValueFlag _max_dist_subpaths(extract_opts, "N", "Maximum distance between subpaths allowed for merging them. " "It reduces the fragmentation of unspecified paths in the input path ranges. " - "Set 0 to disable it.", + "Set 0 to disable it [default: 300000].", {'d', "max-distance-subpaths"}); - args::ValueFlag _num_iterations(mandatory_opts, "N", + args::ValueFlag _num_iterations(extract_opts, "N", "Maximum number of iterations in attempting to merge close subpaths. " - "It stops early if during an iteration no subpaths were merged [default: 3].", + "It stops early if during an iteration no subpaths were merged [default: 6].", {'e', "max-merging-iterations"}); - args::Group extract_opts(parser, "[ Extract Options ]"); args::Flag _split_subgraphs(extract_opts, "split_subgraphs", "Instead of writing the target subgraphs into a single graph, " "write one subgraph per given target to a separate file named path:start-end.og " @@ -84,6 +84,8 @@ namespace odgi { "List of paths to fully retain in the extracted graph. Must " "contain one path name per line and a subset of all paths can be specified.", {'R', "lace-paths"}); + args::Flag _optimize(extract_opts, "optimize", "Compact the node ID space in the extracted graph(s).", + {'O', "optimize"}); args::Group threading_opts(parser, "[ Threading ]"); args::ValueFlag nthreads(threading_opts, "N", "Number of threads to use for parallel operations.", {'t', "threads"}); @@ -116,11 +118,7 @@ namespace odgi { return 1; } - if (!_max_dist_subpaths) { - std::cerr << "[odgi::extract] error: please specify -d/--max-distance-subpaths. Values equal to or greater than 0 are allowed." << std::endl; - return 1; - } - if ((!_max_dist_subpaths || args::get(_max_dist_subpaths) == 0) && _num_iterations) { + if ((_max_dist_subpaths && args::get(_max_dist_subpaths) == 0) && _num_iterations) { std::cerr << "[odgi::extract] error: specified -e/--max-merging-iterations without specifying -d/--max-distance-subpaths greater than 0." << std::endl; return 1; } @@ -136,7 +134,8 @@ namespace odgi { return 1; } - const uint64_t num_iterations = _num_iterations && args::get(_num_iterations) > 0 ? args::get(_num_iterations) : 3; + const uint64_t max_dist_subpaths = _max_dist_subpaths && args::get(_max_dist_subpaths) >= 0 ? args::get(_max_dist_subpaths) : 300000; + const uint64_t num_iterations = _num_iterations && args::get(_num_iterations) > 0 ? args::get(_num_iterations) : 6; if (_split_subgraphs) { if (og_out_file) { @@ -268,25 +267,46 @@ namespace odgi { std::vector input_path_ranges; - // handle targets from BED - if (_path_bed_file && !args::get(_path_bed_file).empty()) { - std::ifstream bed_in(args::get(_path_bed_file)); - std::string line; - while (std::getline(bed_in, line)) { - add_bed_range(input_path_ranges, graph, line); + { + // handle targets from BED + if (_path_bed_file && !args::get(_path_bed_file).empty()) { + std::ifstream bed_in(args::get(_path_bed_file)); + std::string line; + while (std::getline(bed_in, line)) { + add_bed_range(input_path_ranges, graph, line); + } } - } - // handle targets from command line - if (_path_range) { - Region region; - parse_region(args::get(_path_range), region); + // handle targets from command line + if (_path_range) { + Region region; + parse_region(args::get(_path_range), region); - // no coordinates given, we do whole thing (0,-1) - if (region.start < 0 || region.end < 0) { - add_bed_range(input_path_ranges, graph, region.seq); - } else { - add_bed_range(input_path_ranges, graph, region.seq + "\t" + std::to_string(region.start) + "\t" + std::to_string(region.end)); + // no coordinates given, we do whole thing (0,-1) + if (region.start < 0 || region.end < 0) { + add_bed_range(input_path_ranges, graph, region.seq); + } else { + add_bed_range(input_path_ranges, graph, region.seq + "\t" + std::to_string(region.start) + "\t" + std::to_string(region.end)); + } + } + + // Check duplicates + std::vector copy_ranges = input_path_ranges; // Create a copy of the vector to avoid sorting the original one + + auto compare_path_range = [](const path_range_t& a, const path_range_t& b) -> bool { + if (a.begin.path != b.begin.path) return a.begin.path < b.begin.path; + if (a.begin.offset != b.begin.offset) return a.begin.offset < b.begin.offset; + if (a.end.path != b.end.path) return a.end.path < b.end.path; + return a.end.offset < b.end.offset; + }; // Lambda function to compare two path_range_t objects + + std::sort(copy_ranges.begin(), copy_ranges.end(), compare_path_range); // Sort the copied vector using the lambda function + + for (size_t i = 1; i < copy_ranges.size(); i++) { + if (!compare_path_range(copy_ranges[i-1], copy_ranges[i])) { + std::cerr << "[odgi::extract] error: " << graph.get_path_name(copy_ranges[i].begin.path) << ":" << copy_ranges[i].begin.offset << "-" << copy_ranges[i].end.offset << " is a duplicated path range" << std::endl; + return 1; + } } } @@ -384,6 +404,7 @@ namespace odgi { const bool show_progress = args::get(_show_progress); + const bool optimize = args::get(_optimize); const uint64_t context_steps = _context_steps ? args::get(_context_steps) : 0; const uint64_t context_bases = _context_bases ? args::get(_context_bases) : 0; @@ -395,7 +416,7 @@ namespace odgi { std::vector path_ranges, std::vector> pangenomic_ranges, const uint64_t context_steps, const uint64_t context_bases, const bool full_range, const bool inverse, const uint64_t max_dist_subpaths, const uint64_t num_iterations, - const uint64_t num_threads, const bool show_progress) { + const uint64_t num_threads, const bool show_progress, const bool optimize) { if (context_steps > 0 || context_bases > 0) { if (show_progress) { std::cerr << "[odgi::extract] expansion and adding connecting edges" << std::endl; @@ -456,11 +477,37 @@ namespace odgi { if (show_progress) { progress->increment(1); } - algorithms::for_handle_in_path_range( - source, path_handle, path_range.begin.offset, path_range.end.offset, - [&](const handle_t& handle) { - keep_bv.set(source.get_id(handle) - shift); - }); + + // The extraction does not cut nodes, so the input path ranges have to be + // extended if their ranges (start, end) fall in the middle of the nodes. + bool first = true; + uint64_t new_start = 0; + uint64_t new_end = 0; + + const uint64_t start = path_range.begin.offset; + const uint64_t end = path_range.end.offset; + + uint64_t walked = 0; + const auto path_end = source.path_end(path_handle); + for (step_handle_t cur_step = source.path_begin(path_handle); + cur_step != path_end && walked < end; cur_step = source.get_next_step(cur_step)) { + const handle_t cur_handle = source.get_handle_of_step(cur_step); + walked += source.get_length(cur_handle); + if (walked > start) { + keep_bv.set(source.get_id(cur_handle) - shift); + + if (first) { + first = false; + new_start = walked - source.get_length(cur_handle); + } + } + } + new_end = walked; + + // Extend path range to entirely include the first and the last node of the range. + // Thi is important to path names with the correct path ranges. + path_range.begin.offset = new_start; + path_range.end.offset = new_end; } if (!pangenomic_ranges.empty()) { uint64_t pos = 0; @@ -610,6 +657,7 @@ namespace odgi { const std::string path_name = source.get_path_name(path_range.begin.path); subpaths_from_path_ranges.push_back( + // The function assumes that every path is new and unique odgi::algorithms::create_subpath( subgraph, odgi::algorithms::make_path_name(path_name, path_range.begin.offset, path_range.end.offset), @@ -717,6 +765,10 @@ namespace odgi { // This should not be necessary, if the extraction works correctly // subgraph.remove_orphan_edges(); + + if (optimize) { + subgraph.optimize(); + } }; auto check_and_create_handle = [&](const graph_t &source, graph_t &subgraph, const nid_t node_id) { @@ -744,7 +796,13 @@ namespace odgi { << path_range.end.offset << std::endl; } - prep_graph(graph, &paths, lace_paths, subgraph, {path_range}, *pangenomic_ranges, context_steps, context_bases, _full_range, false, args::get(_max_dist_subpaths), num_iterations, num_threads, show_progress); + prep_graph( + graph, &paths, + lace_paths, subgraph, + {path_range}, *pangenomic_ranges, + context_steps, context_bases, _full_range, false, + max_dist_subpaths, num_iterations, + num_threads, show_progress, optimize); const string filename = graph.get_path_name(path_range.begin.path) + ":" + to_string(path_range.begin.offset) + "-" + to_string(path_range.end.offset) + ".og"; @@ -771,7 +829,13 @@ namespace odgi { } } - prep_graph(graph, &paths, lace_paths, subgraph, *path_ranges, *pangenomic_ranges, context_steps, context_bases, _full_range, _inverse, args::get(_max_dist_subpaths), num_iterations, num_threads, show_progress); + prep_graph( + graph, &paths, + lace_paths, subgraph, + *path_ranges, *pangenomic_ranges, + context_steps, context_bases, _full_range, _inverse, + max_dist_subpaths, num_iterations, + num_threads, show_progress, optimize); { const std::string outfile = args::get(og_out_file); diff --git a/src/subcommand/paths_main.cpp b/src/subcommand/paths_main.cpp index 63223b36b..f531ea4cb 100644 --- a/src/subcommand/paths_main.cpp +++ b/src/subcommand/paths_main.cpp @@ -45,7 +45,11 @@ int main_paths(int argc, char** argv) { {'D', "delim"}); args::ValueFlag path_delim_pos(path_investigation_opts, "N", "Consider the N-th occurrence of the delimiter specified with **-D, --delim**" " to obtain the group identifier. Specify 1 for the 1st occurrence (default).", - {'p', "delim-pos"}); + {'p', "delim-pos"}); + args::ValueFlag non_reference_nodes(path_investigation_opts, "FILE", "Print to stdout IDs of nodes that are not in the paths listed (by line) in *FILE*.", {"non-reference-nodes"}); + args::ValueFlag non_reference_ranges(path_investigation_opts, "FILE", "Print to stdout (in BED format) path ranges that are not in the paths listed (by line) in *FILE*.", {"non-reference-ranges"}); + args::ValueFlag min_size(path_investigation_opts, "N", "Minimum size (in bps) of nodes (with --non-reference-nodes) or ranges (with --non-reference-ranges).", {"min-size"}); + args::Flag show_step_ranges(path_investigation_opts, "N", "Show steps (that is, node IDs and strands) of --non-reference-ranges.", {"show-step-ranges"}); args::Group path_modification_opts(parser, "[ Path Modification Options ]"); args::ValueFlag keep_paths_file(path_modification_opts, "FILE", "Keep paths listed (by line) in *FILE*.", {'K', "keep-paths"}); args::ValueFlag drop_paths_file(path_modification_opts, "FILE", "Drop paths listed (by line) in *FILE*.", {'X', "drop-paths"}); @@ -92,6 +96,11 @@ int main_paths(int argc, char** argv) { return 1; } + if (non_reference_nodes && non_reference_ranges) { + std::cerr << "[odgi::paths] error: specify --non-reference-nodes or --non-reference-ranges, not both." << std::endl; + return 1; + } + const uint64_t num_threads = args::get(threads) ? args::get(threads) : 1; omp_set_num_threads(num_threads); @@ -106,6 +115,19 @@ int main_paths(int argc, char** argv) { } } + const uint64_t shift = graph.min_node_id(); + if ( + args::get(haplo_matrix) || + (non_reference_nodes && !args::get(non_reference_nodes).empty()) || + (non_reference_ranges && !args::get(non_reference_ranges).empty()) + ) { + // Check if the node IDs are compacted + if (graph.max_node_id() - shift >= graph.get_node_count()){ + std::cerr << "[odgi::paths] error: the node IDs are not compacted. Please run 'odgi sort' using -O, --optimize to optimize the graph." << std::endl; + exit(1); + } + } + if (list_path_start_end && list_names) { std::vector paths; graph.for_each_path_handle([&](const path_handle_t& p) { @@ -195,13 +217,18 @@ int main_paths(int argc, char** argv) { uint64_t path_length = 0; uint64_t path_step_count = 0; std::vector row(graph.get_node_count()); + // Initialize first to avoid possible bugs later + for (uint32_t i = 0; i < graph.get_node_count(); ++i) { + row[i] = 0; + } + graph.for_each_step_in_path( p, [&](const step_handle_t& s) { const handle_t& h = graph.get_handle_of_step(s); path_length += graph.get_length(h); ++path_step_count; - row[graph.get_id(h)-1]++; + row[graph.get_id(h)-shift]++; }); if (delim) { std::cout << group_name << "\t"; @@ -211,7 +238,7 @@ int main_paths(int argc, char** argv) { << path_step_count; if (node_length_scale) { for (uint64_t i = 0; i < row.size(); ++i) { - std::cout << "\t" << row[i] * graph.get_length(graph.get_handle(i+1)); + std::cout << "\t" << row[i] * graph.get_length(graph.get_handle(i+shift)); } } else { for (uint64_t i = 0; i < row.size(); ++i) { @@ -347,6 +374,140 @@ int main_paths(int argc, char** argv) { } } + if ( + (non_reference_nodes && !args::get(non_reference_nodes).empty()) || + (non_reference_ranges && !args::get(non_reference_ranges).empty()) + ) { + const uint64_t min_size_in_bp = min_size ? args::get(min_size) : 0; + + // Read paths to use as reference paths + std::vector reference_paths; + std::string line; + auto& x = non_reference_nodes && !args::get(non_reference_nodes).empty() ? args::get(non_reference_nodes) : args::get(non_reference_ranges); + std::ifstream infile(x); + while (std::getline(infile, line)) { + // This file should contain path names, one per line + auto& name = line; + if (graph.has_path(name)) { + reference_paths.push_back(graph.get_path_handle(name)); + } else { + std::cerr << "[odgi::paths] error: path'" << name + << "' does not exist in graph." << std::endl; + return 1; + } + } + + if (non_reference_nodes && !args::get(non_reference_nodes).empty()){ + // Emit non-reference nodes + + // Set non-reference nodes + atomicbitvector::atomic_bv_t non_reference_nodes(graph.get_node_count()); + for(uint64_t x = 0; x < non_reference_nodes.size(); ++x) { + if (min_size_in_bp == 0 || graph.get_length(graph.get_handle(x + shift)) >= min_size_in_bp) { + non_reference_nodes.set(x); + } + } +#pragma omp parallel for schedule(dynamic,1) + for (auto &path : reference_paths) { + graph.for_each_step_in_path(path, [&](const step_handle_t& step) { + const handle_t handle = graph.get_handle_of_step(step); + non_reference_nodes.reset(graph.get_id(handle) - shift); + }); + } + + // Emit non-reference nodes + std::cout << "#node.id\tnode.len\tpaths" << std::endl; + for (auto x : non_reference_nodes) { + const handle_t handle = graph.get_handle(x + shift); + + // Check paths that go through this node, if any + std::unordered_set unique_path_handles; + graph.for_each_step_on_handle(handle, [&](const step_handle_t& step) { + unique_path_handles.insert(graph.get_path_handle_of_step(step)); + }); + std::string result; + for (const auto& path : unique_path_handles) { + if (!result.empty()) { + result += ","; + } + result += graph.get_path_name(path); + } + + std::cout << graph.get_id(handle) << "\t" << graph.get_length(handle) << "\t" << result << std::endl; + } + } else { + // Emit non-reference ranges + + const bool _show_step_ranges = args::get(show_step_ranges); + + // Set the reference nodes + atomicbitvector::atomic_bv_t reference_nodes(graph.get_node_count()); + #pragma omp parallel for schedule(dynamic,1) + for (auto &path : reference_paths) { + graph.for_each_step_in_path(path, [&](const step_handle_t& step) { + const handle_t handle = graph.get_handle_of_step(step); + reference_nodes.set(graph.get_id(handle) - shift); + }); + } + + // Prepare non reference path handles for parallel processing + std::vector non_reference_paths; + graph.for_each_path_handle([&non_reference_paths](const path_handle_t& path) { + non_reference_paths.push_back(path); + }); + std::sort(non_reference_paths.begin(), non_reference_paths.end()); + std::sort(reference_paths.begin(), reference_paths.end()); + non_reference_paths.erase( + std::remove_if(non_reference_paths.begin(), non_reference_paths.end(), + [&reference_paths](const auto &x) { + return std::binary_search(reference_paths.begin(), reference_paths.end(), x); + }), non_reference_paths.end()); + + // Traverse non reference paths to emit non-reference ranges + if (_show_step_ranges) { + std::cout << "#path.name\tstart\tend\tsteps" << std::endl; + } else { + std::cout << "#path.name\tstart\tend" << std::endl; + } + #pragma omp parallel for schedule(dynamic, 1) + for (auto& path : non_reference_paths) { + uint64_t start = 0, end = 0; + std::vector step_range; + graph.for_each_step_in_path(path, [&](const step_handle_t& step) { + const handle_t handle = graph.get_handle_of_step(step); + const uint64_t index = graph.get_id(handle) - shift; + if (reference_nodes.test(index)) { + // Emit the previous non reference range, if any + if (end > start && (end - start) >= min_size_in_bp) { + if (_show_step_ranges) { + std::string step_range_str = ""; + for (auto& step : step_range) { + const handle_t handle = graph.get_handle_of_step(step); + step_range_str += std::to_string(graph.get_id(handle)) + (graph.get_is_reverse(handle) ? "-" : "+") + ","; + } + #pragma omp critical (cout) + std::cout << graph.get_path_name(path) << "\t" << start << "\t" << end << "\t" << step_range_str.substr(0, step_range_str.size() - 1) << std::endl; // trim the trailing comma from step_range + } else { + #pragma omp critical (cout) + std::cout << graph.get_path_name(path) << "\t" << start << "\t" << end << std::endl; + } + } + end += graph.get_length(handle); + start = end; + if (_show_step_ranges) { + step_range.clear(); + } + } else { + end += graph.get_length(handle); + } + if (_show_step_ranges) { + step_range.push_back(step); + } + }); + } + } + } + return 0; } diff --git a/src/subcommand/pav_main.cpp b/src/subcommand/pav_main.cpp index 02b2b81db..b8eeffbe2 100644 --- a/src/subcommand/pav_main.cpp +++ b/src/subcommand/pav_main.cpp @@ -119,7 +119,7 @@ int main_pav(int argc, char **argv) { auto vals = split(line, '\t'); if (vals.size() != 2) { std::cerr << "[odgi::pav] line does not have a path.name and path.group value:" - << std::endl << line << std::endl; + << std::endl << line << std::endl; return 1; } auto& path_name = vals.front(); @@ -136,8 +136,8 @@ int main_pav(int argc, char **argv) { if (group_2_index.empty()) { std::cerr - << "[odgi::pav] error: 0 path groups were read. Please specify at least one path group via -p/--path-groups." - << std::endl; + << "[odgi::pav] error: 0 path groups were read. Please specify at least one path group via -p/--path-groups." + << std::endl; return 1; } } else if (_group_by_sample) { @@ -185,7 +185,6 @@ int main_pav(int argc, char **argv) { add_bed_range(path_ranges, graph, line); } } - if (path_ranges.empty()) { std::cerr << "[odgi::pav] error: please specify path ranges via -b/--bed-file." @@ -225,7 +224,12 @@ int main_pav(int argc, char **argv) { path_handle_2_index[p2mm.first] = i++; } - // Prepare the interval trees to query target path ranges + // Prepare the interval trees for querying target path ranges + std::unique_ptr operation_progress; + if (show_progress) { + std::string banner = "[odgi::pav] preparing the interval trees for querying target path ranges:"; + operation_progress = std::make_unique(path_handles.size(), banner); + } std::vector> trees; trees.resize(path_handles.size()); #pragma omp parallel for schedule(dynamic, 1) num_threads(num_threads) @@ -250,8 +254,14 @@ int main_pav(int argc, char **argv) { } tree.index(); // index - } + if (show_progress) { + operation_progress->increment(1); + } + } + if (show_progress) { + operation_progress->finish(); + } const bool emit_matrix_else_table = args::get(_matrix_output); // Emit the PAV matrix @@ -296,6 +306,11 @@ int main_pav(int argc, char **argv) { << (emit_binary_values ? pav_ratio >= binary_threshold : pav_ratio) << "\n"; }; + if (show_progress) { + std::string banner = "[odgi::pav] emitting PAV results:"; + operation_progress = std::make_unique(path_ranges.size(), banner); + } + #pragma omp parallel for schedule(dynamic, 1) num_threads(num_threads) for (uint64_t i = 0; i < path_ranges.size(); ++i) { auto &path_range = path_ranges[i]; @@ -386,8 +401,15 @@ int main_pav(int argc, char **argv) { } } } - } + + if (show_progress) { + operation_progress->increment(1); + } + } + if (show_progress) { + operation_progress->finish(); + } return 0; } diff --git a/src/subcommand/similarity_main.cpp b/src/subcommand/similarity_main.cpp index 0d241d2dd..da913e3ce 100644 --- a/src/subcommand/similarity_main.cpp +++ b/src/subcommand/similarity_main.cpp @@ -162,7 +162,7 @@ args::Group threading_opts(parser, "[ Threading ]"); if (using_delim) { bp_count.resize(path_groups.size()); } else { - bp_count.resize(graph.get_path_count()); + bp_count.resize(graph.get_path_count() + 1); } uint32_t path_max = 0; diff --git a/src/subcommand/sort_main.cpp b/src/subcommand/sort_main.cpp index 3347e59d1..109b9e123 100644 --- a/src/subcommand/sort_main.cpp +++ b/src/subcommand/sort_main.cpp @@ -271,6 +271,7 @@ int main_sort(int argc, char** argv) { std::string banner = "[odgi::sort] preparing target path vectors:"; target_paths_progress = std::make_unique(target_paths.size(), banner); } + uint64_t ref_nodes = 0; for (handlegraph::path_handle_t target_path: target_paths) { graph.for_each_step_in_path( target_path, @@ -280,6 +281,7 @@ int main_sort(int argc, char** argv) { if (!is_ref[i]) { is_ref[i] = true; target_order.push_back(handle); + ref_nodes++; } }); if (args::get(progress)) { @@ -289,12 +291,10 @@ int main_sort(int argc, char** argv) { if (args::get(progress)) { target_paths_progress->finish(); } - uint64_t ref_nodes = 0; for (uint64_t i = 0; i < is_ref.size(); i++) { bool ref = is_ref[i]; if (!ref) { target_order.push_back(graph.get_handle(i + 1)); - ref_nodes++; } } graph.apply_ordering(target_order, true); diff --git a/src/subcommand/stats_main.cpp b/src/subcommand/stats_main.cpp index 6ee48b7cb..15d06110d 100644 --- a/src/subcommand/stats_main.cpp +++ b/src/subcommand/stats_main.cpp @@ -167,7 +167,7 @@ int main_stats(int argc, char** argv) { const uint64_t shift = number_bool_packing::unpack_number(graph.get_handle(graph.min_node_id())); - if (args::get(mean_links_length) || args::get(sum_of_path_node_distances) || _multiqc) { + if (args::get(mean_links_length) || args::get(sum_of_path_node_distances)) { if (number_bool_packing::unpack_number(graph.get_handle(graph.max_node_id())) - shift >= graph.get_node_count()){ std::cerr << "[odgi::stats] error: the node IDs are not compacted. Please run 'odgi sort' using -O, --optimize to optimize the graph." << std::endl; exit(1); @@ -371,7 +371,7 @@ int main_stats(int argc, char** argv) { // TODO clear all sets? } - if (args::get(mean_links_length) || args::get(sum_of_path_node_distances) || _multiqc) { + if (args::get(mean_links_length) || args::get(sum_of_path_node_distances)) { // This vector is needed for computing the metrics in 1D and for detecting gap-links std::vector position_map(graph.get_node_count() + 1); @@ -409,7 +409,7 @@ int main_stats(int argc, char** argv) { position_map[position_map.size() - 1] = len; } - if (args::get(mean_links_length) || _multiqc){ + if (args::get(mean_links_length)){ bool _dont_penalize_gap_links = args::get(dont_penalize_gap_links); uint64_t sum_all_node_space = 0; @@ -582,7 +582,7 @@ int main_stats(int argc, char** argv) { } } - if (args::get(sum_of_path_node_distances) || _multiqc){ + if (args::get(sum_of_path_node_distances)){ bool _penalize_diff_orientation = args::get(penalize_diff_orientation); uint64_t sum_all_path_node_dist_node_space = 0;