DOCSEARCHbyAlgolia
Submit Docs

scikit-allel

repository·master·Indexed 43 days ago

https://github.com/cggh/scikit-allel

A Python library for the exploratory analysis of large-scale genetic variation data. It provides tools for reading and converting VCF, GFF3, and FASTA files, as well as specialized data structures like GenotypeChunkedArray and GenotypeDaskArray for handling datasets that exceed memory. The library includes modules for calculating admixture statistics (Patterson f2, f3, and D), performing PCA and randomized PCA, and analyzing Hardy-Weinberg equilibrium and Mendelian inheritance.

Tokens
10.1K
Snippets
37
Records
60
Agent score
56%

What's inside scikit-allel

  2. Analyze Mendelian inheritance with allel.stats.mendel

    master

    The allel.stats.mendel module provides tools for analyzing Mendelian inheritance patterns, specifically for identifying errors in pedigree data and performing phasing based on transmission patterns.

    Key functionalities include:

    • Error detection: Identifying Mendelian inconsistencies using allel.mendel_errors.
    • Transmission painting: Using allel.paint_transmission to track how alleles are passed down.
    • Phasing: Reconstructing haplotypes using transmission logic via allel.phase_by_transmission, allel.phase_progeny_by_transmission, or allel.phase_parents_by_transmission.

  3. Important notice regarding scikit-allel development status

    master

    scikit-allel is currently in maintenance-only mode and is no longer under active development.

    For new projects or those requiring better scaling, performance, and maintainability, the intended successor is sgkit (the Statistical Genetics Toolkit). While sgkit is the recommended path forward, note that it may not yet have full feature parity with scikit-allel.

  4. Use Dask arrays for large-scale genomic data

    master

    For datasets that exceed available memory, scikit-allel provides Dask-backed array implementations. These classes allow you to work with large genomic datasets by leveraging Dask's lazy evaluation and parallel computing capabilities. Instead of loading entire arrays into RAM, these objects represent chunks of data that are processed on demand.

    Supported Dask-backed array types include:

    • GenotypeDaskArray: For genotype data.
    • HaplotypeDaskArray: For haplotype data.
    • AlleleCountsDaskArray: For allele counts data.

  6. Install scikit-allel via pip

    master

    If you have a C compiler on your system, you can install scikit-allel using pip.

    By default, this installs scikit-allel with minimal dependencies (numpy, cython, and dask).

    To install scikit-allel with all optional third-party dependencies (including scipy, matplotlib, seaborn, pandas, scikit-learn, h5py, and zarr), use the [full] extra.

    # Minimal installation
$ pip install scikit-allel

# Full installation with all optional dependencies
$ pip install scikit-allel[full]

  7. Compare different genotype array storage classes

    master

    scikit-allel provides several array classes for storing genotype data, each with different memory and performance characteristics:

    • allel.model.GenotypeArray: A contiguous NumPy-based array. It is fast for computations but uses more memory.
    • allel.bcolz.GenotypeCArray: A compressed array using the bcolz format. It significantly reduces memory footprint.
    • allel.bcolz.GenotypeCArray (with rootdir): A compressed, persistent array that can be stored on disk.
    • GenotypeArray with a .mask: An array that includes a boolean mask to represent missing data.
    • to_sparse(format='csr'): Converts a GenotypeArray into a SciPy sparse matrix (e.g., CSR format) for memory efficiency when data is sparse.

    Use nbytes for NumPy arrays or cbytes for bcolz arrays to check memory usage.

    # Contiguous array (numpy)
g = allel.model.GenotypeArray(data, copy=False)

# Compressed array (bcolz)
gc = allel.bcolz.GenotypeCArray(data)

# Compressed persistent array (bcolz)
import tempfile
rootdir = tempfile.mkdtemp()
gcp = allel.bcolz.GenotypeCArray(data, rootdir=rootdir, mode='w')

# Arrays with mask
gm = g.copy()
gm.mask = mask

# Sparse array (scipy)
m = g.to_sparse(format='csr')

  8. Handle genomic gaps in selection statistics

    master

    Large gaps between consecutive variants can introduce artifacts in selection statistics.

    When using allel.stats.ihs, you can control how gaps are treated using:

    • max_gap: Sets the maximum allowed distance between variants.
    • clip_gap: Sets the distance to which gaps are clipped.
    • is_accessible: Provides a mask to adjust scores based on genomic accessibility, which can mitigate the impact of gaps in non-accessible regions.

  9. Compare GenotypeArray and GenotypeCArray for FST calculations

    master

    In scikit-allel, $F_{ST}$ calculations can be performed using either allel.GenotypeArray or allel.GenotypeCArray.

    • GenotypeArray is the standard representation for genotype data.
    • GenotypeCArray is a compressed representation.

    Profiling indicates that weir_cockerham_fst can be called directly with a GenotypeCArray, and the resulting components (a, b, c) are equivalent to those produced by a GenotypeArray.

    # Using GenotypeArray
g = allel.GenotypeArray(data)
a, b, c = allel.stats.weir_cockerham_fst(g, subpops)

# Using GenotypeCArray
gc = allel.GenotypeCArray(data)
ac, bc, cc = allel.stats.weir_cockerham_fst(gc, subpops)

# Results are equivalent
import numpy as np
assert np.array_equal(a, ac)

  10. Map alleles in a GenotypeArray using Dask for parallel processing

    master

    When working with large datasets that exceed memory, you can perform allele mapping in parallel using dask.array. Since GenotypeArray.map_alleles is a NumPy-based operation, you can use dask.array.map_blocks to apply the mapping to chunks of the data.

    To use this pattern, you must ensure that the mapping array is also chunked such that its first dimension aligns with the chunks of the GenotypeArray.

    Implementation Pattern:

    1. Convert the GenotypeArray to a Dask array (gd).
    2. Convert the mapping array to a Dask array (md).
    3. Define a helper function f(block, bmapping) that wraps the GenotypeArray constructor and calls .map_alleles() on the block.
    4. Use da.map_blocks to apply the function across the chunks.
    import allel
import dask.array as da
import numpy as np

# Setup GenotypeArray and mapping
g = allel.GenotypeArray(data, dtype='i1')
mapping = np.array(mapping_data, dtype='i1')

# Create Dask arrays
chunks_dim0 = 2
gd = da.from_array(g, chunks=(chunks_dim0, 2, None))
md = da.from_array(mapping, chunks=(chunks_dim0, None))

# Mapping function for blocks
def f(block, bmapping):
    return allel.GenotypeArray(block).map_alleles(bmapping[:, 0, :])

# Apply mapping in parallel
gmapped1 = da.map_blocks(f, gd, md[:, None, :], chunks=gd.chunks)

# Compute result back to GenotypeArray
actual = allel.GenotypeArray(gmapped1.compute())

  11. Parallelize computations with Dask and Blosc

    master

    When using GenotypeDaskArray with bcolz data, performance is influenced by two levels of threading:

    1. Blosc threads: Controlled via bcolz.blosc_set_nthreads(n). This manages the decompression parallelism.
    2. Dask workers: Controlled via the num_workers argument in .compute(). This manages the task-level parallelism.

    To optimize performance, you should experiment with different combinations of these two thread counts. For example, using 4 Blosc threads and 4 Dask workers can often result in a total wall time that is less than the sum of decompression and computation times, as Dask can overlap these tasks.

    import bcolz

# Set Blosc decompression threads
bcolz.blosc_set_nthreads(4)

# Execute Dask computation with a specific number of workers
# This allows Dask to parallelize the workload across 'm' workers
result = gd.max().compute(num_workers=4)