Select only the data's contig from a multi-contig BED file

docs/dev/pipeline_map.toml

@@ -131,7 +131,7 @@ notes = "region is samtools style with an inclusive end; chrom_start and chrom_e
 [inputs."HaplotypeMatrix.from_zarr"]
 params = ["path", "region", "accessible_bed", "pop_assignment", "streaming", "chunk_bp", "prefetch", "backend", "fields"]
 returns = "HaplotypeMatrix, or StreamingHaplotypeMatrix when streaming='always', or 'auto' and the matrix does not fit"
-span = "accessible_bed only; there is no include_invariant, so n_total_sites stays None."
+span = "accessible_bed only; there is no include_invariant, so n_total_sites stays None. The BED contig is the region's, else the store's (eager and streaming alike)."
 notes = "Reads VCZ and scikit-allel layouts (zarr_io). The streaming path rejects fields=."
 
 [inputs."HaplotypeMatrix.from_ts"]
@@ -153,7 +153,7 @@ notes = "No region and no samples parameter, unlike HaplotypeMatrix.from_vcf."
 [inputs."GenotypeMatrix.from_zarr"]
 params = ["path", "region", "accessible_bed", "include_invariant", "pop_assignment", "streaming", "chunk_bp", "prefetch", "backend", "fields"]
 returns = "GenotypeMatrix, or StreamingGenotypeMatrix on the streaming path"
-span = "include_invariant=True sets n_total_sites to the variant count on the eager path only; the streaming path ignores it. accessible_bed sets the mask on both paths."
+span = "include_invariant=True sets n_total_sites to the variant count on the eager path only; the streaming path ignores it. accessible_bed sets the mask on both paths, on the region's contig, else the store's."
 
 [inputs."GenotypeMatrix.from_haplotype_matrix"]
 params = ["hap_matrix"]
@@ -563,7 +563,7 @@ missing_data = ["include", "exclude"]
 span = "per engine; see [engines]"
 undefined = "per statistic and engine; the scatter and fused engines give a row for a window with no variant, the fallback gives no row"
 parallel = []
-notes = "The dispatcher. It picks one engine per request ([engines]); the streaming path runs it once per chunk. It raises TypeError for a keyword argument that no engine reads (windowed_analysis._check_windowed_kwargs), ValueError for an unknown window_type or missing_data, for regions without window_type='regions' and the reverse, and NotImplementedError for window_type 'snp' or 'regions' on a StreamingHaplotypeMatrix and for names that no engine computes on those window types."
+notes = "The dispatcher. It picks one engine per request ([engines]); the streaming path runs it once per chunk. accessible_bed uses chrom, else the stream's own contig; an eager matrix needs chrom for a BED with more than one contig. It raises TypeError for a keyword argument that no engine reads (windowed_analysis._check_windowed_kwargs), ValueError for an unknown window_type or missing_data, for regions without window_type='regions' and the reverse, and NotImplementedError for window_type 'snp' or 'regions' on a StreamingHaplotypeMatrix and for names that no engine computes on those window types."
 
 [functions."windowed_analysis.windowed_statistics_fused"]
 # signature: (haplotype_matrix, bp_bins, statistics, population, pop1, pop2, per_base, is_accessible, _win_starts, _win_stops, missing_data, chrom)
@@ -1866,7 +1866,7 @@ missing_data = []
 span = "none"
 undefined = "pairs below a population's sample-size floor add 0.0; an empty bin sums to 0.0"
 parallel = ["HaplotypeMatrix.compute_ld_statistics_gpu_single_pop", "HaplotypeMatrix.compute_ld_statistics_gpu_two_pops", "haplotype_matrix._stream_ld_single_pop", "haplotype_matrix._stream_ld_two_pops"]
-notes = "'other' is a VCF path; two-pop haplotype runs use ld_statistics.dd/dz/pi2, all other cases the fused kernels."
+notes = "'other' is a VCF path, loaded with accessible_bed so the BED contig is the VCF's; two-pop haplotype runs use ld_statistics.dd/dz/pi2, all other cases the fused kernels."
 
 [functions."decomposition.pca"]
 # signature: (haplotype_matrix, n_components, population, missing_data)
@@ -2214,7 +2214,7 @@ inputs = ["other"]
 streaming = "n/a"
 missing_data = []
 span = "source intervals for the accessible mask that get_span('auto') counts"
-undefined = "skips comment, track, browser and short lines; chrom=None returns every contig"
+undefined = "skips comment, track, browser and short lines; chrom=None returns every contig; a chrom the BED lacks matches a BED contig by name without chr and case, with a ContigNameWarning"
 parallel = []
 notes = "Input is a plain-text BED path (open() with no gzip support); coordinates stay 0-based half-open."
 
@@ -2225,7 +2225,7 @@ inputs = ["other"]
 streaming = "n/a"
 missing_data = []
 span = "builds the AccessibleMask that get_span('auto') counts"
-undefined = "all-False mask, with no error, when the BED has no interval on chrom"
+undefined = "all-False mask, with no error, when the BED has no interval on chrom; ValueError for chrom=None on a BED with more than one contig; chrom matches by name as in resolve_accessible_mask"
 parallel = []
 notes = "Input is a BED path; mask index k is 1-based position offset + k."
 
@@ -2236,7 +2236,7 @@ inputs = ["other"]
 streaming = "n/a"
 missing_data = []
 span = "builds the AccessibleMask that get_span('auto') counts"
-undefined = "raises ValueError when the BED has no interval on chrom; chrom=None reads every contig"
+undefined = "raises ValueError, listing the BED contigs, when the BED has no interval on chrom, or for chrom=None on a BED with more than one contig (a one-contig BED needs no chrom); a chrom the BED lacks matches a BED contig by name without chr and case, with a ContigNameWarning"
 parallel = ["windowed_analysis._compute_window_bases"]
 notes = "Takes a BED path, a numpy bool array (offset chrom_start, or 0) or an AccessibleMask; a cupy array is not caught and is opened as a path (TypeError)."
 

docs/source/changelog.rst

@@ -514,6 +514,15 @@ Bug fixes
   ``regions`` is not passed alongside ``window_type='regions'`` and vice-versa.
   An error is raised when ``window_type='snp'`` or ``'regions'`` for a
   ``StreamingHaplotypeMatrix``.
+* A BED file with more than one contig was merged into one mask when no
+  contig was named, so other contigs' intervals counted as accessible
+  bases of the data. Eager ``from_zarr`` now takes the store's contig (as
+  streaming does), streaming ``windowed_analysis`` the stream's, and
+  ``compute_ld_statistics`` the VCF's. Elsewhere (``set_accessible_mask``,
+  eager ``windowed_analysis``), a multi-contig BED with no ``chrom`` raises
+  ``ValueError``. A BED that names the contig differently, such as ``chr1``
+  for the data's ``1``, is matched by name without ``chr`` and case, with a
+  ``ContigNameWarning``.
 * ``windowed_analysis`` ignored an unknown keyword argument when the
   scatter or fused engine served the request, so a misspelled name, or
   ``population=`` in place of ``populations=``, changed nothing and gave

pg_gpu/__init__.py

@@ -70,10 +70,10 @@
 from ._warnings import (
     BadlyChunkedWarning, HaploidDataWarning, KvikioUnavailableWarning,
     MemoryLimitedWarning, MultiallelicCapWarning, BiallelicOnlyWarning,
-    UnpairedRowsWarning, FewCompleteHaplotypesWarning,
+    UnpairedRowsWarning, FewCompleteHaplotypesWarning, ContigNameWarning,
 )
 
-__all__ = ['ld_statistics', 'diversity', 'divergence', 'windowed_analysis', 'selection', 'sfs', 'admixture', 'decomposition', 'plotting', 'distance_stats', 'resampling', 'HaplotypeMatrix', 'GenotypeMatrix', 'WindowedAnalyzer', 'windowed_analysis', 'AccessibleMask', 'bed_to_mask', 'parse_bed', 'LocalPCAResult', 'LostructResult', 'pca', 'randomized_pca', 'pca_dosage', 'randomized_pca_dosage', 'local_pca', 'local_pca_jackknife', 'lostruct', 'pc_dist', 'corners', 'block_jackknife', 'block_bootstrap', 'MemoryLimitedWarning', 'HaploidDataWarning', 'BadlyChunkedWarning', 'KvikioUnavailableWarning', 'MultiallelicCapWarning', 'BiallelicOnlyWarning', 'UnpairedRowsWarning', 'FewCompleteHaplotypesWarning']
+__all__ = ['ld_statistics', 'diversity', 'divergence', 'windowed_analysis', 'selection', 'sfs', 'admixture', 'decomposition', 'plotting', 'distance_stats', 'resampling', 'HaplotypeMatrix', 'GenotypeMatrix', 'WindowedAnalyzer', 'windowed_analysis', 'AccessibleMask', 'bed_to_mask', 'parse_bed', 'LocalPCAResult', 'LostructResult', 'pca', 'randomized_pca', 'pca_dosage', 'randomized_pca_dosage', 'local_pca', 'local_pca_jackknife', 'lostruct', 'pc_dist', 'corners', 'block_jackknife', 'block_bootstrap', 'MemoryLimitedWarning', 'HaploidDataWarning', 'BadlyChunkedWarning', 'KvikioUnavailableWarning', 'MultiallelicCapWarning', 'BiallelicOnlyWarning', 'UnpairedRowsWarning', 'FewCompleteHaplotypesWarning', 'ContigNameWarning']
 
 # Version is derived from the git tag at build time (hatch-vcs) and read here
 # from the installed package metadata, so there is no hardcoded string to keep

pg_gpu/_warnings.py

@@ -273,6 +273,21 @@ def _warn_biallelic_only(n_dropped, *, context, stacklevel=3, action="dropped"):
     )
 
 
+class ContigNameWarning(UserWarning):
+    """Emitted when a BED contig is matched to the data's contig by name.
+
+    BED files from UCSC name contigs ``chr1``, ``chr2``, ..., while many VCFs
+    and stores name them ``1``, ``2``, .... When the data's contig is not in
+    the BED, pg_gpu compares the names without a leading ``chr`` and without
+    case; if exactly one BED contig matches, it uses that contig and emits
+    this warning. Silence with::
+
+        import warnings
+        from pg_gpu import ContigNameWarning
+        warnings.filterwarnings("ignore", category=ContigNameWarning)
+    """
+
+
 class FewCompleteHaplotypesWarning(UserWarning):
     """Emitted when many windows have few haplotypes with no missing call.
 

pg_gpu/accessible.py

@@ -4,9 +4,12 @@
 boolean array with offset-aware, O(1) range queries via prefix sums.
 """
 
+import warnings
+
 import cupy as cp
 import numpy as np
 
+from ._warnings import ContigNameWarning, _format_sample_list
 from .zarr_io import parse_region
 
 
@@ -141,15 +144,34 @@ def parse_bed(path, chrom=None):
         Path to a BED file (tab-delimited: chrom, start, end, ...).
         Lines starting with '#', 'track', or 'browser' are skipped.
     chrom : str, optional
-        If provided, only return intervals on this chromosome.
+        If provided, only return intervals on this chromosome. A name the
+        BED lacks matches a BED contig by name without a leading 'chr' and
+        case ('1' and 'chr1'), with a ContigNameWarning.
 
     Returns
     -------
     list of (str, int, int)
         Parsed intervals as (chromosome, start, end) with 0-based half-open
         coordinates.
     """
-    intervals = []
+    return _read_bed(path, chrom)[0]
+
+
+def _contig_key(name):
+    """A contig name without a leading 'chr' and without case, so the UCSC
+    name ('chr1') and the Ensembl or NCBI name ('1') of a contig match."""
+    key = str(name).lower()
+    return key[3:] if key.startswith('chr') else key
+
+
+def _read_bed(path, chrom=None):
+    """``(intervals, contigs)`` from one pass over a BED file: the intervals
+    on ``chrom`` (every interval when it is None), and every contig the BED
+    names, in file order. ``chrom`` matches by name as ``parse_bed`` says."""
+    key = None if chrom is None else _contig_key(chrom)
+    contigs = {}
+    kept = []
+    last, keep = None, False
     with open(path) as f:
         for line in f:
             line = line.strip()
@@ -161,10 +183,39 @@ def parse_bed(path, chrom=None):
                 fields = line.split()
             if len(fields) < 3:
                 continue
-            c, s, e = fields[0], int(fields[1]), int(fields[2])
-            if chrom is None or c == chrom:
-                intervals.append((c, s, e))
-    return intervals
+            c = fields[0]
+            # BED lines come grouped by contig, so the name test runs once
+            # per contig, not once per line.
+            if c != last:
+                last = c
+                contigs[c] = None
+                keep = key is None or _contig_key(c) == key
+            if keep:
+                kept.append((c, int(fields[1]), int(fields[2])))
+    if chrom is not None:
+        chrom = str(chrom)
+        names = list(dict.fromkeys(c for c, _, _ in kept))
+        if chrom not in names and len(names) == 1:
+            warnings.warn(
+                f"BED contig {names[0]!r} does not match data contig "
+                f"{chrom!r} exactly; inferring that they are the same contig.",
+                ContigNameWarning, stacklevel=3)
+            chrom = names[0]
+        kept = [iv for iv in kept if iv[0] == chrom]
+    return kept, list(contigs)
+
+
+def _contig_intervals(path, chrom):
+    """``_read_bed(path, chrom)``, but a BED with more than one contig raises
+    ValueError when ``chrom`` is None: merging them would count other
+    contigs' intervals as accessible bases of the data."""
+    intervals, contigs = _read_bed(path, chrom)
+    if chrom is None and len(contigs) > 1:
+        raise ValueError(
+            f"BED file {path} has intervals on {len(contigs)} contigs "
+            f"({_format_sample_list(contigs, range(len(contigs)), limit=5)}); "
+            "name the contig of the data with chrom= or region=.")
+    return intervals, contigs
 
 
 def bed_to_mask(path, chrom, length, offset=0):
@@ -174,8 +225,9 @@ def bed_to_mask(path, chrom, length, offset=0):
     ----------
     path : str or path-like
         Path to a BED file defining accessible regions.
-    chrom : str
-        Chromosome name to extract from the BED file.
+    chrom : str or None
+        Chromosome name to extract from the BED file. None is allowed only
+        for a BED with one contig.
     length : int
         Total length of the mask array (number of genomic positions to cover).
     offset : int
@@ -193,7 +245,12 @@ def bed_to_mask(path, chrom, length, offset=0):
     positions s+1 through e (inclusive). Mask index k represents 1-based
     position offset + k.
     """
-    intervals = parse_bed(path, chrom=chrom)
+    return _intervals_to_mask(_contig_intervals(path, chrom)[0], length, offset)
+
+
+def _intervals_to_mask(intervals, length, offset):
+    """AccessibleMask of ``length`` bases from ``offset`` with ``intervals``
+    (0-based half-open BED records) set True."""
     mask = np.zeros(length, dtype=bool)
     for _, s, e in intervals:
         # BED (s, e) [0-based half-open] = 1-based positions s+1..e
@@ -218,7 +275,8 @@ def resolve_accessible_mask(mask_or_path, chrom_start, chrom_end, chrom=None):
         silently dropping BED accessible bases that fall outside the
         matrix's variant range (relevant for variants-only VCFs).
     chrom : str, optional
-        Chromosome name (required for BED file input).
+        Chromosome name. Required for a BED file with more than one contig;
+        a one-contig BED needs none.
 
     Returns
     -------
@@ -232,10 +290,12 @@ def resolve_accessible_mask(mask_or_path, chrom_start, chrom_end, chrom=None):
     # Treat as file path. Build the mask over the union of [chrom_start,
     # chrom_end] and the BED's own extent so no accessible bases are lost
     # when the matrix's variant range is narrower than the BED.
-    intervals = parse_bed(mask_or_path, chrom=chrom)
+    intervals, contigs = _contig_intervals(mask_or_path, chrom)
     if not intervals:
+        where = "" if chrom is None else f" for contig {str(chrom)!r}"
         raise ValueError(
-            f"No BED intervals found for chrom={chrom!r} in {mask_or_path}")
+            f"No BED intervals found{where} in {mask_or_path}; its contigs "
+            f"are: {_format_sample_list(contigs, range(len(contigs)), limit=5) or 'none'}.")
     # 1-based inclusive bounds of the BED (BED is 0-based half-open: a record
     # (s, e) covers 1-based positions s+1..e).
     bed_min = min(s + 1 for _, s, _ in intervals)
@@ -248,7 +308,7 @@ def resolve_accessible_mask(mask_or_path, chrom_start, chrom_end, chrom=None):
     if length <= 0:
         raise ValueError(
             "Could not determine mask range from BED and chromosome bounds")
-    return bed_to_mask(mask_or_path, chrom=chrom, length=length, offset=offset)
+    return _intervals_to_mask(intervals, length, offset)
 
 
 def slice_fields(fields, keep_idx, accessible_idx=None):

pg_gpu/genotype_matrix.py

@@ -219,7 +219,9 @@ def set_accessible_mask(self, mask_or_path, chrom=None):
         mask_or_path : str, path-like, numpy.ndarray, or AccessibleMask
             BED file path, boolean array, or AccessibleMask instance.
         chrom : str, optional
-            Chromosome name (required for BED file input).
+            Chromosome name to extract from the BED file. Required when the
+            BED file has more than one contig (otherwise ValueError); a
+            one-contig BED needs none.
         """
         self.accessible_mask = resolve_accessible_mask(
             mask_or_path, self.chrom_start, self.chrom_end, chrom)
@@ -636,7 +638,7 @@ def _build_eager(cls, path, *, region, accessible_bed,
 
         n_total_sites = gt.shape[0] if include_invariant else None
 
-        chrom, start, stop = parse_region(region)
+        _, start, stop = parse_region(region)
         start = start if start is not None else int(positions[0])
         stop = stop - 1 if stop is not None else int(positions[-1])
 
@@ -658,7 +660,7 @@ def _build_eager(cls, path, *, region, accessible_bed,
                 region=region,
             )
         if accessible_bed is not None:
-            gm.set_accessible_mask(accessible_bed, chrom=chrom)
+            gm.set_accessible_mask(accessible_bed, chrom=data['chrom'])
 
         try:
             import zarr

pg_gpu/haplotype_matrix.py

@@ -435,8 +435,9 @@ def set_accessible_mask(self, mask_or_path, chrom=None):
             equal to chrom_start. If already an AccessibleMask, assigned
             directly.
         chrom : str, optional
-            Chromosome name to extract from the BED file. Required when
-            mask_or_path is a BED file path.
+            Chromosome name to extract from the BED file. Required when the
+            BED file has more than one contig (otherwise ValueError); a
+            one-contig BED needs none.
 
         Returns
         -------
@@ -787,7 +788,7 @@ def _build_eager(cls, path, *, region, accessible_bed, pop_assignment,
         check_diploid_encoding(gt, sample_names=sample_names,
                                source=f"zarr store '{path}'")
 
-        chrom, start, stop = parse_region(region)
+        _, start, stop = parse_region(region)
         start = start if start is not None else int(positions[0])
         stop = stop - 1 if stop is not None else int(positions[-1])
         hm = build_haplotype_matrix(
@@ -809,7 +810,7 @@ def _build_eager(cls, path, *, region, accessible_bed, pop_assignment,
             )
 
         if accessible_bed is not None:
-            hm.set_accessible_mask(accessible_bed, chrom=chrom)
+            hm.set_accessible_mask(accessible_bed, chrom=data['chrom'])
 
         from .zarr_io import normalize_pop_input
         try:

pg_gpu/moments_ld.py

@@ -173,10 +173,10 @@ def compute_ld_statistics(
                 raise ValueError("pop_file is required when loading from VCF")
             if report:
                 print(f"Loading {vcf_file} (genotypes) ...")
-            gm = GenotypeMatrix.from_vcf(vcf_file)
+            # from_vcf takes the contig from the VCF, so a multi-contig
+            # BED selects the right intervals.
+            gm = GenotypeMatrix.from_vcf(vcf_file, accessible_bed=accessible_bed)
             gm.load_pop_file(pop_file, pops=pops)
-            if accessible_bed is not None and not gm.has_accessible_mask:
-                gm.set_accessible_mask(accessible_bed)
             gm.transfer_to_gpu()
         mat = gm
         if report:
@@ -194,10 +194,8 @@ def compute_ld_statistics(
                 raise ValueError("pop_file is required when loading from VCF")
             if report:
                 print(f"Loading {vcf_file} ...")
-            hm = HaplotypeMatrix.from_vcf(vcf_file)
+            hm = HaplotypeMatrix.from_vcf(vcf_file, accessible_bed=accessible_bed)
             hm.load_pop_file(pop_file, pops=pops)
-            if accessible_bed is not None and not hm.has_accessible_mask:
-                hm.set_accessible_mask(accessible_bed)
             hm.transfer_to_gpu()
         mat = hm
         if report: