Merge branch 'master' into GH985-fix-stray-f-in-error-messages

Author: jonbrenas

.codecov.yml

@@ -0,0 +1,9 @@
+coverage:
+  status:
+    project:
+      default:
+        target: auto
+    patch:
+      default:
+        target: 80%
+        threshold: 0%

malariagen_data/adir1.py

@@ -19,6 +19,10 @@
     "dirus": TAXON_PALETTE[0],
 }
 
+XPEHH_GWSS_CACHE_NAME = "adir1_xpehh_gwss_v1"
+IHS_GWSS_CACHE_NAME = "adir1_ihs_gwss_v1"
+ROH_HMM_CACHE_NAME = "adir1_roh_hmm_v1"
+
 
 class Adir1(AnophelesDataResource):
     """Provides access to data from Adir1.0 releases.
@@ -71,6 +75,10 @@ class Adir1(AnophelesDataResource):
 
     """
 
+    _xpehh_gwss_cache_name = XPEHH_GWSS_CACHE_NAME
+    _ihs_gwss_cache_name = IHS_GWSS_CACHE_NAME
+    _roh_hmm_cache_name = ROH_HMM_CACHE_NAME
+
     def __init__(
         self,
         url=None,

malariagen_data/af1.py

@@ -21,6 +21,10 @@
     "funestus": TAXON_PALETTE[0],
 }
 
+XPEHH_GWSS_CACHE_NAME = "af1_xpehh_gwss_v1"
+IHS_GWSS_CACHE_NAME = "af1_ihs_gwss_v1"
+ROH_HMM_CACHE_NAME = "af1_roh_hmm_v1"
+
 
 class Af1(AnophelesDataResource):
     """Provides access to data from Af1.x releases.
@@ -75,6 +79,7 @@ class Af1(AnophelesDataResource):
 
     _xpehh_gwss_cache_name = XPEHH_GWSS_CACHE_NAME
     _ihs_gwss_cache_name = IHS_GWSS_CACHE_NAME
+    _roh_hmm_cache_name = ROH_HMM_CACHE_NAME
 
     def __init__(
         self,

malariagen_data/ag3.py

@@ -95,6 +95,10 @@ def _setup_aim_palettes():
     "aim_species": "object",
 }
 
+XPEHH_GWSS_CACHE_NAME = "ag3_xpehh_gwss_v1"
+IHS_GWSS_CACHE_NAME = "ag3_ihs_gwss_v1"
+ROH_HMM_CACHE_NAME = "ag3_roh_hmm_v1"
+
 
 class Ag3(AnophelesDataResource):
     """Provides access to data from Ag3.x releases.
@@ -153,6 +157,7 @@ class Ag3(AnophelesDataResource):
 
     _xpehh_gwss_cache_name = XPEHH_GWSS_CACHE_NAME
     _ihs_gwss_cache_name = IHS_GWSS_CACHE_NAME
+    _roh_hmm_cache_name = ROH_HMM_CACHE_NAME
 
     def __init__(
         self,

malariagen_data/amin1.py

@@ -19,6 +19,10 @@
     "dirus": TAXON_PALETTE[0],
 }
 
+XPEHH_GWSS_CACHE_NAME = "amin1_xpehh_gwss_v1"
+IHS_GWSS_CACHE_NAME = "amin1_ihs_gwss_v1"
+ROH_HMM_CACHE_NAME = "amin1_roh_hmm_v1"
+
 
 class Amin1(AnophelesDataResource):
     """Provides access to data from Amin1.0 releases.
@@ -71,8 +75,9 @@ class Amin1(AnophelesDataResource):
 
     """
 
-    #    _xpehh_gwss_cache_name = XPEHH_GWSS_CACHE_NAME
-    #    _ihs_gwss_cache_name = IHS_GWSS_CACHE_NAME
+    _xpehh_gwss_cache_name = XPEHH_GWSS_CACHE_NAME
+    _ihs_gwss_cache_name = IHS_GWSS_CACHE_NAME
+    _roh_hmm_cache_name = ROH_HMM_CACHE_NAME
 
     def __init__(
         self,

malariagen_data/anoph/base.py

@@ -562,6 +562,13 @@ def _sample_set_has_unrestricted_use(self, *, sample_set: str):
         release_manifest_df = self._read_sample_sets_manifest(
             single_release=sample_set_release
         )
+
+        if "unrestricted_use" not in release_manifest_df.columns:
+            raise ValueError(
+                f"Column 'unrestricted_use' missing from manifest for sample set '{sample_set}'. "
+                "This indicates a data integrity issue in the release manifest."
+            )
+
         sample_set_records_srs = release_manifest_df.loc[
             release_manifest_df["sample_set"] == sample_set, "unrestricted_use"
         ]
@@ -824,12 +831,19 @@ def lookup_study_info(self, sample_set: base_params.sample_set) -> dict:
     def lookup_terms_of_use_info(self, sample_set: base_params.sample_set) -> dict:
         if self._cache_sample_set_to_terms_of_use_info is None:
             df_sample_sets = self._available_sample_sets().set_index("sample_set")
+            expected_cols = [
+                "terms_of_use_expiry_date",
+                "terms_of_use_url",
+                "unrestricted_use",
+            ]
+            missing_cols = [c for c in expected_cols if c not in df_sample_sets.columns]
+            if missing_cols:
+                raise ValueError(
+                    f"Terms-of-use columns missing from manifest: {missing_cols}. "
+                    "This indicates a data integrity issue in the release manifest."
+                )
             self._cache_sample_set_to_terms_of_use_info = df_sample_sets[
-                [
-                    "terms_of_use_expiry_date",
-                    "terms_of_use_url",
-                    "unrestricted_use",
-                ]
+                expected_cols
             ].to_dict(orient="index")
         try:
             return self._cache_sample_set_to_terms_of_use_info[sample_set]

malariagen_data/anoph/distance.py

@@ -217,6 +217,9 @@ def _biallelic_diplotype_pairwise_distances(
         n_snps = gn.shape[0]
 
         # Prepare data for pairwise distance calculation.
+        # Mask missing calls (-127) before computing distances.
+        gn = gn.astype(float)
+        gn[gn == -127] = np.nan
         X = np.ascontiguousarray(gn.T)
 
         # Look up distance function.

malariagen_data/anoph/g123.py

@@ -133,6 +133,13 @@ def _g123_gwss(
             chunks=chunks,
         )
 
+        if gt.shape[0] < window_size:
+            raise ValueError(
+                f"Not enough sites ({gt.shape[0]}) for window size "
+                f"({window_size}). Please reduce the window size or "
+                f"use different site selection criteria."
+            )
+
         with self._spinner("Compute G123"):
             g123 = allel.moving_statistic(gt, statistic=_garud_g123, size=window_size)
             x = allel.moving_statistic(pos, statistic=np.mean, size=window_size)
@@ -240,6 +247,12 @@ def _g123_calibration(
 
         calibration_runs: Dict[str, np.ndarray] = dict()
         for window_size in self._progress(window_sizes, desc="Compute G123"):
+            if gt.shape[0] < window_size:
+                raise ValueError(
+                    f"Not enough sites ({gt.shape[0]}) for window size "
+                    f"({window_size}). Please reduce the window size or "
+                    f"use different site selection criteria."
+                )
             g123 = allel.moving_statistic(gt, statistic=_garud_g123, size=window_size)
             calibration_runs[str(window_size)] = g123
 

malariagen_data/anoph/h12.py

@@ -558,6 +558,8 @@ def plot_h12_gwss_multi_overlay_track(
         show: gplt_params.show = True,
         x_range: Optional[gplt_params.x_range] = None,
         output_backend: gplt_params.output_backend = gplt_params.output_backend_default,
+        chunks: base_params.chunks = base_params.native_chunks,
+        inline_array: base_params.inline_array = base_params.inline_array_default,
     ) -> gplt_params.optional_figure:
         cohort_queries = self._setup_cohort_queries(
             cohorts=cohorts,
@@ -585,8 +587,11 @@ def plot_h12_gwss_multi_overlay_track(
                 min_cohort_size=min_cohort_size,
                 max_cohort_size=max_cohort_size,
                 sample_query=cohort_query,
+                sample_query_options=sample_query_options,
                 sample_sets=sample_sets,
                 random_seed=random_seed,
+                chunks=chunks,
+                inline_array=inline_array,
             )
 
         # Determine X axis range.
@@ -625,7 +630,7 @@ def plot_h12_gwss_multi_overlay_track(
         )
 
         # Plot H12.
-        for i, (cohort_label, (x, h12, contig)) in enumerate(res.items()):
+        for i, (cohort_label, (x, h12, contig_idx)) in enumerate(res.items()):
             fig.scatter(
                 x=x,
                 y=h12,
@@ -679,6 +684,8 @@ def plot_h12_gwss_multi_overlay(
         output_backend: gplt_params.output_backend = gplt_params.output_backend_default,
         gene_labels: Optional[gplt_params.gene_labels] = None,
         gene_labelset: Optional[gplt_params.gene_labelset] = None,
+        chunks: base_params.chunks = base_params.native_chunks,
+        inline_array: base_params.inline_array = base_params.inline_array_default,
     ) -> gplt_params.optional_figure:
         # Plot GWSS track.
         fig1 = self.plot_h12_gwss_multi_overlay_track(
@@ -700,6 +707,8 @@ def plot_h12_gwss_multi_overlay(
             height=track_height,
             show=False,
             output_backend=output_backend,
+            chunks=chunks,
+            inline_array=inline_array,
         )
 
         fig1.xaxis.visible = False
@@ -764,6 +773,8 @@ def plot_h12_gwss_multi_panel(
         output_backend: gplt_params.output_backend = gplt_params.output_backend_default,
         gene_labels: Optional[gplt_params.gene_labels] = None,
         gene_labelset: Optional[gplt_params.gene_labelset] = None,
+        chunks: base_params.chunks = base_params.native_chunks,
+        inline_array: base_params.inline_array = base_params.inline_array_default,
     ) -> gplt_params.optional_figure:
         cohort_queries = self._setup_cohort_queries(
             cohorts=cohorts,
@@ -789,6 +800,7 @@ def plot_h12_gwss_multi_panel(
                 window_size=window_size[cohort_label],
                 sample_sets=sample_sets,
                 sample_query=cohort_query,
+                sample_query_options=sample_query_options,
                 cohort_size=cohort_size,
                 min_cohort_size=min_cohort_size,
                 max_cohort_size=max_cohort_size,
@@ -799,6 +811,8 @@ def plot_h12_gwss_multi_panel(
                 height=track_height,
                 show=False,
                 output_backend=output_backend,
+                chunks=chunks,
+                inline_array=inline_array,
             )
             if i > 0:
                 track = self.plot_h12_gwss_track(x_range=figs[0].x_range, **params)

malariagen_data/anoph/pca.py

@@ -44,6 +44,14 @@ def __init__(
             `random_seed`.
 
         """,
+        parameters=dict(
+            imputation_method="""
+                Method to use for imputing missing genotype calls. Options are
+                'most_common' (replace missing calls with the most common genotype at each site,
+                the default), 'mean' (replace missing calls with the
+                mean value at each site), or 'zero' (replace missing calls with zero).
+            """,
+        ),
         returns=("df_pca", "evr"),
         notes="""
             This computation may take some time to run, depending on your computing
@@ -69,6 +77,7 @@ def pca(
         max_missing_an: Optional[
             base_params.max_missing_an
         ] = pca_params.max_missing_an_default,
+        imputation_method: pca_params.imputation_method = pca_params.imputation_method_default,
         cohort_size: Optional[base_params.cohort_size] = None,
         min_cohort_size: Optional[base_params.min_cohort_size] = None,
         max_cohort_size: Optional[base_params.max_cohort_size] = None,
@@ -80,7 +89,7 @@ def pca(
     ) -> Tuple[pca_params.df_pca, pca_params.evr]:
         # Change this name if you ever change the behaviour of this function, to
         # invalidate any previously cached data.
-        name = "pca_v5"
+        name = "pca_v8"
 
         # Check that either sample_query xor sample_indices are provided.
         base_params._validate_sample_selection_params(
@@ -121,6 +130,7 @@ def pca(
             site_class=site_class,
             min_minor_ac=min_minor_ac,
             max_missing_an=max_missing_an,
+            imputation_method=imputation_method,
             n_components=n_components,
             cohort_size=cohort_size,
             min_cohort_size=min_cohort_size,
@@ -152,7 +162,7 @@ def pca(
         # df_pca.index = df_pca.index.astype(str)
 
         # Name the DataFrame's columns as PC1, PC2, etc.
-        df_pca.columns = pd.Index([f"PC{i+1}" for i in range(coords.shape[1])])
+        df_pca.columns = pd.Index([f"PC{i + 1}" for i in range(coords.shape[1])])
 
         # Load the sample metadata.
         df_samples = self.sample_metadata(
@@ -185,6 +195,7 @@ def _pca(
         site_class,
         min_minor_ac,
         max_missing_an,
+        imputation_method="most_common",
         n_components,
         cohort_size,
         min_cohort_size,
@@ -231,6 +242,50 @@ def _pca(
                 loc_keep_fit = np.ones(len(samples), dtype=bool)
                 gn_fit = gn
 
+            # Impute missing calls (-127) using the chosen imputation method.
+            if max_missing_an is not None and max_missing_an > 0:
+                gn_fit = gn_fit.astype(float)
+                gn = gn.astype(float)
+                for arr in [gn_fit, gn]:
+                    missing_mask = arr == -127
+
+                    if imputation_method == "most_common":
+                        # For each site, find the most common non-missing value.
+                        site_modes = []
+                        for row in arr:
+                            non_missing = row[row != -127]
+                            if len(non_missing) == 0:
+                                site_modes.append(0)
+                            else:
+                                values, counts = np.unique(
+                                    non_missing, return_counts=True
+                                )
+                                site_modes.append(values[np.argmax(counts)])
+                        site_modes = np.array(site_modes)
+                        fill_values = np.take(site_modes, np.where(missing_mask)[0])
+                    elif imputation_method == "mean":
+                        site_means = np.where(
+                            np.all(missing_mask, axis=1, keepdims=True),
+                            0,
+                            np.nanmean(
+                                np.where(missing_mask, np.nan, arr),
+                                axis=1,
+                                keepdims=True,
+                            ),
+                        )
+                        fill_values = np.take(
+                            site_means.flatten(), np.where(missing_mask)[0]
+                        )
+                    elif imputation_method == "zero":
+                        fill_values = 0
+                    else:
+                        raise ValueError(
+                            f"Unknown imputation_method: {imputation_method!r}. "
+                            "Choose from 'most_common', 'mean' or 'zero'."
+                        )
+
+                    arr[missing_mask] = fill_values
+
             # Remove any sites where all genotypes are identical.
             loc_var = np.any(gn_fit != gn_fit[:, 0, np.newaxis], axis=1)
             gn_fit_var = np.compress(loc_var, gn_fit, axis=0)