diff --git a/cnvlib/cnary.py b/cnvlib/cnary.py index fc11000e..9d7d20be 100644 --- a/cnvlib/cnary.py +++ b/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 ( diff --git a/cnvlib/core.py b/cnvlib/core.py index f8e1d610..2986a501 100644 --- a/cnvlib/core.py +++ b/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) diff --git a/cnvlib/plots.py b/cnvlib/plots.py index 3b31a206..8279b375 100644 --- a/cnvlib/plots.py +++ b/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(): diff --git a/cnvlib/scatter.py b/cnvlib/scatter.py index a15e4bb6..a339e3a7 100644 --- a/cnvlib/scatter.py +++ b/cnvlib/scatter.py @@ -2,6 +2,7 @@ from __future__ import annotations +import collections import logging from typing import TYPE_CHECKING @@ -9,9 +10,10 @@ 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), diff --git a/doc/plots.rst b/doc/plots.rst index a860ae65..2a6e4a5e 100644 --- a/doc/plots.rst +++ b/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 ``; 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``. diff --git a/test/test_plots.py b/test/test_plots.py index d43ec561..5ba017d5 100644 --- a/test/test_plots.py +++ b/test/test_plots.py @@ -257,6 +257,154 @@ def test_gene_coords_by_name(self): names_seen.update(label.split(",")) self.assertEqual(names_seen, {"ERBB2", "MIR4728"}) + def test_recurrent_name_is_one_region_per_locus(self): + """A name recurring at several loci gets one region per locus. + + The grouping is `CopyNumArray.by_gene`'s, so the regions must equal that + method's groups bin for bin: a placeholder bin is absorbed into the run + it interrupts, while a bin naming another gene ends it. The superseded + rule reported a single span from the first occurrence to the last, + swallowing MIDGENE. + """ + cnarr = cnary.CopyNumArray.from_rows( + [ + ["chr1", 1000, 2000, "REPEAT", 0.0], + ["chr1", 2000, 2500, "-", 0.0], + ["chr1", 2500, 3000, "REPEAT", 0.0], + ["chr1", 3000, 4000, "MIDGENE", 0.0], + ["chr1", 9000, 9500, "REPEAT", 0.0], + ] + ) + coords = plots.gene_coords_by_name(cnarr, ["REPEAT"]) + self.assertEqual( + coords["chr1"], [(1000, 3000, "REPEAT"), (9000, 9500, "REPEAT")] + ) + by_gene_spans = [ + (rows.start.iat[0], rows.end.iat[-1]) + for gene, rows in cnarr.by_gene() + if gene == "REPEAT" + ] + self.assertEqual([(s, e) for s, e, _name in coords["chr1"]], by_gene_spans) + + def test_recurrent_name_on_wgs_fixture(self): + """Every locus of a real repeat family is reported, and only those. + + `Y_RNA` occurs at 29 loci across chr17 of this fixture; the superseded + rule reported one 76 Mb region reaching across all of them, which is + what `scatter -g Y_RNA` highlighted. + """ + cnarr = cnvlib.read("formats/wgs-chr17.cnr") + coords = plots.gene_coords_by_name(cnarr, ["Y_RNA"]) + by_gene_spans = [ + (rows.start.iat[0], rows.end.iat[-1]) + for gene, rows in cnarr.by_gene() + if gene == "Y_RNA" + ] + self.assertEqual(len(by_gene_spans), 29) + self.assertEqual([(s, e) for s, e, _name in coords["chr17"]], by_gene_spans) + + def test_name_on_two_chromosomes_reports_both(self): + """A name reused on another chromosome yields a region on each. + + The superseded rule asserted the name was unique to one chromosome and + died before its caller could report the ambiguity. + """ + cnarr = cnary.CopyNumArray.from_rows( + [ + ["chr1", 1000, 2000, "DUP", 0.0], + ["chr2", 5000, 6000, "DUP", 0.0], + ] + ) + self.assertEqual( + plots.gene_coords_by_name(cnarr, ["DUP"]), + {"chr1": [(1000, 2000, "DUP")], "chr2": [(5000, 6000, "DUP")]}, + ) + + def test_unknown_gene_name_raises(self): + """A requested name absent from the bins is an error, not an omission.""" + cnarr = cnary.CopyNumArray.from_rows([["chr1", 1000, 2000, "REAL", 0.0]]) + with self.assertRaises(ValueError): + plots.gene_coords_by_name(cnarr, ["REAL", "ABSENT"]) + + +class ScatterSelectionTests(unittest.TestCase): + """`scatter`'s region and gene selection (scatter.select_range_genes).""" + + @staticmethod + def _cnarr(): + """Bins with REPEAT at two chr1 loci and DUP on chr1 and chr2.""" + return cnary.CopyNumArray.from_rows( + [ + ["chr1", 1000, 2000, "REPEAT", 0.0], + ["chr1", 2000, 3000, "REPEAT", 0.0], + ["chr1", 3000, 4000, "DUP", 0.0], + ["chr1", 9000, 9500, "REPEAT", 0.0], + ["chr2", 5000, 6000, "DUP", 0.0], + ["chr2", 7000, 8000, "OTHER", 0.0], + ] + ) + + def _select(self, show_range, show_gene, window_width=500): + _probes, _segs, _snvs, _loh, _som, window, gene_ranges, chrom = ( + scatter.select_range_genes( + self._cnarr(), None, None, show_range, show_gene, window_width + ) + ) + return window, gene_ranges, chrom + + def test_each_locus_is_highlighted(self): + """Each locus of a repeated name is highlighted, and the window spans all.""" + with self.assertLogs(level="WARNING") as cm: + window, gene_ranges, chrom = self._select(None, "REPEAT") + self.assertEqual(chrom, "chr1") + self.assertEqual(gene_ranges, [(1000, 3000, "REPEAT"), (9000, 9500, "REPEAT")]) + self.assertEqual(window, (500, 10000)) + self.assertIn("2 separate loci", "".join(cm.output)) + + def test_region_keeps_only_the_loci_inside_it(self): + """An explicit region selects among a repeated name's loci.""" + window, gene_ranges, chrom = self._select("chr1:1-5000", "REPEAT") + self.assertEqual(chrom, "chr1") + self.assertEqual(gene_ranges, [(1000, 3000, "REPEAT")]) + self.assertEqual(window, (0, 5000)) + + def test_region_without_any_locus_of_the_gene_raises(self): + """A named gene with no locus in the region is the user's mistake.""" + with self.assertRaises(ValueError) as cm: + self._select("chr1:4500-8000", "REPEAT") + self.assertIn("outside specified region", str(cm.exception)) + + def test_name_on_two_chromosomes_needs_a_chromosome(self): + """One panel plots one chromosome, so the ambiguity must be resolved.""" + with self.assertRaises(ValueError) as cm: + self._select(None, "DUP") + message = str(cm.exception) + self.assertIn("chr1", message) + self.assertIn("chr2", message) + self.assertIn("-c", message) + + def test_selected_chromosome_disambiguates_a_reused_name(self): + """`-c` picks which of a reused name's chromosomes is plotted.""" + for chrom_sel, expected in (("chr1", (3000, 4000)), ("chr2", (5000, 6000))): + with self.subTest(chromosome=chrom_sel): + _window, gene_ranges, chrom = self._select(chrom_sel, "DUP") + self.assertEqual(chrom, chrom_sel) + self.assertEqual(gene_ranges, [(*expected, "DUP")]) + + def test_gene_absent_from_selected_chromosome_raises(self): + """A gene named with `-c` elsewhere is reported, not silently dropped.""" + with self.assertRaises(ValueError) as cm: + self._select("chr1", "OTHER") + message = str(cm.exception) + self.assertIn("OTHER", message) + self.assertIn("chr2", message) + + def test_empty_gene_list_highlights_nothing(self): + """`-g ''` plots the region with nothing highlighted (doc/plots.rst).""" + window, gene_ranges, chrom = self._select("chr1:1-5000", "") + self.assertEqual((chrom, window), ("chr1", (0, 5000))) + self.assertEqual(gene_ranges, []) + class DiagramGeneLabelTests(unittest.TestCase): """`diagram` --gene and directional --threshold-low/-high (#248)."""