From 336ae557bf2bcd43a64836c4fd900f3de406d139 Mon Sep 17 00:00:00 2001 From: Charlie Vaske Date: Wed, 19 Aug 2026 13:11:23 -0700 Subject: [PATCH 1/2] [OC-832] Fix dropped sample/mapping rows on merge; tighten mergesqlite checks - fix mergesqlite() dropping sample and mapping rows for any variant that already exists in the merged output (uid_dic was never populated for the skip-duplicate-variant case, so later sample/mapping inserts for that uid were silently skipped) - extend the pre-merge consistency check beyond variant/gene_header column names to also cover sample_header/mapping_header columns and (name, version) pairs for variant_annotator/gene_annotator across all input dbs; error out naming the mismatch before doing any merge work - recompute "Number of unique input variants" and "Result modified at" in the merged db's info table instead of leaving db 1's stale values - detect colliding sample_id values across input dbs and error out by default - add an opt-in `path:label` argument suffix to rename a file's sample_ids to `label__sample_id` to resolve a collision - add tests/test_mergesqlite.py covering all of the above --- cravat/cravat_util.py | 137 ++++++++++--- tests/test_mergesqlite.py | 399 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 510 insertions(+), 26 deletions(-) create mode 100644 tests/test_mergesqlite.py diff --git a/cravat/cravat_util.py b/cravat/cravat_util.py index b4c257d1..23f6f869 100644 --- a/cravat/cravat_util.py +++ b/cravat/cravat_util.py @@ -761,32 +761,92 @@ def showsqliteinfo(args): c.close() conn.close() +def mergesqlite_check_info(dbpath): + """Collects the header columns, annotator module versions, and sample + ids used by a result db, for the pre-merge consistency checks in + mergesqlite().""" + conn = sqlite3.connect(dbpath) + c = conn.cursor() + info = {} + for table in ["variant", "gene", "sample", "mapping"]: + c.execute(f'select col_name from {table}_header') + info[table] = sorted([r[0] for r in c.fetchall()]) + for annot_table, sql_table in [("variant_annotators", "variant_annotator"), + ("gene_annotators", "gene_annotator")]: + c.execute(f'select name, version from {sql_table}') + info[annot_table] = {r[0]: r[1] for r in c.fetchall()} + c.execute('select distinct base__sample_id from sample') + info["sample_ids"] = sorted({r[0] for r in c.fetchall()}) + c.close() + conn.close() + return info + +def mergesqlite_parse_path_arg(raw): + """Parses a `path` or `path:label` positional arg for mergesqlite. + A label is only recognized when the part before the last ':' exists + as a file and the raw string as a whole does not (so plain paths + with no label, including Windows drive letters, pass through as-is). + Returns (path, label), label is None when no label was given.""" + raw = str(raw) + if ':' in raw and not os.path.exists(raw): + maybe_path, maybe_label = raw.rsplit(':', 1) + if maybe_label and os.path.exists(maybe_path): + return maybe_path, maybe_label + return raw, None + # For now, only jobs with same annotators are allowed. def mergesqlite(args): - dbpaths = args.path - if len(dbpaths) < 2: + raw_paths = args.path + if len(raw_paths) < 2: exit("Multiple sqlite file paths should be given") + dbpaths = [] + labels = {} + for raw in raw_paths: + dbpath, label = mergesqlite_parse_path_arg(raw) + dbpaths.append(dbpath) + if label is not None: + labels[dbpath] = label outpath = args.outpath if outpath.endswith('.sqlite') == False: outpath = outpath + '.sqlite' - # Checks columns being the same. - conn = sqlite3.connect(dbpaths[0]) - c = conn.cursor() - c.execute('select col_name from variant_header') - v_cols = sorted([r[0] for r in c.fetchall()]) - c.execute('select col_name from gene_header') - g_cols = sorted([r[0] for r in c.fetchall()]) - c.close() - conn.close() + # Checks columns and annotator modules being the same, and that no + # sample_id collides across inputs (after any :label rename). + all_info = {dbpath: mergesqlite_check_info(dbpath) for dbpath in dbpaths} + base_info = all_info[dbpaths[0]] for dbpath in dbpaths[1:]: - conn = sqlite3.connect(dbpath) - c = conn.cursor() - c.execute('select col_name from variant_header') - if v_cols != sorted([r[0] for r in c.fetchall()]): - exit("Annotation columns mismatch (variant table)") - c.execute('select col_name from gene_header') - if g_cols != sorted([r[0] for r in c.fetchall()]): - exit("Annotation columns mismatch (gene table)") + info = all_info[dbpath] + for table in ["variant", "gene", "sample", "mapping"]: + if base_info[table] != info[table]: + exit( + f'Annotation columns mismatch ({table} table) between ' + f'{dbpaths[0]} and {dbpath}' + ) + for annot_table in ["variant_annotators", "gene_annotators"]: + base_annots = base_info[annot_table] + annots = info[annot_table] + for name in sorted(set(base_annots) | set(annots)): + base_version = base_annots.get(name) + version = annots.get(name) + if base_version != version: + exit( + f'Annotator module mismatch ({annot_table.replace("_annotators", "")} ' + f'annotator "{name}"): version {base_version} in {dbpaths[0]} vs ' + f'version {version} in {dbpath}' + ) + sample_id_sources = {} + for dbpath in dbpaths: + label = labels.get(dbpath) + for sid in all_info[dbpath]["sample_ids"]: + eff_sid = f'{label}__{sid}' if label else sid + sample_id_sources.setdefault(eff_sid, []).append(dbpath) + collisions = {sid: paths for sid, paths in sample_id_sources.items() if len(paths) > 1} + if collisions: + lines = [f' "{sid}": {", ".join(paths)}' for sid, paths in sorted(collisions.items())] + exit( + "Sample ID collision(s) across input files. Give the colliding " + "file(s) a path:label suffix to disambiguate (e.g. " + "job1.sqlite:cohortA):\n" + "\n".join(lines) + ) # Copies the first db. print(f'Copying {dbpaths[0]} to {outpath}...') shutil.copy(dbpaths[0], outpath) @@ -805,12 +865,19 @@ def mergesqlite(args): outc.execute('select col_name from sample_header order by rowid') cols = [r[0] for r in outc.fetchall()] s_uid_colno = cols.index('base__uid') + s_sampleid_colno = cols.index('base__sample_id') outc.execute('select col_name from mapping_header order by rowid') cols = [r[0] for r in outc.fetchall()] m_uid_colno = cols.index('base__uid') m_fileno_colno = cols.index('base__fileno') outc.execute('select max(base__uid) from variant') new_uid = outc.fetchone()[0] + 1 + # Renames db 1's own sample_ids if it was given a :label. + if dbpaths[0] in labels: + outc.execute( + 'update sample set base__sample_id = ? || base__sample_id', + (f'{labels[dbpaths[0]]}__',) + ) # Input paths outc.execute('select colkey, colval from info where colkey="_input_paths"') input_paths = json.loads(outc.fetchone()[1].replace("'", '"')) @@ -818,13 +885,14 @@ def mergesqlite(args): rev_input_paths = {} for fileno, filepath in input_paths.items(): rev_input_paths[filepath] = fileno - # Makes initial hugo and variant id lists. + # Makes initial hugo and variant id -> uid lists. outc.execute('select base__hugo from gene') genes = {r[0] for r in outc.fetchall()} - outc.execute('select base__chrom, base__pos, base__ref_base, base__alt_base from variant') - variants = {variant_id(r[0], r[1], r[2] ,r[3]) for r in outc.fetchall()} + outc.execute('select base__uid, base__chrom, base__pos, base__ref_base, base__alt_base from variant') + vid_to_uid = {variant_id(r[1], r[2], r[3], r[4]): r[0] for r in outc.fetchall()} for dbpath in dbpaths[1:]: print(f'Merging {dbpath}...') + label = labels.get(dbpath) conn = sqlite3.connect(dbpath) c = conn.cursor() # Gene @@ -841,16 +909,21 @@ def mergesqlite(args): c.execute('select * from variant order by rowid') for r in c.fetchall(): vid = variant_id(r[v_chrom_colno], r[v_pos_colno], r[v_ref_colno], r[v_alt_colno]) - if vid in variants: - continue old_uid = r[0] + if vid in vid_to_uid: + # Variant already present in the merged output (annotation + # is identical, so the redundant insert is skipped) - but + # the uid mapping still needs recording so this variant's + # sample/mapping rows get merged in below. + uid_dic[old_uid] = vid_to_uid[vid] + continue r = list(r) r[0] = new_uid uid_dic[old_uid] = new_uid + vid_to_uid[vid] = new_uid new_uid += 1 q = f'insert into variant values ({",".join(["?" for v in range(len(r))])})' outc.execute(q, r) - variants.add(vid) # Sample c.execute('select * from sample order by rowid') for r in c.fetchall(): @@ -859,6 +932,8 @@ def mergesqlite(args): new_uid = uid_dic[uid] r = list(r) r[s_uid_colno] = new_uid + if label: + r[s_sampleid_colno] = f'{label}__{r[s_sampleid_colno]}' q = f'insert into sample values ({",".join(["?" for v in range(len(r))])})' outc.execute(q, r) # File numbers @@ -887,6 +962,13 @@ def mergesqlite(args): q = 'update info set colval=? where colkey="Input file name"' v = ';'.join([input_paths[str(v)] for v in sorted(input_paths.keys(), key=lambda v: int(v))]) outc.execute(q, [v]) + outc.execute('select count(*) from variant') + n_variants = outc.fetchone()[0] + q = 'update info set colval=? where colkey="Number of unique input variants"' + outc.execute(q, [str(n_variants)]) + modified = datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S') + q = 'update info set colval=? where colkey="Result modified at"' + outc.execute(q, [modified]) outconn.commit() @@ -1193,7 +1275,10 @@ def jobtopackage(args): parser_mergesqlite = subparsers.add_parser( "mergesqlite", help="Merge SQLite result files" ) -parser_mergesqlite.add_argument("path", nargs='+', help="Path to result database", type=Path) +parser_mergesqlite.add_argument("path", nargs='+', + help="Path to result database. Optionally 'path:label' to rename that " + "db's sample_ids to 'label__sample_id' on merge, to resolve " + "sample_id collisions with other input dbs.") parser_mergesqlite.add_argument("-o", dest="outpath", required=True, help="Output SQLite file path") parser_mergesqlite.set_defaults(func=mergesqlite) diff --git a/tests/test_mergesqlite.py b/tests/test_mergesqlite.py new file mode 100644 index 00000000..99a73279 --- /dev/null +++ b/tests/test_mergesqlite.py @@ -0,0 +1,399 @@ +import json +import os +import shutil +import sqlite3 +import tempfile +import unittest +from types import SimpleNamespace + +from cravat.cravat_util import mergesqlite + +# Column shapes mirror what real oc-produced result dbs carry: a handful of +# "base__*" key columns plus one arbitrary annotator column per table, so +# the header/annotator consistency checks in mergesqlite() have something +# real to compare. +VARIANT_COLS = [ + ("base__uid", "integer"), + ("base__chrom", "text"), + ("base__pos", "integer"), + ("base__ref_base", "text"), + ("base__alt_base", "text"), + ("test__score", "real"), +] +GENE_COLS = [("base__hugo", "text"), ("test__gscore", "real")] +SAMPLE_COLS = [ + ("base__uid", "integer"), + ("base__sample_id", "text"), + ("base__zygosity", "text"), +] +MAPPING_COLS = [ + ("base__uid", "integer"), + ("base__fileno", "integer"), + ("base__transcript", "text"), +] + + +def _col_def(col_name): + return json.dumps({"title": col_name, "type": "string"}) + + +def build_db( + path, + variants, + samples, + mappings, + genes, + input_paths, + variant_annotator_version="1.0.0", + gene_annotator_version="1.0.0", + variant_cols=VARIANT_COLS, + gene_cols=GENE_COLS, + sample_cols=SAMPLE_COLS, + mapping_cols=MAPPING_COLS, +): + """Builds a minimal but real-shaped cravat result sqlite db. + + variants: rows of (uid, chrom, pos, ref, alt, score) + samples: rows of (uid, sample_id, zygosity) + mappings: rows of (uid, fileno, transcript) + genes: rows of (hugo, gscore) + input_paths: {str(fileno): filepath} + """ + conn = sqlite3.connect(path) + c = conn.cursor() + + for table, cols, rows in [ + ("variant", variant_cols, variants), + ("gene", gene_cols, genes), + ("sample", sample_cols, samples), + ("mapping", mapping_cols, mappings), + ]: + col_sql = ", ".join(f"{name} {sqltype}" for name, sqltype in cols) + c.execute(f"create table {table} ({col_sql})") + placeholders = ",".join("?" * len(cols)) + c.executemany(f"insert into {table} values ({placeholders})", rows) + + for header_table, cols in [ + ("variant_header", variant_cols), + ("gene_header", gene_cols), + ("sample_header", sample_cols), + ("mapping_header", mapping_cols), + ]: + c.execute(f"create table {header_table} (col_name text, col_def text)") + c.executemany( + f"insert into {header_table} values (?, ?)", + [(name, _col_def(name)) for name, _ in cols], + ) + + c.execute("create table variant_annotator (name text, displayname text, version text)") + c.execute( + "insert into variant_annotator values (?, ?, ?)", + ("test", "Test Annotator", variant_annotator_version), + ) + c.execute("create table gene_annotator (name text, displayname text, version text)") + c.execute( + "insert into gene_annotator values (?, ?, ?)", + ("test_gene", "Test Gene Annotator", gene_annotator_version), + ) + + c.execute("create table info (colkey text primary key, colval text)") + input_paths_str = json.dumps(input_paths).replace('"', "'") + c.executemany( + "insert into info values (?, ?)", + [ + ("_input_paths", input_paths_str), + ("Input file name", ";".join(input_paths.values())), + ("Number of unique input variants", str(len(variants))), + ("Result created at", "2026-01-01 00:00:00"), + ("Result modified at", "2026-01-01 00:00:00"), + ], + ) + + conn.commit() + conn.close() + + +class MergeSqliteTestBase(unittest.TestCase): + def setUp(self): + self.tmpdir = tempfile.mkdtemp() + self.addCleanup(shutil.rmtree, self.tmpdir, ignore_errors=True) + self.db1 = os.path.join(self.tmpdir, "db1.sqlite") + self.db2 = os.path.join(self.tmpdir, "db2.sqlite") + self.outpath = os.path.join(self.tmpdir, "merged.sqlite") + + def run_merge(self, paths): + args = SimpleNamespace(path=paths, outpath=self.outpath) + mergesqlite(args) + + def query(self, sql, params=()): + conn = sqlite3.connect(self.outpath) + c = conn.cursor() + c.execute(sql, params) + rows = c.fetchall() + conn.close() + return rows + + +class TestSharedVariant(MergeSqliteTestBase): + def test_shared_variant_merges_sample_and_mapping_rows(self): + build_db( + self.db1, + variants=[(1, "chr1", 100, "A", "T", 0.5)], + samples=[(1, "sample1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + build_db( + self.db2, + variants=[(1, "chr1", 100, "A", "T", 0.5)], + samples=[(1, "sample2", "hom")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db2.vcf"}, + ) + + self.run_merge([self.db1, self.db2]) + + variant_rows = self.query("select base__uid from variant") + self.assertEqual(len(variant_rows), 1, "shared variant must not be duplicated") + uid = variant_rows[0][0] + + sample_rows = self.query( + "select base__uid, base__sample_id from sample order by base__sample_id" + ) + self.assertEqual( + sample_rows, [(uid, "sample1"), (uid, "sample2")], + "both samples for the shared variant must be present, under the same uid", + ) + + mapping_rows = self.query("select base__uid from mapping") + self.assertEqual( + len(mapping_rows), 2, + "both mapping rows for the shared variant must be present", + ) + self.assertTrue(all(r[0] == uid for r in mapping_rows)) + + +class TestUniqueVariants(MergeSqliteTestBase): + def test_unique_variants_get_noncolliding_uids(self): + build_db( + self.db1, + variants=[(1, "chr1", 100, "A", "T", 0.5)], + samples=[(1, "sample1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + build_db( + self.db2, + variants=[(1, "chr2", 200, "C", "G", 0.7)], + samples=[(1, "sample2", "hom")], + mappings=[(1, 0, "NM_002")], + genes=[("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + ) + + self.run_merge([self.db1, self.db2]) + + variant_rows = self.query( + "select base__uid, base__chrom from variant order by base__chrom" + ) + self.assertEqual(len(variant_rows), 2) + uids = [r[0] for r in variant_rows] + self.assertEqual(len(set(uids)), 2, "uids must not collide") + + # Every sample/mapping row must reference one of the real merged uids. + sample_uids = {r[0] for r in self.query("select base__uid from sample")} + mapping_uids = {r[0] for r in self.query("select base__uid from mapping")} + self.assertEqual(sample_uids, set(uids)) + self.assertEqual(mapping_uids, set(uids)) + + +class TestConsistencyChecks(MergeSqliteTestBase): + def test_mismatched_variant_annotator_version_errors_no_output(self): + build_db( + self.db1, + variants=[(1, "chr1", 100, "A", "T", 0.5)], + samples=[(1, "sample1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + variant_annotator_version="1.0.0", + ) + build_db( + self.db2, + variants=[(1, "chr2", 200, "C", "G", 0.7)], + samples=[(1, "sample2", "hom")], + mappings=[(1, 0, "NM_002")], + genes=[("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + variant_annotator_version="2.0.0", + ) + + with self.assertRaises(SystemExit): + self.run_merge([self.db1, self.db2]) + self.assertFalse(os.path.exists(self.outpath)) + + def test_mismatched_sample_header_columns_errors_no_output(self): + build_db( + self.db1, + variants=[(1, "chr1", 100, "A", "T", 0.5)], + samples=[(1, "sample1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + extra_sample_cols = SAMPLE_COLS + [("base__extra", "text")] + build_db( + self.db2, + variants=[(1, "chr2", 200, "C", "G", 0.7)], + samples=[(1, "sample2", "hom", "x")], + mappings=[(1, 0, "NM_002")], + genes=[("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + sample_cols=extra_sample_cols, + ) + + with self.assertRaises(SystemExit): + self.run_merge([self.db1, self.db2]) + self.assertFalse(os.path.exists(self.outpath)) + + def test_mismatched_mapping_header_columns_errors_no_output(self): + build_db( + self.db1, + variants=[(1, "chr1", 100, "A", "T", 0.5)], + samples=[(1, "sample1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + extra_mapping_cols = MAPPING_COLS + [("base__extra", "text")] + build_db( + self.db2, + variants=[(1, "chr2", 200, "C", "G", 0.7)], + samples=[(1, "sample2", "hom")], + mappings=[(1, 0, "NM_002", "x")], + genes=[("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + mapping_cols=extra_mapping_cols, + ) + + with self.assertRaises(SystemExit): + self.run_merge([self.db1, self.db2]) + self.assertFalse(os.path.exists(self.outpath)) + + +class TestMergedInfo(MergeSqliteTestBase): + def test_merged_info_has_recomputed_variant_count_and_input_paths(self): + build_db( + self.db1, + variants=[(1, "chr1", 100, "A", "T", 0.5)], + samples=[(1, "sample1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + build_db( + self.db2, + variants=[ + (1, "chr1", 100, "A", "T", 0.5), # shared with db1 + (2, "chr2", 200, "C", "G", 0.7), # unique to db2 + ], + samples=[(1, "sample2", "hom"), (2, "sample2", "hom")], + mappings=[(1, 0, "NM_001"), (2, 0, "NM_002")], + genes=[("GENE1", 0.9), ("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + ) + + self.run_merge([self.db1, self.db2]) + + n_variant_rows = self.query("select count(*) from variant")[0][0] + self.assertEqual(n_variant_rows, 2) + + colval = self.query( + 'select colval from info where colkey="Number of unique input variants"' + )[0][0] + self.assertEqual(int(colval), n_variant_rows) + + input_paths_raw = self.query( + 'select colval from info where colkey="_input_paths"' + )[0][0] + input_paths = json.loads(input_paths_raw.replace("'", '"')) + self.assertEqual(set(input_paths.values()), {"/in/db1.vcf", "/in/db2.vcf"}) + + +class TestSampleIdCollision(MergeSqliteTestBase): + def test_collision_no_label_errors_no_output(self): + build_db( + self.db1, + variants=[(1, "chr1", 100, "A", "T", 0.5)], + samples=[(1, "SAMPLE1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + build_db( + self.db2, + variants=[(1, "chr2", 200, "C", "G", 0.7)], + samples=[(1, "SAMPLE1", "hom")], + mappings=[(1, 0, "NM_002")], + genes=[("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + ) + + with self.assertRaises(SystemExit): + self.run_merge([self.db1, self.db2]) + self.assertFalse(os.path.exists(self.outpath)) + + def test_collision_with_label_renames_only_labeled_db(self): + build_db( + self.db1, + variants=[(1, "chr1", 100, "A", "T", 0.5)], + samples=[(1, "SAMPLE1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + build_db( + self.db2, + variants=[(1, "chr2", 200, "C", "G", 0.7)], + samples=[(1, "SAMPLE1", "hom")], + mappings=[(1, 0, "NM_002")], + genes=[("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + ) + + self.run_merge([self.db1, f"{self.db2}:cohortB"]) + + sample_ids = {r[0] for r in self.query("select base__sample_id from sample")} + self.assertEqual(sample_ids, {"SAMPLE1", "cohortB__SAMPLE1"}) + + def test_collision_with_label_on_first_db_still_renames(self): + # Exercises the rename applied to the already-copied base db, + # rather than to a db merged in later. + build_db( + self.db1, + variants=[(1, "chr1", 100, "A", "T", 0.5)], + samples=[(1, "SAMPLE1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + build_db( + self.db2, + variants=[(1, "chr2", 200, "C", "G", 0.7)], + samples=[(1, "SAMPLE1", "hom")], + mappings=[(1, 0, "NM_002")], + genes=[("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + ) + + self.run_merge([f"{self.db1}:cohortA", self.db2]) + + sample_ids = {r[0] for r in self.query("select base__sample_id from sample")} + self.assertEqual(sample_ids, {"cohortA__SAMPLE1", "SAMPLE1"}) + + +if __name__ == "__main__": + unittest.main() From 0703d35d4068a14db78fd58ca8f1b0db57616833 Mon Sep 17 00:00:00 2001 From: Charlie Vaske Date: Thu, 20 Aug 2026 09:24:49 -0700 Subject: [PATCH 2/2] Fix data merge bugs in mergesqlite() Four correctness bugs in mergesqlite() could silently drop, duplicate, or misattribute merged sample/mapping rows: - new_uid running counter was clobbered by the Sample/Mapping blocks reusing it as a scratch variable, causing duplicate base__uid values and misattributed sample rows when merging 3+ dbs. - variant_id() concatenated chrom+pos+ref+alt with no delimiter, so distinct variants (e.g. chr=1/pos=234 vs chr=12/pos=34) could collide and get silently merged onto the wrong variant. Now delimited with ':'. - fileno_dic only recorded entries for newly-seen input filepaths, raising KeyError when a later db repeated an earlier db's input path; now maps repeated filepaths onto their existing fileno. - labels was a dict keyed by dbpath, so passing the same physical file twice with different :label suffixes silently dropped the first label; now a list parallel to dbpaths. - mergesqlite_check_info()'s column consistency check now also checks order of columns, which matehes the merging logic used later on. Added regression tests for all five, each confirmed to fail against the prior behavior before the fix. --- cravat/cravat_util.py | 41 ++++--- tests/test_mergesqlite.py | 231 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 256 insertions(+), 16 deletions(-) diff --git a/cravat/cravat_util.py b/cravat/cravat_util.py index 23f6f869..5dff8bc2 100644 --- a/cravat/cravat_util.py +++ b/cravat/cravat_util.py @@ -723,7 +723,7 @@ def result2gui(args): def variant_id(chrom, pos, ref, alt): - return chrom + str(pos) + ref + alt + return chrom + ':' + str(pos) + ':' + ref + ':' + alt def showsqliteinfo(args): dbpaths = args.paths @@ -769,8 +769,12 @@ def mergesqlite_check_info(dbpath): c = conn.cursor() info = {} for table in ["variant", "gene", "sample", "mapping"]: - c.execute(f'select col_name from {table}_header') - info[table] = sorted([r[0] for r in c.fetchall()]) + # Order matters here (no sorted()): the merge loop in mergesqlite() + # reads and writes rows positionally, in this same rowid order, so + # two dbs with identical column names but a different physical + # order must be treated as a mismatch, not silently accepted. + c.execute(f'select col_name from {table}_header order by rowid') + info[table] = [r[0] for r in c.fetchall()] for annot_table, sql_table in [("variant_annotators", "variant_annotator"), ("gene_annotators", "gene_annotator")]: c.execute(f'select name, version from {sql_table}') @@ -800,12 +804,14 @@ def mergesqlite(args): if len(raw_paths) < 2: exit("Multiple sqlite file paths should be given") dbpaths = [] - labels = {} + # Parallel to dbpaths (not a dict keyed by dbpath) so that passing the + # same physical file twice with two different :label suffixes keeps + # both labels instead of the second overwriting the first. + labels = [] for raw in raw_paths: dbpath, label = mergesqlite_parse_path_arg(raw) dbpaths.append(dbpath) - if label is not None: - labels[dbpath] = label + labels.append(label) outpath = args.outpath if outpath.endswith('.sqlite') == False: outpath = outpath + '.sqlite' @@ -834,8 +840,7 @@ def mergesqlite(args): f'version {version} in {dbpath}' ) sample_id_sources = {} - for dbpath in dbpaths: - label = labels.get(dbpath) + for dbpath, label in zip(dbpaths, labels): for sid in all_info[dbpath]["sample_ids"]: eff_sid = f'{label}__{sid}' if label else sid sample_id_sources.setdefault(eff_sid, []).append(dbpath) @@ -873,10 +878,10 @@ def mergesqlite(args): outc.execute('select max(base__uid) from variant') new_uid = outc.fetchone()[0] + 1 # Renames db 1's own sample_ids if it was given a :label. - if dbpaths[0] in labels: + if labels[0]: outc.execute( 'update sample set base__sample_id = ? || base__sample_id', - (f'{labels[dbpaths[0]]}__',) + (f'{labels[0]}__',) ) # Input paths outc.execute('select colkey, colval from info where colkey="_input_paths"') @@ -890,9 +895,8 @@ def mergesqlite(args): genes = {r[0] for r in outc.fetchall()} outc.execute('select base__uid, base__chrom, base__pos, base__ref_base, base__alt_base from variant') vid_to_uid = {variant_id(r[1], r[2], r[3], r[4]): r[0] for r in outc.fetchall()} - for dbpath in dbpaths[1:]: + for dbpath, label in zip(dbpaths[1:], labels[1:]): print(f'Merging {dbpath}...') - label = labels.get(dbpath) conn = sqlite3.connect(dbpath) c = conn.cursor() # Gene @@ -929,9 +933,9 @@ def mergesqlite(args): for r in c.fetchall(): uid = r[s_uid_colno] if uid in uid_dic: - new_uid = uid_dic[uid] + mapped_uid = uid_dic[uid] r = list(r) - r[s_uid_colno] = new_uid + r[s_uid_colno] = mapped_uid if label: r[s_sampleid_colno] = f'{label}__{r[s_sampleid_colno]}' q = f'insert into sample values ({",".join(["?" for v in range(len(r))])})' @@ -946,14 +950,19 @@ def mergesqlite(args): rev_input_paths[filepath] = str(new_fileno) fileno_dic[int(fileno)] = new_fileno new_fileno += 1 + else: + # This db's input filepath was already contributed by an + # earlier db (or is db 1's own) - map its fileno onto the + # fileno already assigned to that filepath. + fileno_dic[int(fileno)] = int(rev_input_paths[filepath]) # Mapping c.execute('select * from mapping order by rowid') for r in c.fetchall(): uid = r[m_uid_colno] if uid in uid_dic: - new_uid = uid_dic[uid] + mapped_uid = uid_dic[uid] r = list(r) - r[m_uid_colno] = new_uid + r[m_uid_colno] = mapped_uid r[m_fileno_colno] = fileno_dic[r[m_fileno_colno]] q = f'insert into mapping values ({",".join(["?" for v in range(len(r))])})' outc.execute(q, r) diff --git a/tests/test_mergesqlite.py b/tests/test_mergesqlite.py index 99a73279..d9fc2825 100644 --- a/tests/test_mergesqlite.py +++ b/tests/test_mergesqlite.py @@ -395,5 +395,236 @@ def test_collision_with_label_on_first_db_still_renames(self): self.assertEqual(sample_ids, {"cohortA__SAMPLE1", "SAMPLE1"}) +class TestNewUidCounterClobbered(MergeSqliteTestBase): + def test_third_dbs_new_variant_does_not_reuse_second_dbs_new_uid(self): + # Regression test for the `new_uid` running-counter variable being + # clobbered: the Sample/Mapping blocks reuse `new_uid` as a loop + # scratch variable (`new_uid = uid_dic[uid]`), stomping on the + # counter the Variant block relies on to hand out fresh, unused + # base__uid values across every remaining dbpath. + # + # db1 contributes one variant (uid 1). + # db2 contributes that *same* variant (shared with db1) plus one + # brand-new variant, and has sample/mapping rows for both - + # processing the shared-variant sample row after the new-variant + # variant row is what clobbers `new_uid` back down. + # db3 then contributes one more brand-new variant. With a correct + # running counter it must get a uid distinct from every uid used + # so far; the bug hands it db2's new variant's uid instead. + db3 = os.path.join(self.tmpdir, "db3.sqlite") + + build_db( + self.db1, + variants=[(1, "1", 100, "A", "T", 0.5)], + samples=[(1, "sample1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + build_db( + self.db2, + variants=[ + (1, "1", 100, "A", "T", 0.5), # shared with db1 + (2, "2", 200, "C", "G", 0.7), # new + ], + samples=[(1, "sample2", "het"), (2, "sample3", "het")], + mappings=[(1, 0, "NM_001"), (2, 0, "NM_002")], + genes=[("GENE1", 0.9), ("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + ) + build_db( + db3, + variants=[(1, "3", 300, "G", "A", 0.9)], # new + samples=[(1, "sample4", "het")], + mappings=[(1, 0, "NM_003")], + genes=[("GENE3", 0.1)], + input_paths={"0": "/in/db3.vcf"}, + ) + + self.run_merge([self.db1, self.db2, db3]) + + variant_rows = self.query("select base__uid, base__chrom from variant") + self.assertEqual( + len(variant_rows), 3, "all three distinct variants must be present" + ) + uids = [r[0] for r in variant_rows] + self.assertEqual( + len(set(uids)), 3, + f"every distinct variant must get its own base__uid, got {uids}", + ) + + uid_by_chrom = {chrom: uid for uid, chrom in variant_rows} + sample4_uid = self.query( + 'select base__uid from sample where base__sample_id="sample4"' + )[0][0] + self.assertEqual( + sample4_uid, uid_by_chrom["3"], + "sample4 must be attached to db3's own new variant, not an " + "earlier variant whose uid got reused", + ) + + +class TestVariantIdDelimiterCollision(MergeSqliteTestBase): + def test_distinct_variants_with_colliding_concatenated_id_stay_distinct(self): + # Regression test for variant_id()'s undelimited concatenation + # (chrom + str(pos) + ref + alt): chrom="1"/pos=234 and + # chrom="12"/pos=34 both concatenate to "1234", so with the same + # ref/alt these two genuinely different variants collide onto the + # same id string. + build_db( + self.db1, + variants=[(1, "1", 234, "A", "G", 0.5)], + samples=[(1, "sample1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + build_db( + self.db2, + variants=[(1, "12", 34, "A", "G", 0.7)], + samples=[(1, "sample2", "hom")], + mappings=[(1, 0, "NM_002")], + genes=[("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + ) + + self.run_merge([self.db1, self.db2]) + + variant_rows = self.query( + "select base__uid, base__chrom, base__pos from variant" + ) + self.assertEqual( + len(variant_rows), 2, + "chr1:234 and chr12:34 are distinct variants and must both be " + "present, despite variant_id() concatenating them the same way", + ) + + uid_by_locus = {(chrom, pos): uid for uid, chrom, pos in variant_rows} + sample_uids = dict( + self.query("select base__sample_id, base__uid from sample") + ) + self.assertEqual( + sample_uids["sample1"], uid_by_locus[("1", 234)], + "sample1 belongs to the chr1:234 variant", + ) + self.assertEqual( + sample_uids["sample2"], uid_by_locus[("12", 34)], + "sample2 belongs to the chr12:34 variant, not chr1:234's uid", + ) + + +class TestRepeatedInputFileAcrossDbs(MergeSqliteTestBase): + def test_merging_dbs_with_shared_input_filepath_does_not_crash(self): + # Regression test for fileno_dic only recording an entry for a + # source db's fileno when that db's input filepath is *new* to the + # merge, while the Mapping block unconditionally looks up every + # row's fileno in fileno_dic - so a later db whose input filepath + # was already contributed by an earlier db raises a KeyError. + build_db( + self.db1, + variants=[(1, "1", 100, "A", "T", 0.5)], + samples=[(1, "sample1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/shared.vcf"}, + ) + build_db( + self.db2, + variants=[(1, "2", 200, "C", "G", 0.7)], + samples=[(1, "sample2", "hom")], + mappings=[(1, 0, "NM_002")], + genes=[("GENE2", 0.3)], + input_paths={"0": "/in/shared.vcf"}, # same filepath as db1 + ) + + self.run_merge([self.db1, self.db2]) + + mapping_rows = self.query( + "select base__uid, base__fileno from mapping order by base__uid" + ) + self.assertEqual(len(mapping_rows), 2, "both mapping rows must be merged in") + filenos = {fileno for _, fileno in mapping_rows} + self.assertEqual( + len(filenos), 1, + "both mapping rows point at the same shared input file and " + "must resolve to the same merged fileno", + ) + + +class TestDuplicateDbpathLabels(MergeSqliteTestBase): + def test_same_physical_db_merged_twice_with_different_labels(self): + # Regression test for `labels` being a dict keyed by the resolved + # dbpath string: passing the same physical file twice with two + # different :label suffixes silently drops the first label, since + # the second assignment overwrites labels[dbpath]. + build_db( + self.db1, + variants=[(1, "1", 100, "A", "T", 0.5)], + samples=[(1, "sample1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + build_db( + self.db2, + variants=[(1, "2", 200, "C", "G", 0.7)], + samples=[(1, "SAMPLE_X", "hom")], + mappings=[(1, 0, "NM_002")], + genes=[("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + ) + + self.run_merge([self.db1, f"{self.db2}:cohortA", f"{self.db2}:cohortB"]) + + sample_ids = {r[0] for r in self.query("select base__sample_id from sample")} + self.assertEqual( + sample_ids, + {"sample1", "cohortA__SAMPLE_X", "cohortB__SAMPLE_X"}, + "each :label suffix on the same physical db must be honored " + "independently, not collapsed to the last one seen", + ) + + +class TestColumnOrderMismatch(MergeSqliteTestBase): + def test_same_column_names_different_order_errors_no_output(self): + # Regression test for mergesqlite_check_info() sorting column names + # before comparing them: two dbs with the same column *set* but a + # different physical *order* used to pass the check even though + # the merge loop reads/writes rows positionally, using db1's + # column order for every db - a same-named-but-reordered column + # (e.g. ref_base/alt_base swapped) would silently merge wrong + # values into the output with no error at all. + build_db( + self.db1, + variants=[(1, "chr1", 100, "A", "T", 0.5)], + samples=[(1, "sample1", "het")], + mappings=[(1, 0, "NM_001")], + genes=[("GENE1", 0.9)], + input_paths={"0": "/in/db1.vcf"}, + ) + # Same column names as VARIANT_COLS, but ref_base/alt_base swapped. + reordered_variant_cols = [ + ("base__uid", "integer"), + ("base__chrom", "text"), + ("base__pos", "integer"), + ("base__alt_base", "text"), + ("base__ref_base", "text"), + ("test__score", "real"), + ] + build_db( + self.db2, + variants=[(1, "chr2", 200, "G", "C", 0.7)], # matches reordered cols + samples=[(1, "sample2", "hom")], + mappings=[(1, 0, "NM_002")], + genes=[("GENE2", 0.3)], + input_paths={"0": "/in/db2.vcf"}, + variant_cols=reordered_variant_cols, + ) + + with self.assertRaises(SystemExit): + self.run_merge([self.db1, self.db2]) + self.assertFalse(os.path.exists(self.outpath)) + + if __name__ == "__main__": unittest.main()