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.

4 files changed, 80 insertions, 15 deletions

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 @@
         <term><option>--path-range=<replaceable>PATH_RANGE</replaceable></option></term>
         <listitem><para>Extract segments in <replaceable>PATH_RANGE</replaceable>, specified in the <replaceable>path[:pos1[-pos2]]</replaceable> format. <replaceable>pos1</replaceable> and <replaceable>pos2</replaceable> are 0-based coordinates. The extracted segments include <replaceable>pos1</replaceable> (inclusive) but not <replaceable>pos2</replaceable> (exclusive).</para></listitem>
       </varlistentry>
+      <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>
+      </varlistentry>
     </variablelist>
   </refsection>
 </refentry>
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")),