feat: add --seed to simphenotype (#162)

CAST-genomics · Jan 10, 2023 · 9f7890a · 9f7890a
1 parent bcf778a
commit 9f7890a
Show file tree

Hide file tree

Showing 4 changed files with 39 additions and 5 deletions.
diff --git a/docs/api/data.rst b/docs/api/data.rst
@@ -580,11 +580,11 @@ To write to a **.bp** file, you must first initialize a :class:`Breakpoints` obj
 	breakpoints = data.Breakpoints('tests/data/example-write.bp')
 	breakpoints.data = {
 	    'HG00096': [
-	        np.array([('YRI','chr1',10114,4.3),('CEU','chr1',10116,5.2)], dtype=data.HapBlock)
-	        np.array([('CEU','chr1',10114,4.3),('YRI','chr1',10116,5.2)], dtype=data.HapBlock)
+	        np.array([('YRI','chr1',10114,4.3),('CEU','chr1',10116,5.2)], dtype=data.HapBlock),
+	        np.array([('CEU','chr1',10114,4.3),('YRI','chr1',10116,5.2)], dtype=data.HapBlock),
 	    ], 'HG00097': [
-	        np.array([('YRI','chr1',10114,4.3),('CEU','chr2',10116,5.2)], dtype=data.HapBlock)
-	        np.array([('CEU','chr1',10114,4.3),('YRI','chr2',10116,5.2)], dtype=data.HapBlock)
+	        np.array([('YRI','chr1',10114,4.3),('CEU','chr2',10116,5.2)], dtype=data.HapBlock),
+	        np.array([('CEU','chr1',10114,4.3),('YRI','chr2',10116,5.2)], dtype=data.HapBlock),
 	    ]
 	}
 	breakpoints.write()
diff --git a/haptools/__main__.py b/haptools/__main__.py
@@ -341,6 +341,13 @@ def simgenotype(
     show_default="all variants",
     help="If using a PGEN file, read genotypes in chunks of X variants; reduces memory",
 )
+@click.option(
+    "--seed",
+    type=int,
+    default=None,
+    show_default="chosen randomly",
+    help="Use this option across executions to make the output reproducible",
+)
 @click.option(
     "-o",
     "--output",
@@ -370,6 +377,7 @@ def simphenotype(
     ids: tuple[str] = tuple(),
     ids_file: Path = None,
     chunk_size: int = None,
+    seed: int = None,
     output: Path = Path("-"),
     verbosity: str = "INFO",
 ):
@@ -421,6 +429,7 @@ def simphenotype(
         samples,
         ids,
         chunk_size,
+        seed,
         output,
         log,
     )

diff --git a/haptools/sim_phenotype.py b/haptools/sim_phenotype.py
@@ -219,6 +219,7 @@ def simulate_pt(
     samples: list[str] = None,
     haplotype_ids: set[str] = None,
     chunk_size: int = None,
+    seed: int = None,
     output: Path = Path("-"),
     log: logging.Logger = None,
 ):
@@ -319,7 +320,7 @@ def simulate_pt(
 
     # Initialize phenotype simulator (haptools simphenotype)
     log.info("Simulating phenotypes")
-    pt_sim = PhenoSimulator(gt, output=output, log=log)
+    pt_sim = PhenoSimulator(gt, output=output, seed=seed, log=log)
     for i in range(num_replications):
         pt_sim.run(hp.data.values(), heritability, prevalence, normalize)
     log.info("Writing phenotypes")

diff --git a/tests/test_simphenotype.py b/tests/test_simphenotype.py
@@ -468,3 +468,27 @@ def test_no_normalize(self, capfd):
         assert result.exit_code == 0
 
         tmp_transform.unlink()
+
+    def test_seed(self, capfd):
+        # first, create a temporary file containing the output of transform
+        tmp_transform = Path("temp-transform.vcf")
+        with open(tmp_transform, "w") as file:
+            file.write(self._get_transform_stdin())
+
+        cmd = f"simphenotype --seed 42 {tmp_transform} tests/data/simple.hap"
+        runner = CliRunner()
+        result = runner.invoke(main, cmd.split(" "), catch_exceptions=False)
+        captured = capfd.readouterr()
+        assert captured.out
+        assert result.exit_code == 0
+
+        captured1 = captured.out
+
+        cmd = f"simphenotype --seed 42 {tmp_transform} tests/data/simple.hap"
+        runner = CliRunner()
+        result = runner.invoke(main, cmd.split(" "), catch_exceptions=False)
+        captured = capfd.readouterr()
+        assert captured.out == captured1
+        assert result.exit_code == 0
+
+        tmp_transform.unlink()