Repository navigation
Description
Since #401, HaplotypeMatrix.from_ts gives phased=True, but a matrix saved with to_zarr and read back with from_zarr has phased=False. Phase-dependent statistics on it then emit UnphasedDataWarning, although the data are the same phased haplotypes. A common workflow (simulate with msprime, save to VCZ, analyze later) hits this.
On main (53ea50b):
import tempfile, warnings
import msprime
from pg_gpu import HaplotypeMatrix, selection
ts = msprime.sim_mutations(
msprime.sim_ancestry(10, sequence_length=1e5, recombination_rate=1e-8,
population_size=1e4, random_seed=1),
rate=1e-8, random_seed=1)
hm = HaplotypeMatrix.from_ts(ts)
print(hm.phased) # True
path = tempfile.mkdtemp() + "/sim.vcz"
hm.to_zarr(path, format="vcz", contig_name="1")
back = HaplotypeMatrix.from_zarr(path)
print(back.phased) # False
with warnings.catch_warnings(record=True) as w:
warnings.simplefilter("always")
selection.garud_h(back)
print([type(x.message).__name__ for x in w]) # ['UnphasedDataWarning']Cause
- The VCZ writer (
zarr_io.write_vcz, called byto_zarr) does not writecall_genotype_phased, the VCZ field that records phase per call. from_zarrdoes not read it. Itsphasedargument defaults toFalse.
Proposed fix
-
to_zarrwritescall_genotype_phased, filled with the matrix'sphasedvalue (all True or all False). bio2zarr writes the same field from a VCF's|and/. -
On load, use the field. Two options:
- a.
from_zarr(phased=None)by default:Nonemeans phased only when the store hascall_genotype_phasedand every call is phased; an explicitTrueorFalsestill wins. - b. Keep
phased=Falseas the default, and read the field only when the user passesphased='auto'.
Option a keeps the round trip without user action. It also makes a VCZ store converted from a fully phased VCF load as phased, which goes a step past Phased attribute to HaplotypeMatrix and warnings in phase-dependent statistics. #401's rule that loaders do not inspect the input's phase (see the discussion on pg-gpu needs to know if HaplotypeMatrix is phased or not #168: partial phasing, for example from GATK, makes inference risky). Requiring every call to be phased keeps a partly phased store unphased. Option b keeps Phased attribute to HaplotypeMatrix and warnings in phase-dependent statistics. #401's rule and needs the user to ask for the store's phase.
- a.
-
A test:
from_ts, thento_zarr, thenfrom_zarr(eager and streaming) keepsphased=True, and the same round trip keepsphased=Falsefor an unphased matrix.
GenotypeMatrix.to_zarr has no flag to write, so it is not affected.
Activity