diff options
| -rw-r--r-- | doc/domagi-extract.dbk | 5 | ||||
| -rw-r--r-- | domagi/domagi.py | 76 | ||||
| -rw-r--r-- | test-data/expected-output/test-extract-path-bridge-0.gfa (renamed from test-data/expected-output/test1-extract-path-no-bridge-10-20.gfa) | 2 | ||||
| -rw-r--r-- | test-data/expected-output/test-extract-path-bridge-1.gfa | 9 | ||||
| -rw-r--r-- | test-data/expected-output/test-extract-path-bridge-2.gfa | 9 | ||||
| -rw-r--r-- | test-data/expected-output/test-extract-path-bridge-3.gfa | 13 | ||||
| -rw-r--r-- | test-data/expected-output/test1-extract-bed-windows.gfa | 46 | ||||
| -rw-r--r-- | test-data/expected-output/test2-extract-bed-windows.gfa | 12 | ||||
| -rw-r--r-- | test-data/expected-output/test3-extract-bed-windows.gfa | 13 | ||||
| -rw-r--r-- | test-data/test-extract-path-bridge.gfa | 18 | ||||
| -rw-r--r-- | tests/test_domagi.py | 28 |
11 files changed, 203 insertions, 28 deletions
diff --git a/doc/domagi-extract.dbk b/doc/domagi-extract.dbk index 9671d33..736703d 100644 --- a/doc/domagi-extract.dbk +++ b/doc/domagi-extract.dbk @@ -35,7 +35,10 @@ <varlistentry> <term><option>-d <replaceable>DISTANCE</replaceable></option></term> <term><option>--max-distance-subpaths=<replaceable>DISTANCE</replaceable></option></term> - <listitem><para>Bridge subpaths that are separated by less than <replaceable>DISTANCE</replaceable>. Default <replaceable>DISTANCE</replaceable> is <literal>300000</literal>. This reduces the fragmentation of paths that are unspecified in the input path ranges. Set <replaceable>DISTANCE</replaceable> to <literal>0</literal> to disable bridging.</para></listitem> + <listitem> + <para>Bridge subpaths that are separated by less than <replaceable>DISTANCE</replaceable>. Default <replaceable>DISTANCE</replaceable> is <literal>300000</literal>. This reduces the fragmentation of paths that are unspecified in the input path ranges. Set <replaceable>DISTANCE</replaceable> to <literal>0</literal> to disable bridging.</para> + <para>In contrast to odgi, when bridging subpaths, domagi does not use the newly extracted segments to extend the subgraph further and recursively bridge more subpaths.</para> + </listitem> </varlistentry> </variablelist> </refsection> diff --git a/domagi/domagi.py b/domagi/domagi.py index f398926..a221a6b 100644 --- a/domagi/domagi.py +++ b/domagi/domagi.py @@ -254,44 +254,82 @@ def extract(con, outfile, segment_name, path_range, steps, max_distance_subpaths con.execute(f""" ATTACH '{outfile}' AS subset_db (READ_WRITE); - INSERT INTO subset_db.segment - SELECT segment.id, name, sequence FROM selected_segment - INNER JOIN segment ON segment.id=selected_segment.id; - - INSERT INTO subset_db.link - SELECT from_segment, from_orientation, to_segment, to_orientation - 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 ( + WITH selected_path_segment_with_gap AS ( + -- Select path segments computing gaps within paths. + SELECT path_segment.path_id, + segment_id, + segment_orientation, + start, + "end", + start-LAG("end") OVER (PARTITION BY path_id ORDER BY start) AS gap + FROM path_segment + INNER JOIN selected_segment ON selected_segment.id=path_segment.segment_id), + gap AS ( + -- Find gaps within paths. + SELECT path_id, + LAG("end") OVER (PARTITION BY path_id ORDER BY start) AS start, + start AS "end" + FROM selected_path_segment_with_gap + QUALIFY gap<=?), + selected_path_segment_with_bridged_gaps AS ( + -- Find segments in gaps, and union them with already selected + -- segments. + SELECT path_segment.path_id, + segment_id, + segment_orientation, + path_segment.start, + path_segment.end + FROM path_segment + INNER JOIN gap + ON gap.path_id=path_segment.path_id + AND gap.start<path_segment.end AND path_segment.start<gap.end + UNION ALL + SELECT path_id, segment_id, segment_orientation, start, "end" + FROM selected_path_segment_with_gap), + 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, + SELECT 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) + THEN 0 ELSE 1 END AS gap + FROM selected_path_segment_with_bridged_gaps) -- 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, + SELECT SUM(gap) OVER (ORDER BY path_id, start)-1 AS path_id, + path_id AS parent_path_id, segment_id, segment_orientation, start, "end" FROM selected_path_segment_with_gap_flag; - + """, + [max_distance_subpaths]) + con.execute(""" 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_path_segment GROUP BY path_id; + INSERT INTO subset_db.segment + SELECT segment.id, name, sequence + FROM segment + INNER JOIN (SELECT DISTINCT segment_id FROM selected_path_segment) AS path_segment + ON path_segment.segment_id=segment.id; + + INSERT INTO subset_db.link + -- Only select links where both the from and to segments are part of the + -- extracted subgraph. + SELECT from_segment, from_orientation, to_segment, to_orientation + FROM link + INNER JOIN subset_db.segment segment1 + ON segment1.id=link.from_segment + INNER JOIN subset_db.segment segment2 + ON segment2.id=link.to_segment; + INSERT INTO subset_db.path_segment -- 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 diff --git a/test-data/expected-output/test1-extract-path-no-bridge-10-20.gfa b/test-data/expected-output/test-extract-path-bridge-0.gfa index 4846f99..96dae01 100644 --- a/test-data/expected-output/test1-extract-path-no-bridge-10-20.gfa +++ b/test-data/expected-output/test-extract-path-bridge-0.gfa @@ -6,4 +6,4 @@ L 8 + 9 + 0M S 9 AAATTTTCTGGAGTTCTAT P x:10-33 6+,8+,9+ * P y:10-13 6+ * -P y:14-33 9+ * +P y:16-35 9+ * diff --git a/test-data/expected-output/test-extract-path-bridge-1.gfa b/test-data/expected-output/test-extract-path-bridge-1.gfa new file mode 100644 index 0000000..96dae01 --- /dev/null +++ b/test-data/expected-output/test-extract-path-bridge-1.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:16-35 9+ * diff --git a/test-data/expected-output/test-extract-path-bridge-2.gfa b/test-data/expected-output/test-extract-path-bridge-2.gfa new file mode 100644 index 0000000..96dae01 --- /dev/null +++ b/test-data/expected-output/test-extract-path-bridge-2.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:16-35 9+ * diff --git a/test-data/expected-output/test-extract-path-bridge-3.gfa b/test-data/expected-output/test-extract-path-bridge-3.gfa new file mode 100644 index 0000000..7ccf7e8 --- /dev/null +++ b/test-data/expected-output/test-extract-path-bridge-3.gfa @@ -0,0 +1,13 @@ +H VN:Z:1.0 +S 6 TTG +L 6 + 7 + 0M +L 6 + 8 + 0M +S 7 AG +L 7 + 10 + 0M +S 8 G +L 8 + 9 + 0M +S 9 AAATTTTCTGGAGTTCTAT +S 10 A +L 10 + 9 + 0M +P x:10-33 6+,8+,9+ * +P y:10-35 6+,7+,10+,9+ * diff --git a/test-data/expected-output/test1-extract-bed-windows.gfa b/test-data/expected-output/test1-extract-bed-windows.gfa new file mode 100644 index 0000000..ba2bd79 --- /dev/null +++ b/test-data/expected-output/test1-extract-bed-windows.gfa @@ -0,0 +1,46 @@ +H VN:Z:1.0 +S 1 CAAATAAG +L 1 + 2 + 0M +L 1 + 3 + 0M +S 2 A +L 2 + 4 + 0M +L 2 + 5 + 0M +S 3 G +L 3 + 4 + 0M +L 3 + 5 + 0M +S 4 T +L 4 + 6 + 0M +S 5 C +L 5 + 6 + 0M +S 6 TTG +L 6 + 7 + 0M +L 6 + 8 + 0M +S 7 A +L 7 + 9 + 0M +S 8 G +L 8 + 9 + 0M +S 9 AAATTTTCTGGAGTTCTAT +L 9 + 10 + 0M +L 9 + 11 + 0M +S 10 A +L 10 + 12 + 0M +S 11 T +L 11 + 12 + 0M +S 12 ATAT +L 12 + 13 + 0M +L 12 + 14 + 0M +S 13 A +L 13 + 15 + 0M +S 14 T +L 14 + 15 + 0M +S 15 CCAACTCTCTG +P x:0-10 1+,3+,5+ * +P x:10-33 6+,8+,9+ * +P x:14-33 9+ * +P x:14-50 9+,11+,12+,14+,15+ * +P x:39-50 15+ * +P y:0-10 1+,2+,4+ * +P y:10-33 6+,7+,9+ * +P y:14-33 9+ * +P y:14-50 9+,10+,12+,13+,15+ * +P y:39-50 15+ * diff --git a/test-data/expected-output/test2-extract-bed-windows.gfa b/test-data/expected-output/test2-extract-bed-windows.gfa new file mode 100644 index 0000000..bd0b466 --- /dev/null +++ b/test-data/expected-output/test2-extract-bed-windows.gfa @@ -0,0 +1,12 @@ +H VN:Z:1.0 +S 1 AGGA +L 1 + 3 + 0M +L 1 + 3 - 0M +S 3 TC +L 3 + 4 + 0M +S 4 TCTCAGG +L 4 - 3 + 0M +P 5+:0-13 1+,3+,4+ * +P 5+:6-13 4+ * +P 5-:0-13 1+,3-,4+ * +P 5-:6-13 4+ * diff --git a/test-data/expected-output/test3-extract-bed-windows.gfa b/test-data/expected-output/test3-extract-bed-windows.gfa new file mode 100644 index 0000000..8fbb0bc --- /dev/null +++ b/test-data/expected-output/test3-extract-bed-windows.gfa @@ -0,0 +1,13 @@ +H VN:Z:1.0 +S 1 AGGA +L 1 + 3 + 0M +L 1 + 3 - 0M +S 3 TC +L 3 + 4 + 0M +S 4 TCTCAGG +L 4 - 3 + 0M +L 4 + 4 + 0M +P 5+:0-13 1+,3+,4+ * +P 5+:6-20 4+,4+ * +P 5-:0-13 1+,3-,4+ * +P 5-:6-13 4+ * diff --git a/test-data/test-extract-path-bridge.gfa b/test-data/test-extract-path-bridge.gfa new file mode 100644 index 0000000..1f2846f --- /dev/null +++ b/test-data/test-extract-path-bridge.gfa @@ -0,0 +1,18 @@ +H VN:Z:1.0 +S 1 CAAATAAG +S 2 A +S 3 G +S 4 T +S 5 C +S 6 TTG +S 7 AG +S 8 G +S 9 AAATTTTCTGGAGTTCTAT +S 10 A +L 6 + 7 + +L 6 + 8 + +L 7 + 10 + +L 8 + 9 + +L 10 + 9 + +P x 1+,3+,5+,6+,8+,9+ +P y 1+,2+,4+,6+,7+,10+,9+ diff --git a/tests/test_domagi.py b/tests/test_domagi.py index ade24cd..d12ce42 100644 --- a/tests/test_domagi.py +++ b/tests/test_domagi.py @@ -80,6 +80,11 @@ def domagi_db_testcrush(tmp_path_factory): return build_db(Path("test-data/test-crush.gfa"), tmp_path_factory.mktemp("db") / "test-crush.db") +@pytest.fixture(scope="session") +def domagi_db_testextractpathbridge(tmp_path_factory): + return build_db(Path("test-data/test-extract-path-bridge.gfa"), + tmp_path_factory.mktemp("db") / "test-extract-path-bridge.db") + @pytest.mark.parametrize("domagi_db_name, chop_to, expected_output", [("domagi_db_test3", 1, @@ -265,18 +270,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): +@pytest.mark.parametrize("domagi_db_name, max_distance_subpaths, expected_output", + [("domagi_db_testextractpathbridge", + 0, + Path("test-data/expected-output/test-extract-path-bridge-0.gfa")), + ("domagi_db_testextractpathbridge", + 1, + Path("test-data/expected-output/test-extract-path-bridge-1.gfa")), + ("domagi_db_testextractpathbridge", + 2, + Path("test-data/expected-output/test-extract-path-bridge-2.gfa")), + ("domagi_db_testextractpathbridge", + 3, + Path("test-data/expected-output/test-extract-path-bridge-3.gfa"))]) +def test_domagi_extract_path_bridge(tmp_path, request, domagi_db_name, max_distance_subpaths, 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", + "--path-range", "x:10-20", + "--max-distance-subpaths", str(max_distance_subpaths), "--out", output_duckdb_path]) assert result.exit_code == 0 result = runner.invoke(main, ["view", |