[andrewkern]

Closes #363.

With no contig named, a BED file with more than one contig went into one mask. So the other contigs' intervals counted as accessible bases of the data, and variants at those positions stayed in the matrix, with no error. A genome-wide BED is the usual input. This affected set_accessible_mask, windowed_analysis(accessible_bed=) (eager and streaming), eager from_zarr(accessible_bed=) with no region, and compute_ld_statistics(vcf_file=, accessible_bed=).

Change

• accessible.py: a BED with more than one contig and chrom=None raises ValueError in resolve_accessible_mask and bed_to_mask (one helper, _contig_intervals). A one-contig BED still needs no chrom.
• zarr_io.py: each reader returns the contig it read ('chrom'): the store's contig for VCZ, the region's contig for the scikit-allel layouts.
• Eager HaplotypeMatrix.from_zarr and GenotypeMatrix.from_zarr use that contig, as the streaming loader already does. So eager and streaming now see the same accessible bases, which missing_data.rst promises.
• Streaming windowed_analysis uses the stream's own contig when chrom is not given. Eager windowed_analysis now documents chrom, and needs it for a multi-contig BED.
• compute_ld_statistics loads the VCF with accessible_bed, so from_vcf applies the VCF's contig.
• resolve_accessible_mask parses the BED once; it parsed it a second time inside bed_to_mask.
• Contig names that differ. A BED that names the contig differently from the data (UCSC chr1 for the store's or VCF's 1) is matched by name without a leading chr and without case. If exactly one BED contig matches, it is used with a new ContigNameWarning (exported from pg_gpu):
  BED contig 'chr1' does not match data contig '1' exactly; inferring that they are the same contig.
  With no match, the error lists the BED's contigs. The rule sits in one reading pass that parse_bed and the mask builders share, so parse_bed(path, chrom=) agrees with them. Without it, the loaders' new contig would break calls that worked on main with a one-contig chr1 BED.
• A flat scikit-allel store read without a region takes its contig from variants/CHROM when the column holds one value.
• Docstrings, the changelog (Bug fixes) and the pipeline map rows are updated.

utils/genome_scan.py builds its own mask from every BED line for tree-sequence input, so it merges contigs the same way. It is a utility script outside the library, so this PR leaves it unchanged.

Tests

tests/test_multi_contig_bed.py, with a BED whose second contig covers what the data's contig leaves out:

• resolve_accessible_mask raises without chrom for two contigs, and works for one contig or with chrom;

• eager from_zarr (both matrix classes) gives the same accessible bases as streaming and as a one-contig BED (50000; main gives 99886);

• eager windowed_analysis raises without chrom, and with chrom matches the one-contig BED;

• streaming windowed_analysis with no chrom matches the eager one-contig result;

• compute_ld_statistics from a VCF gives the same sums with the two-contig and the one-contig BED (genotype and haplotype paths).

• the name rule: parse_bed and resolve_accessible_mask match chr1 to 1 with the warning, and an unmatched name lists the BED's contigs; eager and streaming from_zarr (both classes, one- and two-contig chr BEDs), eager and streaming windowed_analysis, and compute_ld_statistics from a VCF each give the same result with a chr-named BED as with the matching BED;

• a flat scikit-allel store with a one-value CHROM column needs no region.

The tests fail on main and pass here. The masks equal main's byte for byte on 66 edge cases (empty BED, overlapping intervals, intervals outside the range, and others). The full suite passes on a GPU: 2341 passed, 81 skipped, 15 xfailed. Lint passes.

[nspope]

Looks good generally, a couple nits.

One blind spot is that HaplotypeMatrix.from_vcf/GenotypeMatrix.from_vcf infer the chromosome from the first variant's field, for example haplotype_matrix.py:622-623 (there's an equivalent at genotype_matrix.py:535-536):

elif 'variants/CHROM' in vcf:
    chrom = vcf['variants/CHROM'][0]   # just the FIRST variant's contig, never checked

thus there's no inference that the VCF is single-contig. For example,

import numpy as np
from pg_gpu import GenotypeMatrix

vcf_text = """
##fileformat=VCFv4.2
##contig=<ID=1,length=1000>
##contig=<ID=2,length=1000>
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\ts1\ts2\ts3\ts4
1\t100\t.\tA\tT\t.\tPASS\t.\tGT\t0/1\t0/0\t1/1\t0/1
1\t200\t.\tA\tT\t.\tPASS\t.\tGT\t0/0\t0/1\t0/1\t1/1
2\t600\t.\tA\tT\t.\tPASS\t.\tGT\t1/1\t0/1\t0/0\t0/1
2\t700\t.\tA\tT\t.\tPASS\t.\tGT\t0/1\t1/1\t0/1\t0/0
"""
with open("mc.vcf", "w") as f:
    f.write(vcf_text)

# contig 1 accessible only 0-150 (site at 200 is inaccessible);
# contig 2 accessible only 0-650 (site at 700 is inaccessible).
with open("mc.bed", "w") as f:
    f.write("1\t0\t150\n2\t0\t650\n")

gm = GenotypeMatrix.from_vcf("mc.vcf", accessible_bed="mc.bed")
print(gm.chrom_start, gm.chrom_end)        # "100 700" (spans both contigs already, a separate pre-existing quirk of from_vcf unrelated to this PR)
print(gm.accessible_mask.total_accessible) # "150" ... so contig 2's 0-650 interval is silently dropped

One consequence is that compute_ld_statistics routes its VCF through hm/gm.from_vcf(..., accessible_bed=...), avoiding set_accessible_mask, so:

from pg_gpu.moments_ld import compute_ld_statistics

pop_path = "pops.tsv"
with open(pop_path, "w") as f:
    f.write("sample\tpop\ns1\tpop0\ns2\tpop0\ns3\tpop0\ns4\tpop0\n")

# contig 1 accessible only 0-150; contig 2 fully accessible 0-1000.
with open("both.bed", "w") as f:
    f.write("1\t0\t150\n2\t0\t1000\n")

# Same, but contig 2's interval entirely removed.
with open("one.bed", "w") as f:
    f.write("1\t0\t150\n")

kw = dict(vcf_file="mc.vcf", pop_file=pop_path, pops=["pop0"],
          bp_bins=[0, 1000], use_genotypes=True, report=False)

got_both = compute_ld_statistics(accessible_bed="both.bed", **kw)
got_one = compute_ld_statistics(accessible_bed="one.bed", **kw)
print(got_both["sums"])  # [array([0., 0., 0.]), array([0.57142857])]
print(got_one["sums"])   # [array([0., 0., 0.]), array([0.57142857])]
# contig 2's entire accessible interval made no difference to the
# result; it was silently ignored

Note that in this example, the function is not only using sites across both contigs (and ignoring the contig 2 accessibility mask) but is also actually pairing sites across the two contigs... this seems like a deeper issue (recorded in #398) that's outside the scope of the fix here.

[nspope]

Hmm... say you're loading from zarr. Neither hm.from_zarr nor gm.from_zarr has a chrom parameter, so the suggestion in the error message won't work. For example,

from pg_gpu import HaplotypeMatrix

# vcz_store spans 2 contigs; bed has intervals on 2 contigs; no region given.
hm = HaplotypeMatrix.from_zarr(vcz_store, accessible_bed=two_contig_bed)
# ValueError: BED file ... has intervals on 2 contigs (...);
#             name the contig of the data with chrom= or region=.

hm = HaplotypeMatrix.from_zarr(vcz_store, accessible_bed=two_contig_bed, chrom="1")  # following the hint
# TypeError: from_zarr() got an unexpected keyword argument 'chrom'

so that needs special casing to the error or adding a chrom arg for from_zarr

[nspope]

is the stacklevel right here? It's correct for parse_bed(path, chrom=...) -> _read_bed, but there's other entry points that are deeper. For example, resolve_accessible_mask(...) -> _contig_intervals -> read_bed needs stacklevel=4. I'm not clear on what the intended behaviour is here.