Bridge extracted subpaths. HEAD main

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",