From 1783136fa37cb502d07ce2e4bb0474cd3a0ff1f2 Mon Sep 17 00:00:00 2001
From: Arun Isaac
Date: Fri, 28 Aug 2026 00:46:37 +0100
Subject: Split subpaths selected by extract path range into discrete
fragments.
When `domagi extract` is given a path range, it may select
discontinuous fragments of other paths. These fragments must be output
as separate paths, not erroneously merged together.
---
doc/domagi-extract.dbk | 5 ++
domagi/domagi.py | 60 ++++++++++++++++------
.../test1-extract-path-no-bridge-10-20.gfa | 9 ++++
tests/test_domagi.py | 21 ++++++++
4 files changed, 80 insertions(+), 15 deletions(-)
create mode 100644 test-data/expected-output/test1-extract-path-no-bridge-10-20.gfa
diff --git a/doc/domagi-extract.dbk b/doc/domagi-extract.dbk
index d20afd9..9671d33 100644
--- a/doc/domagi-extract.dbk
+++ b/doc/domagi-extract.dbk
@@ -32,6 +32,11 @@
Extract segments in PATH_RANGE, specified in the path[:pos1[-pos2]] format. pos1 and pos2 are 0-based coordinates. The extracted segments include pos1 (inclusive) but not pos2 (exclusive).
+
+
+
+ Bridge subpaths that are separated by less than DISTANCE. Default DISTANCE is 300000. This reduces the fragmentation of paths that are unspecified in the input path ranges. Set DISTANCE to 0 to disable bridging.
+
diff --git a/domagi/domagi.py b/domagi/domagi.py
index c369b38..f398926 100644
--- a/domagi/domagi.py
+++ b/domagi/domagi.py
@@ -233,8 +233,12 @@ def depth(con, graph_depth_table, paths, bed_input, threads, progress):
type=click.STRING,
metavar="PATH[:POS1[-POS2]]",
help="extract segments in path range")
+@click.option("-d", "--max-distance-subpaths",
+ type=click.INT,
+ default=300*10**3,
+ help="bridge subpaths that are separated by less than this distance (default: 300000)")
@common_options
-def extract(con, outfile, segment_name, path_range, steps, threads, progress):
+def extract(con, outfile, segment_name, path_range, steps, max_distance_subpaths, threads, progress):
set_duckdb_threads(con, threads)
with connect_duckdb(outfile, threads) as out_con:
out_con.execute(read_sql("schema.sql"))
@@ -259,34 +263,60 @@ def extract(con, outfile, segment_name, path_range, steps, threads, progress):
FROM selected_segment
INNER JOIN link ON from_segment=selected_segment.id;
+ CREATE TEMPORARY TABLE selected_path_segment AS
+ WITH selected_path_segment_with_gap_flag AS (
+ -- Select path segments marking gaps within paths using a gap flag.
+ SELECT path_segment.path_id AS parent_path_id,
+ segment_id,
+ segment_orientation,
+ start,
+ "end",
+ CASE WHEN LAG("end") OVER (PARTITION BY path_id ORDER BY start)=start
+ THEN 0
+ ELSE 1
+ END AS gap
+ FROM selected_segment
+ INNER JOIN path_segment ON path_segment.segment_id=selected_segment.id)
+ -- Sum over gap flags to get the new path ID.
+ SELECT SUM(gap) OVER (ORDER BY parent_path_id, start)-1 AS path_id,
+ parent_path_id,
+ segment_id,
+ segment_orientation,
+ start,
+ "end"
+ FROM selected_path_segment_with_gap_flag;
+
CREATE TEMPORARY TABLE path_first_start AS
+ -- Find the first start coordinate of each new path.
SELECT path_id, MIN(start) AS first_start
- FROM selected_segment
- INNER JOIN path_segment ON path_segment.segment_id=selected_segment.id
+ FROM selected_path_segment
GROUP BY path_id;
INSERT INTO subset_db.path_segment
- -- There are no gaps in the path segments. So, it's enough to subtract
- -- first_start from the path segment coordinates.
- SELECT path_segment.path_id,
+ -- The gaps in the path segments have been resolved by splitting paths into
+ -- smaller continuous paths. So, it's enough to subtract first_start from
+ -- the path segment coordinates.
+ SELECT selected_path_segment.path_id,
segment_id,
segment_orientation,
start-first_start AS start,
- "end"-first_start AS end
- FROM selected_segment
- INNER JOIN path_segment ON path_segment.segment_id=selected_segment.id
- INNER JOIN path_first_start ON path_first_start.path_id=path_segment.path_id
+ "end"-first_start AS "end"
+ FROM selected_path_segment
+ INNER JOIN path_first_start
+ ON path_first_start.path_id=selected_path_segment.path_id
-- Re-order similar to post-import.sql for optimal access.
ORDER BY path_id, start, "end";
INSERT INTO subset_db.path
- SELECT id, ANY_VALUE(name) || ':' || MIN(start)+ANY_VALUE(first_start) || '-' || MAX("end")+ANY_VALUE(first_start)
- FROM subset_db.path_segment
- INNER JOIN path ON subset_db.path_segment.path_id=path.id
- INNER JOIN path_first_start ON path_first_start.path_id=path.id
- GROUP BY id;
+ -- Write the new paths labelling them based on their parent paths.
+ SELECT selected_path_segment.path_id,
+ ANY_VALUE(name) || ':' || MIN(start) || '-' || MAX("end") AS path_name
+ FROM selected_path_segment
+ INNER JOIN path ON path.id=selected_path_segment.parent_path_id
+ GROUP BY selected_path_segment.path_id;
DROP TABLE path_first_start;
+ DROP TABLE selected_path_segment;
DROP TABLE selected_segment;
DETACH subset_db;
""")
diff --git a/test-data/expected-output/test1-extract-path-no-bridge-10-20.gfa b/test-data/expected-output/test1-extract-path-no-bridge-10-20.gfa
new file mode 100644
index 0000000..4846f99
--- /dev/null
+++ b/test-data/expected-output/test1-extract-path-no-bridge-10-20.gfa
@@ -0,0 +1,9 @@
+H VN:Z:1.0
+S 6 TTG
+L 6 + 8 + 0M
+S 8 G
+L 8 + 9 + 0M
+S 9 AAATTTTCTGGAGTTCTAT
+P x:10-33 6+,8+,9+ *
+P y:10-13 6+ *
+P y:14-33 9+ *
diff --git a/tests/test_domagi.py b/tests/test_domagi.py
index 8f92de6..ade24cd 100644
--- a/tests/test_domagi.py
+++ b/tests/test_domagi.py
@@ -265,6 +265,27 @@ def test_domagi_extract_path(tmp_path, request, domagi_db_name, path_range, expe
with open(expected_output) as file:
assert_gfa_equal(file, io.StringIO(result.stdout))
+@pytest.mark.parametrize("domagi_db_name, path_range, expected_output",
+ [("domagi_db_test1",
+ "x:10-20",
+ Path("test-data/expected-output/test1-extract-path-no-bridge-10-20.gfa"))])
+def test_domagi_extract_path_no_bridge(tmp_path, request, domagi_db_name, path_range, expected_output):
+ domagi_db = request.getfixturevalue(domagi_db_name)
+ output_duckdb_path = tmp_path / f"{domagi_db.stem}-output.db"
+ runner = CliRunner()
+ result = runner.invoke(main, ["extract",
+ "--db", domagi_db,
+ "--path-range", path_range,
+ "--max-distance-subpaths", "0",
+ "--out", output_duckdb_path])
+ assert result.exit_code == 0
+ result = runner.invoke(main, ["view",
+ "--to-gfa",
+ "--db", output_duckdb_path])
+ assert result.exit_code == 0
+ with open(expected_output) as file:
+ assert_gfa_equal(file, io.StringIO(result.stdout))
+
@pytest.mark.parametrize("domagi_db_name, expected_output",
[("domagi_db_test1",
Path("test-data/expected-output/test1-matrix")),
--
cgit 1.4.1