Metadata-Version: 2.4
Name: genoray
Version: 3.4.0
Requires-Dist: seqpro>=0.21.1,<0.23
Requires-Dist: numpy>=1.26
Requires-Dist: pandas>=2.2.3
Requires-Dist: hirola>=0.3.0
Requires-Dist: pgenlib>=0.91.0
Requires-Dist: cyvcf2>=0.31.1
Requires-Dist: pysam>=0.22
Requires-Dist: polars>=1.37.1
Requires-Dist: polars-bio>=0.20.1,<0.34
Requires-Dist: pyranges>=0.1.3
Requires-Dist: rich>=13
Requires-Dist: typing-extensions>=4.14
Requires-Dist: pyarrow>=21
Requires-Dist: tqdm>=4.65
Requires-Dist: phantom-types>=3
Requires-Dist: more-itertools>=10
Requires-Dist: loguru>=0.7.0
Requires-Dist: attrs
Requires-Dist: awkward
Requires-Dist: numba
Requires-Dist: cyclopts
Requires-Dist: zstandard
Requires-Dist: pydantic
Requires-Dist: oxbow>=0.5.1,<0.6
Requires-Dist: joblib>=1.4.2,<2
Requires-Dist: joblib-progress>=1.0.6,<2
Requires-Dist: filelock>3,<4
Requires-Dist: scipy>=1.10
Requires-Dist: pooch>=1.7
License-File: LICENSE.txt
Summary: Add your description here
Author-email: David Laub <dlaub@ucsd.edu>
Requires-Python: >=3.10, <3.15
Description-Content-Type: text/markdown; charset=UTF-8; variant=GFM

# `genoray`

If you want to use NumPy with genetic variant data, `genoray` is for you! `genoray` enables ergonomic and efficient range queries of genotypes and dosages from VCF and PGEN (PLINK 2.0) files. `genoray` is also fully type-safe and has minimal dependencies.

# Summary

The `genoray` API more-or-less boils down to just two classes and up to five methods:

- `VCF` and `PGEN` classes for reading VCF and PGEN files, respectively.
- `read` read variants for a single range.
- `chunk` read variants for a single range in chunks.
- `read_ranges` read multiple ranges of variants at once.
- `chunk_ranges` read multiple ranges of variants in chunks.
- `set_samples` subset and/or re-order the samples.

The other important arguments to know are `mode` (and `phasing` for VCF) to set the return type and `max_mem` for chunking. The modes that are available for each file format are always accessible from the class itself, e.g. `VCF.Genos16`, `PGEN.GenosDosages`, etc. You can also filter variants on the fly using the `filter` argument to class constructors.

Also included:

- `SparseVar` and `SparseVar2` sparse variant stores for compact, range-queryable on-disk representations of genotype data.
- `Reference` for reading reference genome sequence.
- Mutation catalogues and signature refitting via `cosmic_signatures` and `fit_signatures`.
- A `genoray index|write|view` CLI for building indices, converting VCF/PGEN to sparse formats, and inspecting variant files.

See the [`genoray-api` skill](skills/genoray-api/SKILL.md) and the [docs](docs/source/api.md) for details.

# Examples

## VCF

We work with VCFs using the (you guessed it) VCF class:

```python
from genoray import VCF

vcf = VCF("file.vcf.gz")
```

Querying data for a region is as simple as:

```python
# shape: (samples ploidy variants)
genos = vcf.read("1")  # read all variants on chromosome 1
```

You can also change the return type to be either genotypes and/or dosages by providing a `mode` argument:

```python
vcf = VCF("file.vcf.gz", dosage_field="DS")  # need a dosage_field to read dosages

genos, dosages = vcf.read("1", mode=VCF.Genos16Dosages)
```

Dosages have shape `(samples, variants)` and dtype `np.float32`.

> [!NOTE]
> VCFs must also be provided a FORMAT `dosage_field` to read dosages and this field must have `Number=A` in the header, meaning there is one value for each ALT allele. 

A key feature of `genoray` is letting you work with data that is too large to fit into memory. For example:

```python
vcf = VCF("file.vcf.gz", phasing=True)  # include phasing status

# max_mem defaults to "4g", can also be capitalized or be "GB", for example
# Genos8 reduces precision to int8 from the default int16 that cyvcf2 uses
genos = vcf.chunk("1", max_mem="4g", mode=VCF.Genos8)

for chunk in genos:
    # do something with chunk, each chunk is a NumPy array of shape (samples, ploidy+1, variants)
    ...
```
The `chunk` method will automatically chunk the data along the `variants` axis to respect the memory limit, returning a generator of data instead of everything at once.

> [!NOTE]
> We also set `phasing=True` and changed the mode to `VCF.Genos8` to read phased genotypes as `int8`. The `phasing` argument lets us have access to the phase of the genotype in the format the cyvcf2 adheres to: the 3rd entry along the ploidy axis is the phase: 0 for unphased, 1 for phased. Reducing precision to `int8` instead of `int16` reduces the memory per variant by half -- we would only need higher precision if we expected to have more than 128 alleles at a variant site.

## PGEN

```python
from genoray import PGEN

pgen = PGEN("file.pgen")
```

> [!IMPORTANT]
> PGEN files are automatically indexed on construction, creating a `<prefix>.gvi` file. This is a one-time cost to enable fast range queries, but it takes longer for larger files. Don't delete this index file unless you want to re-index the PGEN file.

We can query data for a region in the same way as VCF:

```python
# shape: (samples ploidy variants)
genos = pgen.read("1")  # read all variants on chromosome 1
genos = pgen.chunk("1")  # read all variants on chromosome 1
```

However, PGEN files also support reading multiple ranges at once since this improves throughput substantially:

```python
# shape: (samples, ploidy, variants), shape: (n_ranges+1)
genos, offsets = pgen.read_ranges('1', starts=[1, 1000, 2000], ends=[1000, 2000, 3000])
first_range_genos = genos[..., offsets[0]:offsets[1]]

genos = pgen.chunk_ranges('1', starts=[1, 1000, 2000], ends=[1000, 2000, 3000])
for range_ in genos:
    if range_ is None:
        # no data for this range
        continue
    for chunk in range_:
        # do something with chunk, each chunk is a NumPy array of shape (samples, ploidy, variants)
        ...
```

The `read_ranges` method takes starts and ends and returns data for each range and the offsets to slice out the variants for each range. Since the data is allocated as a single array, the offsets let you slice out the data for each range from the `variants` axis.

> [!NOTE]
> We do not provide an API for multi-range queries of VCFs because benchmarking showed that this provided no benefit to throughput.

Like VCF, methods for PGENs accept a `mode` argument to change the return type to include genotypes, phasing, and/or dosages:

```python
genos, phasing, dosages = pgen.read("1", mode=PGEN.GenosPhasingDosages)
```

The PGEN reader adheres to pgenlib's API, so the phasing information is in a separate boolean array instead of using an extra column like VCF/cyvcf2. The phasing information is a boolean array of shape `(samples, variants)` where `True` indicates that the genotype is phased and `False` indicates that it is unphased.

> [!IMPORTANT]
> PGEN files either store hardcalls (genotypes) or dosages, not both, and dosage PGENs infer hardcalls based on a [hardcall threshold](https://www.cog-genomics.org/plink/2.0/input#dosage_import_settings). Thus, if you want to read hardcalls that do not correspond to inferred hardcalls from a dosage PGEN, you can provide two different PGEN files to the constructor. This will read hardcalls from `hardcalls.pgen` and dosages from `dosage.pgen`. The two PGEN files must have the same samples and variants in the same order. The `dosage_path` argument is optional, and if not provided, both hardcalls and dosages will be sourced from the path argument (`"hardcalls.pgen"` in the example):

```python
pgen = PGEN("hardcalls.pgen", dosage_path="dosage.pgen", ...)
```

## Filtering

You can filter variants from VCF or PGEN files by a providing a function or [polars expression](https://docs.pola.rs/user-guide/concepts/expressions-and-contexts/) to the constructor, respectively.

For VCFs, the function must accept a [cyvcf2.Variant](https://brentp.github.io/cyvcf2/docstrings.html#cyvcf2.cyvcf2.Variant) and return a boolean indicating whether to include the site.

```python
# only include variants that are common in EUR
vcf = VCF("file.vcf.gz", filter=lambda v: v.INFO['AF_EUR'] > 0.05)
```

For PGENs, the expression operates on the `.gvi` index — a polars DataFrame with columns:
- `CHROM` — contig name
- `POS` — 1-based position
- `REF` — reference allele
- `ALT` — list of alternate alleles
- `ILEN` — list of indel lengths (one per ALT: `len(ALT) - len(REF)`, or a signed size for symbolic SVs; `null` for un-sizable symbolic/breakend alleles)

Prefer the ready-made expressions in `genoray.exprs` — `is_snp`, `is_indel`, `is_biallelic`, `is_symbolic`, `is_breakend`, `is_imprecise`, and `ILEN` — and combine them with polars operators. For custom predicates, use `pl.col("CHROM"/"POS"/"REF"/"ALT"/"ILEN")` directly.

```python
import genoray
from genoray import PGEN

# only include SNPs
pgen = PGEN("file.pgen", filter=genoray.exprs.is_snp)

# exclude symbolic alleles and breakends
pgen = PGEN("file.pgen", filter=~genoray.exprs.is_symbolic & ~genoray.exprs.is_breakend)
```

# ⚠️ Important ⚠️

- For the time being, ploidy is 2 for all classes in `genoray`, but this could be more flexible for VCFs in the future. The PGEN format does not support ploidy other than 2.
- Different file formats may use different data types for their respective representations of genotypes, phasing, and dosages.
- Ranges are 0-based, so starts begin at 0 and ends are exclusive.
- Missing genotypes and dosages are encoded as -1 and `np.nan`, respectively.
- Dosages from PGEN files may not exactly match VCF files (up to a fraction of a percent) because PLINK 2.0 must encode dosages with fixed precision which can not match what can be represented by text in a VCF (may also disagree with how BCF encodes dosage).

# Contributing

To contribute to `genoray`, please fork the repository and create a pull request. We welcome contributions of all kinds, including bug fixes, new features, and documentation improvements. Please make sure to run the tests before submitting a pull request. We provide a Pixi environment that includes all development dependencies. To use the environment, install Pixi and run `pixi run prek-install` to activate pre-commit in your clone of the repo, and then run `pixi s` in the repository root directory. `pixi s` will activate the development environment and install all dependencies. You can then run the tests using `pytest`. ❗Note that all commits must adhere to [conventional commits](https://www.conventionalcommits.org/). If you have any questions or suggestions, please open an issue on the repository.

