Merge pull request #1027 from Sharon-codes/issue-919-cnv-variant-query

test(cnv): stabilise simulated high-variance sampling

Author: jonbrenas

tests/anoph/conftest.py

@@ -334,7 +334,7 @@ def simulate_exons(
         # keep things simple for now.
         if strand == "-":
             # Take exons in reverse order.
-            exons == exons[::-1]
+            exons = exons[::-1]
         for exon_ix, exon in enumerate(exons):
             first_exon = exon_ix == 0
             last_exon = exon_ix == len(exons) - 1
@@ -646,8 +646,8 @@ def simulate_cnv_hmm(zarr_path, metadata_path, contigs, contig_sizes, rng):
     # - sample_is_high_variance [1D array] [bool] [True or False for n_samples]
     # - samples [1D array] [str]
 
-    # Get a random probability for a sample being high variance, between 0 and 1.
-    p_variance = rng.random()
+    # Keep high variance sample prevalence stable for deterministic tests.
+    p_variance = 0.1
 
     # Open a zarr at the specified path.
     root = zarr.open(zarr_path, mode="w")
@@ -862,8 +862,8 @@ def simulate_cnv_discordant_read_calls(
     # - sample_is_high_variance [1D array] [bool] [True or False for n_samples]
     # - samples [1D array] [str for n_samples]
 
-    # Get a random probability for a sample being high variance, between 0 and 1.
-    p_variance = rng.random()
+    # Keep high variance sample prevalence stable for deterministic tests.
+    p_variance = 0.1
 
     # Get a random probability for choosing allele 1, between 0 and 1.
     p_allele = rng.random()

tests/anoph/test_simulated_cnv_data.py

@@ -0,0 +1,81 @@
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+import zarr
+
+from .conftest import (
+    Gff3Simulator,
+    simulate_cnv_discordant_read_calls,
+    simulate_cnv_hmm,
+)
+
+
+def _write_sample_metadata(path: Path, n_samples: int = 100) -> None:
+    df_samples = pd.DataFrame({"sample_id": [f"S{i:04d}" for i in range(n_samples)]})
+    df_samples.to_csv(path, index=False)
+
+
+def test_simulate_cnv_hmm_limits_high_variance_fraction(tmp_path):
+    zarr_path = tmp_path / "cnv_hmm.zarr"
+    metadata_path = tmp_path / "samples.csv"
+    _write_sample_metadata(metadata_path)
+
+    simulate_cnv_hmm(
+        zarr_path=zarr_path,
+        metadata_path=metadata_path,
+        contigs=("2L",),
+        contig_sizes={"2L": 10_000},
+        rng=np.random.default_rng(0),
+    )
+
+    root = zarr.open(zarr_path, mode="r")
+    high_variance_fraction = np.mean(root["sample_is_high_variance"][:])
+    assert high_variance_fraction < 0.3
+
+
+def test_simulate_cnv_discordant_read_calls_limits_high_variance_fraction(tmp_path):
+    zarr_path = tmp_path / "cnv_discordant.zarr"
+    metadata_path = tmp_path / "samples.csv"
+    _write_sample_metadata(metadata_path)
+
+    simulate_cnv_discordant_read_calls(
+        zarr_path=zarr_path,
+        metadata_path=metadata_path,
+        contigs=("2L",),
+        contig_sizes={"2L": 10_000},
+        rng=np.random.default_rng(0),
+    )
+
+    root = zarr.open(zarr_path, mode="r")
+    high_variance_fraction = np.mean(root["sample_is_high_variance"][:])
+    assert high_variance_fraction < 0.3
+
+
+def test_simulate_exons_on_minus_strand_reverses_feature_order():
+    sim = Gff3Simulator(
+        contig_sizes={"2L": 10_000},
+        rng=np.random.default_rng(0),
+        n_exons_low=3,
+        n_exons_high=3,
+        intron_size_low=10,
+        intron_size_high=10,
+        exon_size_low=100,
+        exon_size_high=100,
+    )
+    rows = list(
+        sim.simulate_exons(
+            contig="2L",
+            strand="-",
+            gene_ix=0,
+            transcript_ix=0,
+            transcript_id="transcript-2L-0-0",
+            transcript_start=1,
+            transcript_end=1_000,
+        )
+    )
+    cds_and_utrs = [
+        row for row in rows if row[2] in {sim.utr5_type, sim.utr3_type, sim.cds_type}
+    ]
+    starts = [row[3] for row in cds_and_utrs]
+    assert starts == sorted(starts, reverse=True)