From d2282f6b7453a023af92424309e8713f5c8f7ca7 Mon Sep 17 00:00:00 2001 From: Arun Isaac Date: Tue, 18 Aug 2026 14:44:46 +0100 Subject: Add chop subcommand. --- doc/domagi-chop.dbk | 22 +++++ doc/domagi.dbk | 1 + domagi/chop-1.sql | 21 +++++ domagi/chop-2.sql | 124 +++++++++++++++++++++++++++++ domagi/chop.sql | 12 --- domagi/domagi.py | 30 +++++++ meson.build | 2 + test-data/expected-output/test3-chop-1.gfa | 33 ++++++++ test-data/expected-output/test3-chop-2.gfa | 21 +++++ test-data/expected-output/test3-chop-3.gfa | 19 +++++ test-data/expected-output/test3-chop-4.gfa | 15 ++++ tests/test_domagi.py | 33 ++++++++ 12 files changed, 321 insertions(+), 12 deletions(-) create mode 100644 doc/domagi-chop.dbk create mode 100644 domagi/chop-1.sql create mode 100644 domagi/chop-2.sql delete mode 100644 domagi/chop.sql create mode 100644 test-data/expected-output/test3-chop-1.gfa create mode 100644 test-data/expected-output/test3-chop-2.gfa create mode 100644 test-data/expected-output/test3-chop-3.gfa create mode 100644 test-data/expected-output/test3-chop-4.gfa diff --git a/doc/domagi-chop.dbk b/doc/domagi-chop.dbk new file mode 100644 index 0000000..dc308e1 --- /dev/null +++ b/doc/domagi-chop.dbk @@ -0,0 +1,22 @@ + + + + + Description + Divide segments into smaller pieces while preserving the graph topology. + + Options + + + + + + Divide nodes that are longer than N base pairs into nodes no longer than N while preserving the graph topology. + + + + + + diff --git a/doc/domagi.dbk b/doc/domagi.dbk index 53f25e2..0888ac1 100644 --- a/doc/domagi.dbk +++ b/doc/domagi.dbk @@ -8,6 +8,7 @@ Reference + diff --git a/domagi/chop-1.sql b/domagi/chop-1.sql new file mode 100644 index 0000000..51e7203 --- /dev/null +++ b/domagi/chop-1.sql @@ -0,0 +1,21 @@ +-- Chop up segments longer than $1 while preserving the graph +-- topology. All segments are renamed. This operation is split across +-- the files chop-*.sql. +CREATE TEMPORARY TABLE segment_chop AS + WITH segment_chop_start AS ( + -- Chop up segments. + SELECT id AS parent_id, + generate_series(1, len(sequence), $1) AS starts, + unnest(starts) AS start, + sequence[start:start+$1-1] AS sequence, + generate_subscripts(starts, 1) AS chop_index + FROM segment) + -- Re-number all segment IDs and names with sequential integers. + -- We start id from 0 and name from 1. + SELECT row_number() over ()::INTEGER-1 AS id, + row_number() over ()::VARCHAR AS name, + sequence, + parent_id, + chop_index, + length(starts) as last_chop + FROM segment_chop_start; diff --git a/domagi/chop-2.sql b/domagi/chop-2.sql new file mode 100644 index 0000000..2be1e8a --- /dev/null +++ b/domagi/chop-2.sql @@ -0,0 +1,124 @@ +-- Chop up segments longer than $1 while preserving the graph +-- topology. All segments are renamed. This operation is split across +-- the files chop-*.sql. + +-- Write output segment table. +INSERT INTO output_db.segment +SELECT id, name, sequence FROM segment_chop; + +CREATE TEMPORARY MACRO invert_orientation(orientation) AS + CASE WHEN orientation='+' THEN '-' + ELSE '+' + END; + +WITH link_chop_first AS ( + -- Rewrite existing links to go from the first/last chop of the + -- "from segment" to the first/last chop of the "to segment" + -- based on the orientation. + SELECT from_chop.id AS from_segment, + from_orientation, + to_chop.id AS to_segment, + to_orientation + FROM link + INNER JOIN segment_chop AS from_chop ON from_chop.parent_id=link.from_segment + INNER JOIN segment_chop AS to_chop ON to_chop.parent_id=link.to_segment + WHERE (from_orientation='+' AND to_orientation='+' + AND from_chop.chop_index=from_chop.last_chop + AND to_chop.chop_index=1) + OR (from_orientation='+' AND to_orientation='-' + AND from_chop.chop_index=from_chop.last_chop + AND to_chop.chop_index=to_chop.last_chop) + OR (from_orientation='-' AND to_orientation='+' + AND from_chop.chop_index=1 + AND to_chop.chop_index=1) + OR (from_orientation='-' AND to_orientation='-' + AND from_chop.chop_index=1 + AND to_chop.chop_index=to_chop.last_chop)), + link_segment_with_orientation AS ( + -- Derive a list of (segment, orientation) tuples used in the + -- link table. + SELECT from_segment AS segment, + from_orientation AS orientation + FROM link + UNION + SELECT to_segment AS segment, + to_orientation AS orientation + FROM link), + link_chop_internal AS ( + -- Create new links for the other chops—links that were + -- internal to the original unchopped segment. + SELECT lag(id) OVER (PARTITION BY parent_id ORDER by chop_index) AS from_segment, + id AS to_segment, + parent_id + FROM segment_chop + QUALIFY from_segment IS NOT NULL), + link_chop_internal_with_orientation AS ( + -- Add orientation to the new links. + SELECT from_segment, + orientation AS from_orientation, + to_segment, + orientation AS to_orientation + FROM link_chop_internal + INNER JOIN link_segment_with_orientation + ON link_segment_with_orientation.segment=link_chop_internal.parent_id + WHERE orientation='+' + UNION ALL + -- Reverse the from/to for negatively oriented segments. + SELECT to_segment, + orientation AS from_orientation, + from_segment, + orientation AS to_orientation + FROM link_chop_internal + INNER JOIN link_segment_with_orientation + ON link_segment_with_orientation.segment=link_chop_internal.parent_id + WHERE orientation='-'), + link_uncanonical AS ( + -- Union all links—new and old—to new table. + SELECT * FROM link_chop_first + UNION ALL + SELECT * FROM link_chop_internal_with_orientation), + link_canonical AS ( + -- Canonicalize links ensuring that from_segment <= to_segment. + SELECT * + FROM link_uncanonical + WHERE from_segment <= to_segment + UNION ALL + SELECT to_segment AS from_segment, + invert_orientation(to_orientation), + from_segment AS to_segment, + invert_orientation(from_orientation) + FROM link_uncanonical + WHERE from_segment > to_segment) + INSERT INTO output_db.link + -- Write only the distinct links to the output link table. + SELECT DISTINCT * + FROM link_canonical; + +-- Copy path table. +INSERT INTO output_db.path +SELECT * from path; + +-- Split up path segments to reflect the chops. +WITH path_segment_chop AS ( + SELECT path_id, + id AS segment_id, + segment_orientation, + start+$1*(chop_index-1) AS start, + least(start+$1*chop_index, "end") AS "end", + FROM path_segment + INNER JOIN segment_chop ON segment_chop.parent_id=path_segment.segment_id + WHERE segment_orientation='+' + UNION ALL + -- Reverse the order of chops in negatively oriented segments. + SELECT path_id, + id AS segment_id, + segment_orientation, + start+greatest(0, "end"-start-$1*chop_index) AS start, + "end"-$1*(chop_index-1) AS "end" + FROM path_segment + INNER JOIN segment_chop ON segment_chop.parent_id=path_segment.segment_id + WHERE segment_orientation='-') + INSERT INTO output_db.path_segment + SELECT * FROM path_segment_chop + -- Re-order similar to post-import.sql for optimal access. + ORDER BY path_id, start, "end"; diff --git a/domagi/chop.sql b/domagi/chop.sql deleted file mode 100644 index 571a8ef..0000000 --- a/domagi/chop.sql +++ /dev/null @@ -1,12 +0,0 @@ -WITH segment_chops AS ( - SELECT id, name, sequence, range(0, len(sequence), 3) AS starts - FROM segment), - segment_chop AS ( - SELECT id, name, sequence, 1 + unnest(starts) AS start, generate_subscripts(starts, 1)-1 AS chop_index - FROM segment_chops) - SELECT id, CASE WHEN chop_index=0 THEN name ELSE NULL END, array_slice(sequence, start, start + 3) - FROM segment_chop; - --- chop index is computed as ceil(a/b) = (a+b-1)//b --- SELECT id, name, unnest(range(0, len(sequence), 3)) AS starts, unnest(range(0, (len(sequence)+3-1)//3)) AS chop_index --- FROM segment; diff --git a/domagi/domagi.py b/domagi/domagi.py index adf3844..d3321b7 100644 --- a/domagi/domagi.py +++ b/domagi/domagi.py @@ -96,6 +96,36 @@ def build(gfa, db, threads): with connect_duckdb(db, threads) as con: con.execute(read_sql("post-import.sql")) +@main.command(short_help="Divide segments into smaller pieces") +@click.option("-i", "--db", "--idx", "con", + type=DuckDBParamType(), + required=True, + help="pangenome duckdb database") +@click.option("-o", "--out", "outfile", + type=click.Path(), + required=True, + help="path to output pangenome duckdb database") +@click.option("-c", "--chop-to", + type=click.INT, + metavar="N", + required=True, + help="divide segments longer than N") +@common_options +def chop(con, outfile, chop_to, threads): + set_duckdb_threads(con, threads) + with connect_duckdb(outfile, threads) as out_con: + out_con.execute(read_sql("schema.sql")) + con.execute(f"ATTACH '{outfile}' AS output_db (READ_WRITE)") + # At the moment, DuckDB only supports prepared parameters in the last + # statement. Hence, we have to split up the statements into separate execute + # calls. + con.execute(read_sql("chop-1.sql"), [chop_to]) + con.execute(read_sql("chop-2.sql"), [chop_to]) + con.execute(""" + DROP TABLE segment_chop; + DETACH output_db; + """) + @main.command(short_help="Crush runs of Ns") @click.option("-i", "--db", "--idx", "con", type=DuckDBParamType(), diff --git a/meson.build b/meson.build index 7b86e87..fe2e8ae 100644 --- a/meson.build +++ b/meson.build @@ -11,6 +11,8 @@ py.install_sources('domagi/domagi.py', install_data('domagi/matrix.sql', 'domagi/schema.sql', + 'domagi/chop-1.sql', + 'domagi/chop-2.sql', 'domagi/bed-depth.sql', 'domagi/overlap.sql', 'domagi/path-depth.sql', diff --git a/test-data/expected-output/test3-chop-1.gfa b/test-data/expected-output/test3-chop-1.gfa new file mode 100644 index 0000000..fa89dce --- /dev/null +++ b/test-data/expected-output/test3-chop-1.gfa @@ -0,0 +1,33 @@ +H VN:Z:1.0 +S 1 A +L 1 + 2 + 0M +S 2 G +L 2 + 3 + 0M +S 3 G +L 3 + 4 + 0M +S 4 A +L 4 + 5 + 0M +L 4 + 6 + 0M +S 5 A +L 5 + 8 + 0M +S 6 T +L 6 + 7 + 0M +L 6 - 8 + 0M +S 7 C +L 7 + 4 - 0M +L 7 + 8 + 0M +S 8 T +L 8 + 9 + 0M +S 9 C +L 9 + 10 + 0M +S 10 T +L 10 + 11 + 0M +S 11 C +L 11 + 12 + 0M +S 12 A +L 12 + 13 + 0M +S 13 G +L 13 + 14 + 0M +S 14 G +P 5+ 1+,2+,3+,4+,6+,7+,8+,9+,10+,11+,12+,13+,14+,8+,9+,10+,11+,12+,13+,14+ * +P 5- 1+,2+,3+,4+,7-,6-,8+,9+,10+,11+,12+,13+,14+ * diff --git a/test-data/expected-output/test3-chop-2.gfa b/test-data/expected-output/test3-chop-2.gfa new file mode 100644 index 0000000..bdc77d3 --- /dev/null +++ b/test-data/expected-output/test3-chop-2.gfa @@ -0,0 +1,21 @@ +H VN:Z:1.0 +S 1 AG +L 1 + 2 + 0M +S 2 GA +L 2 + 3 + 0M +L 2 + 4 + 0M +L 2 + 4 - 0M +S 3 A +L 3 + 5 + 0M +S 4 TC +L 4 - 5 + 0M +L 4 + 5 + 0M +S 5 TC +L 5 + 6 + 0M +S 6 TC +L 6 + 7 + 0M +S 7 AG +L 7 + 8 + 0M +S 8 G +P 5+ 1+,2+,4+,5+,6+,7+,8+,5+,6+,7+,8+ * +P 5- 1+,2+,4-,5+,6+,7+,8+ * diff --git a/test-data/expected-output/test3-chop-3.gfa b/test-data/expected-output/test3-chop-3.gfa new file mode 100644 index 0000000..0e16301 --- /dev/null +++ b/test-data/expected-output/test3-chop-3.gfa @@ -0,0 +1,19 @@ +H VN:Z:1.0 +S 1 AGG +L 1 + 2 + 0M +S 2 A +L 2 + 3 + 0M +L 2 + 4 + 0M +L 2 + 4 - 0M +S 3 A +L 3 + 5 + 0M +S 4 TC +L 4 - 5 + 0M +L 4 + 5 + 0M +S 5 TCT +L 5 + 6 + 0M +S 6 CAG +L 6 + 7 + 0M +S 7 G +P 5+ 1+,2+,4+,5+,6+,7+,5+,6+,7+ * +P 5- 1+,2+,4-,5+,6+,7+ * diff --git a/test-data/expected-output/test3-chop-4.gfa b/test-data/expected-output/test3-chop-4.gfa new file mode 100644 index 0000000..f043b82 --- /dev/null +++ b/test-data/expected-output/test3-chop-4.gfa @@ -0,0 +1,15 @@ +H VN:Z:1.0 +S 1 AGGA +L 1 + 2 + 0M +L 1 + 3 + 0M +L 1 + 3 - 0M +S 2 A +L 2 + 4 + 0M +S 3 TC +L 3 - 4 + 0M +L 3 + 4 + 0M +S 4 TCTC +L 4 + 5 + 0M +S 5 AGG +P 5+ 1+,3+,4+,5+,4+,5+ * +P 5- 1+,3-,4+,5+ * diff --git a/tests/test_domagi.py b/tests/test_domagi.py index 38c7c69..61b8e0b 100644 --- a/tests/test_domagi.py +++ b/tests/test_domagi.py @@ -52,6 +52,39 @@ def assert_gfa_equal(expected, actual): assert read_gfa_file(expected) == read_gfa_file(actual) +@pytest.mark.parametrize("test_data_file, chop_to, expected_output", + [(Path("test-data/test3.gfa"), + 1, + Path("test-data/expected-output/test3-chop-1.gfa")), + (Path("test-data/test3.gfa"), + 2, + Path("test-data/expected-output/test3-chop-2.gfa")), + (Path("test-data/test3.gfa"), + 3, + Path("test-data/expected-output/test3-chop-3.gfa")), + (Path("test-data/test3.gfa"), + 4, + Path("test-data/expected-output/test3-chop-4.gfa"))]) +def test_domagi_chop(tmp_path, test_data_file, chop_to, expected_output): + duckdb_path = tmp_path / f"{test_data_file.stem}.db" + output_duckdb_path = tmp_path / f"{test_data_file.stem}-output.db" + runner = CliRunner() + result = runner.invoke(main, ["build", + "--gfa", test_data_file, + "--out", duckdb_path]) + assert result.exit_code == 0 + result = runner.invoke(main, ["chop", + "--db", duckdb_path, + "--chop-to", chop_to, + "--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("test_data_file, expected_output", [(Path("test-data/test-crush.gfa"), Path("test-data/expected-output/test-crush.gfa"))]) -- cgit 1.4.1