Dosage

bgen_reader.compute_dosage(expec, alt=None)

Compute dosage from allele expectation.

Parameters:

• expec (array_like) – Allele expectations encoded as a samples-by-alleles matrix.

• alt (array_like, optional) – Alternative allele index. If None, the allele having the minor allele frequency for the provided expec is used as the alternative. Defaults to None.

Returns:

Dosage encoded as an array of size equal to the number of samples.

Return type:

numpy.ndarray

Examples

First a quick-start example.

>>> from bgen_reader import allele_expectation, compute_dosage
>>> from bgen_reader import example_filepath, read_bgen
>>>
>>> filepath = example_filepath("example.32bits.bgen")
>>>
>>> # Read the example.
>>> bgen = read_bgen(filepath, verbose=False)
>>>
>>> # Extract the allele expectations of the fourth variant.
>>> variant_idx = 3
>>> e = allele_expectation(bgen, variant_idx)
>>>
>>> # Compute the dosage when considering the first allele
>>> # as the reference/alternative one.
>>> alt_allele_idx = 1
>>> d = compute_dosage(e, alt=alt_allele_idx)
>>>
>>> # Print the dosage of the first five samples only.
>>> print(d[:5])
[1.96185308 0.00982666 0.01745552 1.00347899 1.01153563]

Genotype probabilities, allele expectations and frequencies.

>>> from bgen_reader import (
...     allele_expectation,
...     allele_frequency,
...     compute_dosage,
...     example_filepath,
...     read_bgen,
... )
>>> from pandas import DataFrame
>>> from xarray import DataArray
>>>
>>> # Download an example
>>> filepath = example_filepath("example.32bits.bgen")
>>>
>>> # Open the bgen file.
>>> bgen = read_bgen(filepath, verbose=False)
>>> variants = bgen["variants"]
>>> genotype = bgen["genotype"]
>>> samples = bgen["samples"]
>>>
>>> variant_idx = 3
>>> variant = variants.loc[variant_idx].compute()
>>> # Print the metadata of the fourth variant.
>>> print(variant)
        id    rsid chrom   pos  nalleles allele_ids  vaddr
3  SNPID_5  RSID_5    01  5000         2        A,G  16034

>>> geno = bgen["genotype"][variant_idx].compute()
>>> metageno = DataFrame({k: geno[k] for k in ["ploidy", "missing"]},
...                      index=samples)
>>> metageno.index.name = "sample"
>>> print(metageno)  
            ploidy  missing
sample
sample_001       2    False
sample_002       2    False
sample_003       2    False
sample_004       2    False
...            ...      ...
sample_497       2    False
sample_498       2    False
sample_499       2    False
sample_500       2    False

[500 rows x 2 columns]
>>> p = DataArray(
...     geno["probs"],
...     name="probability",
...     coords={"sample": samples},
...     dims=["sample", "genotype"],
... )
>>> # Print the genotype probabilities.
>>> print(p.to_series().unstack(level=-1))  
genotype          0        1        2
sample
sample_001  0.00488  0.02838  0.96674
sample_002  0.99045  0.00928  0.00027
sample_003  0.98932  0.00391  0.00677
sample_004  0.00662  0.98328  0.01010
...             ...      ...      ...
sample_497  0.00137  0.01312  0.98550
sample_498  0.00552  0.99423  0.00024
sample_499  0.01266  0.01154  0.97580
sample_500  0.00021  0.98431  0.01547

[500 rows x 3 columns]
>>> alleles = variant["allele_ids"].values[0].split(",")
>>> e = DataArray(
...     allele_expectation(bgen, variant_idx),
...     name="expectation",
...     coords={"sample": samples, "allele": alleles},
...     dims=["sample", "allele"],
... )
>>> # Print the allele expectations.
>>> print(e.to_series().unstack(level=-1))  
allele            A        G
sample
sample_001  0.03815  1.96185
sample_002  1.99017  0.00983
sample_003  1.98254  0.01746
sample_004  0.99652  1.00348
...             ...      ...
sample_497  0.01587  1.98413
sample_498  1.00528  0.99472
sample_499  0.03687  1.96313
sample_500  0.98474  1.01526

[500 rows x 2 columns]
>>> rsid = variant["rsid"].values[0]
>>> chrom = variant["chrom"].values[0]
>>> variant_name = f"{chrom}:{rsid}"
>>> f = DataFrame(allele_frequency(e), columns=[variant_name], index=alleles)
>>> f.index.name = "allele"
>>> # Allele frequencies.
>>> print(f) 
        01:RSID_5
allele
A       305.97218
G       194.02782
>>> alt = f.idxmin().values[0]
>>> alt_idx = alleles.index(alt)
>>> d = compute_dosage(e, alt=alt_idx).to_series()
>>> d = DataFrame(d.values, columns=[f"alt={alt}"], index=d.index)
>>> # Dosages when considering G as the alternative allele.
>>> print(d) 
              alt=G
sample
sample_001  1.96185
sample_002  0.00983
sample_003  0.01746
sample_004  1.00348
...             ...
sample_497  1.98413
sample_498  0.99472
sample_499  1.96313
sample_500  1.01526

[500 rows x 1 columns]