Highlight each locus of a gene name in scatter --gene

cnvlib/cnary.py

@@ -36,7 +36,7 @@ def is_female_default(is_xx: bool | bool_ | None) -> bool:
     return True if is_xx is None else bool(is_xx)
 
 
-def _gene_runs(gene_col: Series, ignore: tuple[str, ...]) -> list[tuple[int, int, str]]:
+def gene_runs(gene_col: Series, ignore: tuple[str, ...]) -> list[tuple[int, int, str]]:
     """Locate each gene's maximal runs of bins within one chromosome.
 
     A run for gene *g* opens at the first bin naming *g*, survives bins that
@@ -258,7 +258,7 @@ def by_gene(
         ignore = (*ignore, *params.ANTITARGET_ALIASES)
         for _chrom, subgary in self.by_chromosome():
             prev_pos = 0
-            for start_pos, end_pos, gene in _gene_runs(subgary["gene"], ignore):
+            for start_pos, end_pos, gene in gene_runs(subgary["gene"], ignore):
                 if prev_pos < start_pos:
                     # Include intergenic regions
                     yield (

cnvlib/core.py

@@ -105,16 +105,6 @@ def assert_equal(msg: str, **values) -> None:
         raise ValueError(msg)
 
 
-def check_unique(items, title):
-    """Ensure all items in an iterable are identical; return that one item."""
-    its = set(items)
-    assert len(its) == 1, "Inconsistent %s keys: %s" % (
-        title,
-        " ".join(map(str, sorted(its))),
-    )
-    return its.pop()
-
-
 def fbase(fname: str) -> str:
     """Strip directory and all extensions from a filename."""
     base = os.path.basename(fname)

cnvlib/plots.py

@@ -12,7 +12,7 @@
 from skgenome.intersect import require_ascending_starts
 from skgenome.rangelabel import Region, unpack_range
 
-from . import core, params
+from . import cnary, params
 
 if TYPE_CHECKING:
     from matplotlib.axes._axes import Axes
@@ -377,45 +377,43 @@ def cvg2rgb(cvg: float, desaturate: bool) -> tuple[float, float, float]:
     return rgb
 
 
-# XXX should this be a CopyNumArray method?
-# or: use by_genes internally
-# or: have by_genes use this internally
 def gene_coords_by_name(probes, names):
     """Find the chromosomal position of each named gene in probes.
 
+    Each locus of a requested name gets its own region, grouped exactly as
+    `CopyNumArray.by_gene` groups bins, via `cnvlib.cnary.gene_runs`. A name
+    that recurs at several loci -- a repeat family, or a gene name reused on
+    another chromosome -- therefore yields one region per locus instead of one
+    span reaching from its first occurrence to its last.
+
     Returns
     -------
     dict
         Of: {chromosome: [(start, end, gene name), ...]}
     """
-    names = list(filter(None, set(names)))
+    names = set(filter(None, names))
     if not names:
         return {}
 
-    # Create an index of gene names
-    gene_index = collections.defaultdict(set)
-    for i, gene in enumerate(probes["gene"]):
-        for gene_name in gene.split(","):
-            if gene_name in names:
-                gene_index[gene_name].add(i)
-    # Retrieve coordinates by name
+    ignore = (*params.IGNORE_GENE_NAMES, *params.ANTITARGET_ALIASES)
+    # A locus spans its run's first and last bin -- the same start and end
+    # 'genemetrics' publishes for that locus (reports.group_by_genes).
+    # Regions are keyed so that requested names sharing a locus, e.g. the
+    # co-binned pair "ERBB2,MIR4728", are labeled once; co-binned names the
+    # caller did not ask for are never surfaced (#458).
     all_coords: dict = collections.defaultdict(lambda: collections.defaultdict(set))
-    for name in names:
-        gene_probes = probes.data.take(sorted(gene_index.get(name, [])))
-        if not len(gene_probes):
-            raise ValueError(f"No targeted gene named {name!r} found")
-        # Find the genomic range of this gene's probes
-        start = gene_probes["start"].min()
-        end = gene_probes["end"].max()
-        chrom = core.check_unique(gene_probes["chromosome"], name)
-        # Deduce the unique set of *requested* gene names for this region.
-        # A bin label may pack several genes (e.g. "ERBB2,MIR4728"); only the
-        # names the caller asked for should be surfaced, never co-binned
-        # neighbors (#458).
-        uniq_names: set[str] = set()
-        for oname in gene_probes["gene"].dropna().unique():
-            uniq_names.update(g for g in oname.split(",") if g in names)
-        all_coords[chrom][start, end].update(uniq_names)
+    found: set[str] = set()
+    for chrom, subarr in probes.by_chromosome():
+        starts = subarr["start"].to_numpy()
+        ends = subarr["end"].to_numpy()
+        for start_pos, end_pos, gene in cnary.gene_runs(subarr["gene"], ignore):
+            if gene in names:
+                found.add(gene)
+                all_coords[chrom][starts[start_pos], ends[end_pos - 1]].add(gene)
+    if unfound := names - found:
+        raise ValueError(
+            "No targeted gene named " + ", ".join(map(repr, sorted(unfound))) + " found"
+        )
     # Consolidate each region's gene names into a string
     uniq_coords = {}
     for chrom, hits in all_coords.items():

cnvlib/scatter.py

@@ -2,16 +2,18 @@
 
 from __future__ import annotations
 
+import collections
 import logging
 from typing import TYPE_CHECKING
 
 import numpy as np
 from matplotlib import pyplot
 
 from skgenome.chromnames import detect_chr_prefix, diagnose_missing_chromosome
+from skgenome.chromsort import sorter_chrom
 from skgenome.rangelabel import unpack_range
 
-from . import core, params, plots
+from . import params, plots
 from .cnary import CopyNumArray as CNA
 from .plots import MB
 
@@ -607,6 +609,15 @@ def chromosome_scatter(
     return axis.get_figure()
 
 
+def _region_gene_names(regions) -> set[str]:
+    """The individual gene names labeling `regions`.
+
+    A region's label packs the requested names sharing that locus, e.g. the
+    co-binned pair "ERBB2,MIR4728".
+    """
+    return {name for _start, _end, label in regions for name in label.split(",")}
+
+
 def select_range_genes(
     cnarr,
     segments,
@@ -623,11 +634,18 @@ def select_range_genes(
     Behaviors::
 
         start/end   show_gene
-           +           +       given region + genes; err if any gene outside it
-           -           +       window +/- around genes
+           +           +       given region + the genes' loci within it; err if
+                               a named gene has no locus there
+           -           +       window +/- around the named genes' loci
            +           -       given region, highlighting any genes within it
            -           -       whole chromosome, no genes
 
+    A gene name is highlighted at each locus where it occurs, one region per
+    locus, so a repeat-family name yields many highlights rather than one span
+    reaching across all of them. A name occurring on several chromosomes needs
+    a chromosome selected by `show_range`, since one panel plots one
+    chromosome.
+
     If `show_range` is a chromosome name only, no start/end positions, then the
     whole chromosome will be shown.
 
@@ -677,41 +695,78 @@ def select_range_genes(
                 window_coords = (max(0, start - window_width), end + window_width)
 
     else:
-        gene_names = filter(None, show_gene.split(","))
+        gene_names = list(filter(None, show_gene.split(",")))
         if gene_names:
             # Scan for probes matching the specified gene(s)
             gene_coords = plots.gene_coords_by_name(cnarr or segments, gene_names)
-            if len(gene_coords) > 1:
+            if chrom:
+                # A name can occur on more than one chromosome; -c picks among
+                # them, so only its loci there are plotted.
+                gene_ranges = gene_coords.get(chrom, [])
+                absent = sorted(set(gene_names) - _region_gene_names(gene_ranges))
+                if absent:
+                    elsewhere = sorted(
+                        (
+                            other
+                            for other, regions in gene_coords.items()
+                            if _region_gene_names(regions).intersection(absent)
+                        ),
+                        key=sorter_chrom,
+                    )
+                    raise ValueError(
+                        f"Gene(s) {', '.join(absent)} not found on chromosome "
+                        f"{chrom} selected by region; found on "
+                        f"{', '.join(elsewhere)}"
+                    )
+            elif len(gene_coords) > 1:
                 raise ValueError(
                     f"Genes {show_gene} are split across chromosomes "
-                    f"{list(gene_coords.keys())}"
-                )
-            g_chrom, gene_ranges = gene_coords.popitem()
-            if chrom:
-                # Confirm that the selected chromosomes match
-                core.assert_equal(
-                    "Chromosome also selected by region (-c) does not match",
-                    **{"chromosome": chrom, "gene(s)": g_chrom},
+                    f"{sorted(gene_coords, key=sorter_chrom)}; "
+                    "select one with -c/--chromosome"
                 )
             else:
-                chrom = g_chrom
+                chrom, gene_ranges = gene_coords.popitem()
 
             gene_ranges.sort()
             if window_coords:
-                # Verify all genes fit in the given window
-                for gene_start, gene_end, gene_name in gene_ranges:
-                    if not (start <= gene_start and gene_end <= end):
-                        raise ValueError(
-                            f"Selected gene {gene_name} "
-                            + f"({chrom}:{gene_start}-{gene_end}) "
-                            + f"is outside specified region {show_range}"
-                        )
+                # Highlight the loci inside the given region; a gene with none
+                # there is the mistake this check exists to catch.
+                inside = [
+                    region
+                    for region in gene_ranges
+                    if start <= region[0] and region[1] <= end
+                ]
+                inside_names = _region_gene_names(inside)
+                for gene_start, gene_end, label in gene_ranges:
+                    for gene_name in label.split(","):
+                        if gene_name not in inside_names:
+                            raise ValueError(
+                                f"Selected gene {gene_name} "
+                                + f"({chrom}:{gene_start}-{gene_end}) "
+                                + f"is outside specified region {show_range}"
+                            )
+                if dropped := len(gene_ranges) - len(inside):
+                    logging.info(
+                        "Not highlighting %d gene loci outside the selected region",
+                        dropped,
+                    )
+                gene_ranges = inside
             elif not show_range:
                 # Set the display window to the selected genes +/- a margin
                 window_coords = (
                     max(0, gene_ranges[0][0] - window_width),
                     gene_ranges[-1][1] + window_width,
                 )
+            for gene_name, n_loci in collections.Counter(
+                name for _s, _e, label in gene_ranges for name in label.split(",")
+            ).items():
+                if n_loci > 1:
+                    logging.warning(
+                        "Gene %s occurs at %d separate loci on %s; highlighting each",
+                        gene_name,
+                        n_loci,
+                        chrom,
+                    )
 
     # Prune plotted elements to the selected region
     sel_probes = cnarr.in_range(chrom, *window_coords) if cnarr else CNA([])
@@ -724,7 +779,7 @@ def select_range_genes(
         somatic_variants.in_range(chrom, *window_coords) if somatic_variants else None
     )
     logging.info(
-        "Showing %d probes and %d selected genes in region %s",
+        "Showing %d probes and %d highlighted gene regions in region %s",
         len(sel_probes),
         len(gene_ranges),
         (chrom + ":{}-{}".format(*window_coords) if window_coords else chrom),

doc/plots.rst

@@ -93,7 +93,15 @@ Chromosome-level views are controlled with the ``--chromosome``/``-c`` and
   (e.g. ``-g CDK4,MDM2``) will plot the genomic around that gene, or genes, and
   highlight the gene or genes with a vertical gold stripe.
 
-    - If multiple genes, they must all be on the same chromosome.
+    - Each locus where a name occurs is highlighted separately, grouped as the
+      ``genemetrics`` command groups bins into genes. A name shared by many
+      loci -- a repeat family such as ``Y_RNA``, say -- therefore yields one
+      stripe per locus rather than one stripe reaching from the first
+      occurrence to the last.
+    - If multiple genes, they must all be on the same chromosome. A name
+      occurring on more than one chromosome likewise needs one chosen with
+      ``-c``, e.g. ``-c chr7 -g <name>``; only its loci on that chromosome are
+      plotted.
     - The ``--width``/``-w`` argument determines the size of the plotted
       genomic region, in terms of basepairs flanking the selected region.
     - Any other genes in the plotted region will not be shown unless also
@@ -115,7 +123,9 @@ Chromosome-level views are controlled with the ``--chromosome``/``-c`` and
       will be shown, with no padding.
     - The ``-g`` option overrides the default behavior of showing all genes in
       the selection -- only the genes specified with ``-g`` will be highlighted
-      and labeled. To not show any genes, specify an empty string: ``-g ''``
+      and labeled. To not show any genes, specify an empty string: ``-g ''``.
+      Only the loci within the selected region are highlighted; if a named gene
+      has no locus there at all, that is reported as an error.
     - Special behavior occurs if there are no genes in the selected region:
       Instead, the selection itself is treated as a "gene", highlighted and
       labeled with the string "Selection", with padding controlled by ``-w``.