diff --git a/dinf/dinf_model.py b/dinf/dinf_model.py index 9b96aeb8..5b403713 100644 --- a/dinf/dinf_model.py +++ b/dinf/dinf_model.py @@ -14,9 +14,10 @@ from .misc import pytree_equal, pytree_shape, pytree_dtype FeatureCollection = Union[np.ndarray, Dict[str, np.ndarray]] -FeatureCollection.__doc__ = ( - """A feature array or a labelled collection of feature arrays.""" -) +#This bit is commented out to be compatible with Python 3.14: +#FeatureCollection.__doc__ = ( +# """A feature array or a labelled collection of feature arrays.""" +#) def _sim_shim(args, *, func, keys): diff --git a/dinf/feature_extractor.py b/dinf/feature_extractor.py index e11fc0bf..62777b1b 100644 --- a/dinf/feature_extractor.py +++ b/dinf/feature_extractor.py @@ -714,8 +714,11 @@ def from_vcf( require_phased=self._global_phased, rng=rng, ) + sample_indices = { + pop: [vb.sample_index_map[id] for id in vb.samples[pop]] + for pop in list(vb.samples) + } num_samples = [len(v) for v in vb.samples.values()] - offsets = np.cumsum([0] + num_samples[:-1]) labelled_indexes = {} ac0 = np.zeros(len(positions)) @@ -728,8 +731,8 @@ def from_vcf( ) # Subsample individuals. - idx = offsets[j] + rng.choice( - num_samples[j], size=self._num_individuals[label], replace=False + idx = rng.choice( + sample_indices[label], size=self._num_individuals[label], replace=False ) labelled_indexes[label] = idx H = G[:, idx, : self._ploidy[label]] diff --git a/dinf/plot.py b/dinf/plot.py index 8b9f07d5..6a326090 100644 --- a/dinf/plot.py +++ b/dinf/plot.py @@ -1189,7 +1189,7 @@ def __call__(self, args: argparse.Namespace): ax.set_xlabel("Iteration") ax.xaxis.set_major_locator(matplotlib.ticker.MaxNLocator(integer=True)) - if parameters is not None and x_param != "_Pr": + if parameters is not None and x_param != "_Pr" and parameters[x_param].truth is not None: truth = parameters[x_param].truth ax.axhline(truth, c="red", ls="-", alpha=0.7) @@ -1200,7 +1200,7 @@ def __call__(self, args: argparse.Namespace): def main(args_list=None): top_parser = argparse.ArgumentParser( - prog="dinf.plot", description="Dinf plotting tools." + prog="dinf-plot", description="Dinf plotting tools." ) top_parser.add_argument( "-V", "--version", action="version", version=dinf.__version__ diff --git a/dinf/tabulate.py b/dinf/tabulate.py index d988f0d2..033fd0a5 100644 --- a/dinf/tabulate.py +++ b/dinf/tabulate.py @@ -217,7 +217,7 @@ def __call__(self, args: argparse.Namespace): def main(args_list=None): top_parser = argparse.ArgumentParser( - prog="dinf.tabulate", description="Tabulate Dinf output." + prog="dinf-tabulate", description="Tabulate Dinf output." ) top_parser.add_argument( "-V", "--version", action="version", version=dinf.__version__ diff --git a/dinf/vcf.py b/dinf/vcf.py index b47ee2cf..a84f2d1c 100644 --- a/dinf/vcf.py +++ b/dinf/vcf.py @@ -24,7 +24,7 @@ def get_samples_from_1kgp_metadata(filename: str, /, *, populations: list) -> di :return: A dict mapping population names to a list of sample IDs. """ - data = np.recfromtxt(filename, names=True, encoding="ascii") + data = np.genfromtxt(filename, names=True, dtype=[' 2. # https://github.com/brentp/cyvcf2/issues/227 @@ -128,7 +129,16 @@ def get_genotype_matrix( if not variant.is_snp: continue - a = variant.genotype.array() + try: + a = variant.genotype.array() + except AttributeError as e: + n_empty += 1 + if n_empty >= 10: + raise ValueError( + f"At least 10 variants from region {chrom}:{start:d}-{end:d} failed to return valid genotypes, cannot proceed. Did any samples from your metadata get detected in the VCFs?" + ) from e + else: + continue gt = a[:, :-1] # Check phasing. For ploidy == 1, genotypes are reported as unphased. @@ -209,6 +219,17 @@ class BagOfVcf(collections.abc.Mapping): that all individuals in the VCF will be treated as exchangeable. """ + sample_index_map: collections.abc.Mapping[str, int] + """ + A dictionary that maps an individual name to a genotype matrix column index. + + The individual names correspond to the VCF columns for which genotypes + will be sampled, and the genotype matrix column indices correspond to + the genotype matrix containing only the samples requested. This map + enables random subsampling of individuals from a population where the + individuals may not be in contiguous population blocks in the VCF. + """ + def __init__( self, files: Iterable[str | pathlib.Path], @@ -311,6 +332,10 @@ def _fill_bag( f"Requested individuals not found in {file}: " f"{', '.join(individuals_not_found)}." ) + #This gets overwritten several times, but that's okay: + self.sample_index_map = { + id: i for i, id in enumerate(vcf.samples) + } if contigs is None: try: vcf.seqlens @@ -360,6 +385,9 @@ def _fill_bag( # The VCF objects are not valid in child processes, so we set a flag # before forking to ensure they're reopened when used in the child. + #Note: Given 'spawn' is used instead of 'fork', this doesn't actually + # have any impact, but could be part of the code path if switched + # back to 'fork'. self._needs_reopen = False os.register_at_fork(before=self._close) diff --git a/examples/gutenkunst2009/model.py b/examples/gutenkunst2009/model.py index bc1ba12a..9f7e9446 100644 --- a/examples/gutenkunst2009/model.py +++ b/examples/gutenkunst2009/model.py @@ -14,7 +14,7 @@ ) contig_lengths = dinf.get_contig_lengths( "GRCh38_full_analysis_set_plus_decoy_hla.fa.fai", - keep_contigs={f"chr{c + 1}" for c in range(21)}, # Exclude chrX, etc. + keep_contigs={f"chr{c + 1}" for c in range(22)}, # Exclude chrX, etc. ) num_individuals = 64 recombination_rate = 1.25e-8 diff --git a/examples/one_population_growth/model.py b/examples/one_population_growth/model.py index cd1bfd99..60fd62d1 100644 --- a/examples/one_population_growth/model.py +++ b/examples/one_population_growth/model.py @@ -12,11 +12,11 @@ population_name = "CEU" samples = dinf.get_samples_from_1kgp_metadata( - "20130606_g1k_3202_samples_ped_population.txt", [population_name] + "20130606_g1k_3202_samples_ped_population.txt", populations=[population_name] ) contig_lengths = dinf.get_contig_lengths( "GRCh38_full_analysis_set_plus_decoy_hla.fa.fai", - keep_contigs={f"chr{c + 1}" for c in range(21)}, # Exclude chrX, etc. + keep_contigs={f"chr{c + 1}" for c in range(22)}, # Exclude chrX, etc. ) recombination_rate = 1.25e-8 mutation_rate = 1.25e-8 diff --git a/examples/post_ooa/model.py b/examples/post_ooa/model.py index 76cacf97..5807839e 100644 --- a/examples/post_ooa/model.py +++ b/examples/post_ooa/model.py @@ -12,11 +12,11 @@ populations = ["CEU", "CHB"] samples = dinf.get_samples_from_1kgp_metadata( - "20130606_g1k_3202_samples_ped_population.txt", populations + "20130606_g1k_3202_samples_ped_population.txt", populations=populations ) contig_lengths = dinf.get_contig_lengths( "GRCh38_full_analysis_set_plus_decoy_hla.fa.fai", - keep_contigs={f"chr{c + 1}" for c in range(21)}, # Exclude chrX, etc. + keep_contigs={f"chr{c + 1}" for c in range(22)}, # Exclude chrX, etc. ) recombination_rate = 1.25e-8 mutation_rate = 1.25e-8 diff --git a/examples/two_pop_ooa/model.py b/examples/two_pop_ooa/model.py index cabbf0e6..955b0f86 100644 --- a/examples/two_pop_ooa/model.py +++ b/examples/two_pop_ooa/model.py @@ -12,11 +12,11 @@ populations = ["YRI", "CEU"] samples = dinf.get_samples_from_1kgp_metadata( - "20130606_g1k_3202_samples_ped_population.txt", populations + "20130606_g1k_3202_samples_ped_population.txt", populations=populations ) contig_lengths = dinf.get_contig_lengths( "GRCh38_full_analysis_set_plus_decoy_hla.fa.fai", - keep_contigs={f"chr{c + 1}" for c in range(21)}, # Exclude chrX, etc. + keep_contigs={f"chr{c + 1}" for c in range(22)}, # Exclude chrX, etc. ) recombination_rate = 1.25e-8 mutation_rate = 1.25e-8 diff --git a/tests/conftest.py b/tests/conftest.py index ea044953..94647a26 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -5,7 +5,7 @@ @pytest.fixture(scope="session") -@pytest.mark.usefixtures("tmp_path_factory") +#@pytest.mark.usefixtures("tmp_path_factory") def discriminator_file(tmp_path_factory): discriminator_file = tmp_path_factory.mktemp("discr") / "discriminator.nn" ex = "examples/bottleneck/model.py" @@ -28,8 +28,8 @@ def discriminator_file(tmp_path_factory): @pytest.fixture(scope="session") -@pytest.mark.usefixtures("tmp_path_factory") -@pytest.mark.usefixtures("discriminator_file") +#@pytest.mark.usefixtures("tmp_path_factory") +#@pytest.mark.usefixtures("discriminator_file") def data_file(tmp_path_factory, discriminator_file): data_file = tmp_path_factory.mktemp("data") / "data.npz" ex = "examples/bottleneck/model.py" @@ -51,7 +51,7 @@ def data_file(tmp_path_factory, discriminator_file): @pytest.fixture(scope="session") -@pytest.mark.usefixtures("tmp_path_factory") +#@pytest.mark.usefixtures("tmp_path_factory") def mc_outdir(tmp_path_factory): output_folder = tmp_path_factory.mktemp("out") ex = "examples/bottleneck/model.py" diff --git a/tests/test_vcf.py b/tests/test_vcf.py index 3b80c8bb..00002d97 100644 --- a/tests/test_vcf.py +++ b/tests/test_vcf.py @@ -19,6 +19,10 @@ FamilyID SampleID FatherID MotherID Sex Population Superpopulation HG00096 HG00096 0 0 1 GBR EUR HG00097 HG00097 0 0 2 GBR EUR +HG00180 HG00180 0 0 2 FIN EUR +HG00181 HG00181 0 0 1 FIN EUR +HG00234 HG00234 0 0 1 GBR EUR +HG00235 HG00235 0 0 2 GBR EUR SH001 HG00403 0 0 1 CHS EAS SH001 HG00404 0 0 2 CHS EAS SH001 HG00405 HG00403 HG00404 2 CHS EAS @@ -32,13 +36,15 @@ def test_simple(self, tmp_path): with open(filename, "w") as f: f.write(_1kgp_test_metadata) gbr = dinf.get_samples_from_1kgp_metadata(filename, populations=["GBR"]) - assert gbr == {"GBR": ["HG00096", "HG00097"]} + assert gbr == {"GBR": ["HG00096", "HG00097", "HG00234", "HG00235"]} + fin = dinf.get_samples_from_1kgp_metadata(filename, populations=["FIN"]) + assert fin == {"FIN": ["HG00180", "HG00181"]} chs = dinf.get_samples_from_1kgp_metadata(filename, populations=["CHS"]) assert chs == {"CHS": ["HG00403", "HG00404"]} samples = dinf.get_samples_from_1kgp_metadata( - filename, populations=["GBR", "CHS"] + filename, populations=["GBR", "FIN", "CHS"] ) - assert samples == dict(**gbr, **chs) + assert samples == dict(**gbr, **fin, **chs) class TestGetContigLengths: