Merge branch 'master' into feature/775-heterozygosity-support

Author: Tanisha127

CONTRIBUTING.md

@@ -12,9 +12,10 @@ This package provides Python tools for accessing and analyzing genomic data from
 
 You'll need:
 
-- Python 3.10.x (CI-tested version)
-- [Poetry](https://python-poetry.org/) for dependency management
-- [Git](https://git-scm.com/) for version control
+- [pipx](https://python-poetry.org/) for installing Python tools
+- [git](https://git-scm.com/) for version control
+
+Both of these can be installed using your distribution's package manager or [Homebrew](https://brew.sh/) on Mac.
 
 ### Initial setup
 
@@ -33,16 +34,25 @@ You'll need:
    git remote add upstream https://github.com/malariagen/malariagen-data-python.git
    ```
 
-3. **Install Poetry** (if not already installed)
+3. **Install Poetry**
 
    ```bash
    pipx install poetry
    ```
 
-4. **Install the project and its dependencies**
+4. **Install Python 3.12**
+
+   Python 3.12 is tested in the CI-system and is the recommended version to use.
+
+   ```bash
+   poetry python install 3.12
+   ```
+
+5. **Install the project and its dependencies**
 
    ```bash
-   poetry install
+   poetry env use 3.12
+   poetry install --extras dev
    ```
 
    **Recommended**: Use `poetry run` to run commands inside the virtual environment:
@@ -71,7 +81,7 @@ You'll need:
    python script.py
    ```
 
-5. **Install pre-commit hooks**
+6. **Install pre-commit hooks**
 
    ```bash
    pipx install pre-commit
@@ -107,16 +117,40 @@ You'll need:
 
 4. **Run tests locally**
 
-   Fast unit tests (no external data access):
+   Fast unit tests using simulated data (no external data access):
 
    ```bash
-   poetry run pytest -v tests/anoph
+   poetry run pytest -v tests --ignore tests/integration
    ```
 
-   All unit tests (requires setting up credentials for legacy tests):
+   To run integration tests which read data from GCS, you'll need to [request access to MalariaGEN data on GCS](https://malariagen.github.io/vector-data/vobs/vobs-data-access.html).
+
+   Once access has been granted, [install the Google Cloud CLI](https://cloud.google.com/sdk/docs/install). E.g., if on Linux:
 
    ```bash
-   poetry run pytest -v tests --ignore tests/integration
+   ./install_gcloud.sh
+   ```
+
+   You'll then need to obtain application-default credentials, e.g.:
+
+   ```bash
+   ./google-cloud-sdk/bin/gcloud auth application-default login
+   ```
+
+   Once this is done, you can run integration tests:
+
+   ```bash
+   poetry run pytest -v tests/integration
+   ```
+
+   Tests will run slowly the first time, as data required for testing will be read from GCS. Subsequent runs will be faster as data will be cached locally in the "gcs_cache" folder.
+
+6. **Run typechecking**
+
+   Run static typechecking with mypy:
+
+   ```bash
+   poetry run mypy malariagen_data tests --ignore-missing-imports
    ```
 
 5. **Check code quality**
@@ -150,7 +184,7 @@ ruff format .
 - **Fast tests**: Unit tests should use simulated data when possible (see `tests/anoph/`)
 - **Integration tests**: Tests requiring GCS data access are slower and run separately
 
-Run type checking with:
+Run dynamic type checking with:
 
 ```bash
 poetry run pytest -v tests --typeguard-packages=malariagen_data,malariagen_data.anoph

malariagen_data/anoph/distance.py

@@ -365,24 +365,43 @@ def _njt(
         from scipy.spatial.distance import squareform  # type: ignore
 
         # Compute pairwise distances.
-        dist, samples, n_snps = self.biallelic_diplotype_pairwise_distances(
-            region=region,
-            n_snps=n_snps,
-            metric=metric,
-            sample_sets=sample_sets,
-            sample_indices=sample_indices,
-            site_mask=site_mask,
-            site_class=site_class,
-            inline_array=inline_array,
-            chunks=chunks,
-            cohort_size=cohort_size,
-            min_cohort_size=min_cohort_size,
-            max_cohort_size=max_cohort_size,
-            random_seed=random_seed,
-            max_missing_an=max_missing_an,
-            min_minor_ac=min_minor_ac,
-            thin_offset=thin_offset,
-        )
+        try:
+            dist, samples, n_snps_used = self.biallelic_diplotype_pairwise_distances(
+                region=region,
+                n_snps=n_snps,
+                metric=metric,
+                sample_sets=sample_sets,
+                sample_indices=sample_indices,
+                site_mask=site_mask,
+                site_class=site_class,
+                inline_array=inline_array,
+                chunks=chunks,
+                cohort_size=cohort_size,
+                min_cohort_size=min_cohort_size,
+                max_cohort_size=max_cohort_size,
+                random_seed=random_seed,
+                max_missing_an=max_missing_an,
+                min_minor_ac=min_minor_ac,
+                thin_offset=thin_offset,
+            )
+
+        except ValueError as e:
+            raise ValueError(
+                f"Unable to construct neighbour-joining tree. {e} "
+                f"This could be because the selected region does not "
+                f"contain enough polymorphic SNPs for the given sample "
+                f"sets and query parameters."
+            ) from e
+
+        # Validate enough samples for a tree.
+        n_samples = len(samples)
+        if n_samples < 3:
+            raise ValueError(
+                f"Not enough samples to construct a neighbour-joining tree. "
+                f"A minimum of 3 samples is required, but only {n_samples} "
+                f"were found for the given region and sample sets."
+            )
+
         D = squareform(dist)
 
         # anjl supports passing in a progress bar function to get progress on the

malariagen_data/anoph/frq_base.py

@@ -622,13 +622,27 @@ def plot_frequencies_interactive_map(
         variants = ds["variant_label"].values
         taxa = ds["cohort_taxon"].to_pandas().dropna().unique()  # type: ignore
         periods = ds["cohort_period"].to_pandas().dropna().unique()  # type: ignore
+
+        if len(variants) == 0:
+            raise ValueError("No variants available in dataset.")
+        if len(taxa) == 0:
+            raise ValueError("No taxons available in dataset.")
+        if len(periods) == 0:
+            raise ValueError("No periods available in dataset.")
+
         controls = ipywidgets.interactive(
             self.plot_frequencies_map_markers,
             m=ipywidgets.fixed(freq_map),
             ds=ipywidgets.fixed(ds),
-            variant=ipywidgets.Dropdown(options=variants, description="Variant: "),
-            taxon=ipywidgets.Dropdown(options=taxa, description="Taxon: "),
-            period=ipywidgets.Dropdown(options=periods, description="Period: "),
+            variant=ipywidgets.Dropdown(
+                options=variants, value=variants[0], description="Variant: "
+            ),
+            taxon=ipywidgets.Dropdown(
+                options=taxa, value=taxa[0], description="Taxon: "
+            ),
+            period=ipywidgets.Dropdown(
+                options=periods, value=periods[0], description="Period: "
+            ),
             clear=ipywidgets.fixed(True),
         )
 

malariagen_data/anoph/snp_data.py

@@ -1903,7 +1903,12 @@ def biallelic_snp_calls(
                     ds_out = ds_out.isel(variants=loc_thin)
 
                 elif ds_out.sizes["variants"] < n_snps:
-                    raise ValueError("Not enough SNPs.")
+                    raise ValueError(
+                        f"Not enough SNPs. Requested {n_snps} SNPs but only "
+                        f"{ds_out.sizes['variants']} were found in the selected "
+                        f"region after applying filters. Try using a larger region, "
+                        f"relaxing site filters, or reducing the n_snps parameter."
+                    )
 
         return ds_out
 

tests/anoph/test_cnv_frq.py

@@ -643,7 +643,7 @@ def test_gene_cnv_frequencies_advanced_with_variant_query(
     area_by = "admin1_iso"
     period_by = "year"
     region = random.choice(api.contigs)
-    variant_query = "cnv_type == '{variant_query_option}'"
+    variant_query = f"cnv_type == '{variant_query_option}'"
 
     check_gene_cnv_frequencies_advanced(
         api=api,
@@ -743,7 +743,8 @@ def check_gene_cnv_frequencies_advanced(
     check_plot_frequencies_time_series(api, ds)
     check_plot_frequencies_time_series_with_taxa(api, ds)
     check_plot_frequencies_time_series_with_areas(api, ds)
-    check_plot_frequencies_interactive_map(api, ds)
+    if variant_query is None:
+        check_plot_frequencies_interactive_map(api, ds)
     assert set(ds.dims) == {"cohorts", "variants"}
 
     # Check variant variables.

tests/anoph/test_distance.py

@@ -269,3 +269,20 @@ def test_plot_njt(fixture, api: AnophelesDistanceAnalysis):
             **data_params,
         )
         assert isinstance(fig, go.Figure)
+
+
+@parametrize_with_cases("fixture,api", cases=".")
+def test_njt_not_enough_snps(fixture, api: AnophelesDistanceAnalysis):
+    all_sample_sets = api.sample_sets()["sample_set"].to_list()
+    with pytest.raises(
+        ValueError,
+        match="Unable to construct neighbour-joining tree|Not enough SNPs",
+    ):
+        api.njt(
+            region=random.choice(api.contigs),
+            n_snps=1_000_000_000,  # impossibly high to guarantee failure
+            sample_sets=random.sample(all_sample_sets, 1),
+            site_mask=random.choice((None,) + api.site_mask_ids),
+            min_minor_ac=pca_params.min_minor_ac_default,
+            max_missing_an=pca_params.max_missing_an_default,
+        )