Make the sex-linked contig configurable in CNV frequency analysis

Fix #1313.

gene_cnv_frequencies() and gene_cnv_frequencies_advanced() hardcoded
`region.contig == "X"` when computing expected copy number. That is
correct for every currently supported assembly, but bakes in a
gambiae-complex assumption: any future species whose sex-linked contig
uses a different name would silently fall through to the diploid branch
and report incorrect amplification / deletion frequencies for males.

Changes:

- Add a `sex_contig` constructor parameter to
  `AnophelesGenomeSequenceData` (defaulting to "X") and expose it as a
  public `sex_contig` property, matching the pattern already used for
  `virtual_contigs`.
- Thread `sex_contig` through `AnophelesDataResource.__init__` so
  species subclasses can override it when their assembly differs.
- Replace both hardcoded `"X"` comparisons in cnv_frq.py with a shared
  `_expected_copy_number()` helper that consults `self._sex_contig`.
- Handle missing `sex_call` explicitly: if the sex-linked contig is
  requested but `sex_call` is absent from the sample metadata, emit a
  UserWarning and fall back to diploid for all samples rather than
  raising KeyError or silently returning a misleading scalar.

The default of "X" keeps behaviour unchanged for Ag3, Af1, As1, and
every other existing resource. New species (e.g. An. farauti) can opt
in by passing `sex_contig=` during construction.

Author: khushthecoder

malariagen_data/anoph/cnv_frq.py

@@ -29,6 +29,32 @@
 from .sample_metadata import _locate_cohorts
 
 
+def _expected_copy_number(contig, sex_contig, df_samples):
+    """Compute expected copy number per sample for a given contig.
+
+    On the sex-linked contig, males are hemizygous (CN=1) and females
+    are diploid (CN=2). On all other contigs, samples are diploid (CN=2).
+
+    If the sex-linked contig is requested but the ``sex_call`` column is
+    missing, a warning is emitted and all samples are treated as diploid
+    rather than silently returning a constant that would misrepresent
+    male hemizygous loci.
+    """
+    if contig != sex_contig:
+        return 2
+    if "sex_call" not in df_samples.columns:
+        warnings.warn(
+            "sex_call column not found; defaulting to diploid (CN=2) for "
+            "all samples on the sex-linked contig. Amplification and "
+            "deletion calls for hemizygous males may be incorrect.",
+            UserWarning,
+            stacklevel=3,
+        )
+        return 2
+    is_male = (df_samples["sex_call"] == "M").values
+    return np.where(is_male, 1, 2)[np.newaxis, :]
+
+
 class AnophelesCnvFrequencyAnalysis(AnophelesCnvData, AnophelesFrequencyAnalysis):
     def __init__(
         self,
@@ -290,11 +316,7 @@ def _gene_cnv_frequencies(
         df_samples = df_samples.set_index("sample_id").loc[sample_id].reset_index()
 
         debug("figure out expected copy number")
-        if region.contig == "X":
-            is_male = (df_samples["sex_call"] == "M").values
-            expected_cn = np.where(is_male, 1, 2)[np.newaxis, :]
-        else:
-            expected_cn = 2
+        expected_cn = _expected_copy_number(region.contig, self._sex_contig, df_samples)
 
         debug(
             "setup output dataframe - two rows for each gene, one for amplification and one for deletion"
@@ -557,11 +579,7 @@ def _gene_cnv_frequencies_advanced(
         )
 
         debug("figure out expected copy number")
-        if region.contig == "X":
-            is_male = (df_samples["sex_call"] == "M").values
-            expected_cn = np.where(is_male, 1, 2)[np.newaxis, :]
-        else:
-            expected_cn = 2
+        expected_cn = _expected_copy_number(region.contig, self._sex_contig, df_samples)
 
         debug("set up intermediates")
         cn = ds_cnv["CN_mode"].values

malariagen_data/anoph/genome_sequence.py

@@ -19,7 +19,10 @@ class AnophelesGenomeSequenceData(AnophelesBase):
     _cache_genome: Optional[zarr.hierarchy.Group]
 
     def __init__(
-        self, virtual_contigs: Optional[Mapping[str, Sequence[str]]] = None, **kwargs
+        self,
+        virtual_contigs: Optional[Mapping[str, Sequence[str]]] = None,
+        sex_contig: str = "X",
+        **kwargs,
     ):
         # N.B., this class is designed to work cooperatively, and
         # so it's important that any remaining parameters are passed
@@ -31,6 +34,12 @@ def __init__(
             virtual_contigs = dict()
         self._virtual_contigs = virtual_contigs
 
+        # Store the name of the sex-linked contig. Defaults to "X",
+        # which is correct for all currently supported assemblies, but
+        # can be overridden per species to support assemblies that use
+        # different contig naming conventions.
+        self._sex_contig = sex_contig
+
     @property
     def contigs(self) -> Tuple[str, ...]:
         """Contigs in the reference genome."""
@@ -41,6 +50,11 @@ def virtual_contigs(self) -> Mapping[str, Sequence[str]]:
         """Virtual contigs made by concatenating contigs in the reference genome."""
         return self._virtual_contigs
 
+    @property
+    def sex_contig(self) -> str:
+        """Name of the sex-linked contig used for ploidy-aware analyses."""
+        return self._sex_contig
+
     @property
     def _genome_zarr_path(self) -> str:
         return self.config["GENOME_ZARR_PATH"]

malariagen_data/anopheles.py

@@ -141,6 +141,7 @@ def __init__(
         taxon_colors: Optional[Mapping[str, str]] = None,
         aim_species_colors: Optional[Mapping[str, str]] = None,
         virtual_contigs: Optional[Mapping[str, Sequence[str]]] = None,
+        sex_contig: str = "X",
         gene_names: Optional[Mapping[str, str]] = None,
         inversion_tag_path: Optional[str] = None,
         unrestricted_use_only: Optional[bool] = None,
@@ -179,6 +180,7 @@ def __init__(
             taxon_colors=taxon_colors,
             aim_species_colors=aim_species_colors,
             virtual_contigs=virtual_contigs,
+            sex_contig=sex_contig,
             gene_names=gene_names,
             inversion_tag_path=inversion_tag_path,
             unrestricted_use_only=unrestricted_use_only,