Merge pull request #911 from adilraza99/GH907-n-jack-validation

Guard against invalid jackknife block sizing when n_jack exceeds site count

Author: jonbrenas

malariagen_data/anoph/fst.py

@@ -404,8 +404,14 @@ def average_fst(
         )
 
         # Calculate block length for jackknife.
-        n_sites = ac1.shape[0]  # number of sites
-        block_length = n_sites // n_jack  # number of sites in each block
+        n_sites = ac1.shape[0]
+        block_length = n_sites // n_jack
+
+        if block_length < 1:
+            raise ValueError(
+                f"Not enough sites ({n_sites}) for {n_jack} jackknife blocks. "
+                "Choose a larger region or reduce n_jack."
+            )
 
         # Calculate average Fst.
         fst, se, _, _ = allel.blockwise_hudson_fst(ac1, ac2, blen=block_length)

tests/anoph/test_fst.py

@@ -190,6 +190,27 @@ def test_average_fst_with_min_cohort_size(fixture, api: AnophelesFstAnalysis):
         api.average_fst(**fst_params)
 
 
+@parametrize_with_cases("fixture,api", cases=".")
+def test_average_fst_region_too_small(fixture, api: AnophelesFstAnalysis):
+    """ValueError should be raised when block_length == 0 (n_jack > n_sites)."""
+    all_sample_sets = api.sample_sets()["sample_set"].to_list()
+    all_countries = api.sample_metadata()["country"].dropna().unique().tolist()
+    countries = random.sample(all_countries, 2)
+    cohort1_query = f"country == {countries[0]!r}"
+    cohort2_query = f"country == {countries[1]!r}"
+    fst_params = dict(
+        region=random.choice(api.contigs),
+        sample_sets=all_sample_sets,
+        cohort1_query=cohort1_query,
+        cohort2_query=cohort2_query,
+        site_mask=random.choice(api.site_mask_ids),
+        min_cohort_size=1,
+        n_jack=1_000_000,  # deliberately exceeds available sites
+    )
+    with pytest.raises(ValueError, match="Not enough sites"):
+        api.average_fst(**fst_params)
+
+
 def check_pairwise_average_fst(api: AnophelesFstAnalysis, fst_params):
     # Run main function under test.
     fst_df = api.pairwise_average_fst(**fst_params)