grapp

grapp provides statistical and population genetics analyses that natively use the GRG format. GRG is an order of magnitude smaller than other representations, and can be multiple orders of magnitude faster for calculations, such as PCA.

Statgen features: GWAS, PCA, phenotype simulation (via grg_pheno_sim), linear mixed models, filtering samples and variants.

Popgen features: Polarization, integration with scikit-allel for various SFS and population statistics, conversions between tskit.TreeSequence files and various other formats.

grapp provides a command-line interface (ala plink or GCTA), and also a Python API that integrates with numpy, scipy, and pandas. Of particular interest to method developers, grapp provides linear operators that are compatible with scipy.sparse.linalg and other compatible libraries, which enables iterative numerical algorithms to scale to datasets with a million samples and hundreds of millions of variants. Some notes on usage are below, and you can also browse the Python documentation for grapp.

Check out the command-line cheatsheet for examples on how to perform many tasks. There are also some tutorials in the GRG documentation (and correspond Jupyter notebooks).

See the GRG paper and the main GRG documentation for more details about Genotype Representation Graphs.

Installation

pip install grapp

Functionality

GWAS

Perform association tests between phenotypes and genotypes.

• Command: grapp assoc [-h] [-p PHENOTYPES] [-b] [-c COVARIATES] [-o OUT_FILE] [-r] [-s] [-j JOBS] grg_input [grg_input ...]
• Python: grapp.assoc.linear_assoc_covar(), grapp.assoc.linear_assoc_no_covar()

PCA

Principal Component Analysis (PCA) using iterative methods.

• Command: grapp pca [-h] [-d DIMENSIONS] [-o PCS_OUT] [--normalize] [--pro-pca] grg_input
• Python: grapp.linalg.PCs(), grapp.linalg.eigs()

BOLT-LMM

BOLT-LMM reimplemented on GRG.

• Command: grapp bolt-lmm [-h] -p PHENOTYPES [-c COVARIATES] [--covar-cols [COVAR_COLS ...]] [--q-covar-cols [Q_COVAR_COLS ...]] [--covar-max-levels COVAR_MAX_LEVELS] [--num-calib-snps NUM_CALIB_SNPS] [--mc-trials MC_TRIALS] [--cg-tol CG_TOL] [--max-iter MAX_ITER] [--seed SEED] [-j JOBS] [-o OUT_FILE] grg_input [grg_input ...]
• Python: grapp.assoc.bolt_lmm_inf()

Filtering

Filtering variants and samples, with syntax similar to that of bcftools.

• Command: grapp filter [-h] [-j JOBS] [-S INDIVIDUALS | --hap-samples HAP_SAMPLES | -P POPULATIONS] [-r RANGE] [-c MIN_AC] [-C MAX_AC] [-q MIN_AF] [-Q MAX_AF] [-v TYPES] [-A] [-m MIN_ALLELES] [-M MAX_ALLELES] [-H HWE] [--multi-ref] grg_input grg_output
• Python: grapp.util.filter.grg_save_freq(), grapp.util.filter.grg_save_individuals(), etc..

Phenotype Simulation

Simulate phenotypes, for testing GWAS and association study scenarios.

• Command: grapp pheno [-h] [-e HERITABILITY] [-n NUM_CAUSAL] [--seed SEED] [--no-normalize] [--save-effects SAVE_EFFECTS] [-o OUT_FILE] grg_input
• Python: See grg_pheno_sim

Data Conversion

GRG can be converted back to a tabular (sparse matrix) format, IGD:

• Command: grapp export [-h] [--igd IGD | --vcf VCF] [-f] [-j JOBS] [--temp-dir TEMP_DIR] [--contig CONTIG] grg_input
• Python: grapp.util.igd.export_igd()

grapp also supports converting tskit.TreeSequence to IGD, since for really large simulated TreeSequence files (e.g., from msprime) converting them to VCF can take days (or weeks):

• Command: grapp ts2igd [-h] [-f] [-j JOBS] [--temp-dir TEMP_DIR] ts_input igd_output

You can also convert TreeSequence to GRG using the core GRG library, via grg convert ts_input igd_output. In situations where converting a TreeSequence to VCF takes days, converting to IGD takes hours, and converting to GRG takes minutes.

Simple Computations

The following information can be easily displayed on the command line: GRG summary, HWE p-values, allele counts, allele frequencies, sample identifiers, mapping between samples and populations.

• Command: grapp show [-h] [-j JOBS] [-S | -P | -i | -f | -c | -H | --multi-sites] [--multi-ref] grg_input
• Python: grapp.util.simple.hwe(), grapp.util.simple.get_zygosities(), etc..

LinearOperators

• grapp.linalg.ops_scipy.SciPyXOperator: An operator that performs matrix multiplication against the genotype matrix X (NxM) or its transpose (MxN).
• grapp.linalg.ops_scipy.SciPyXTXOperator: An operator that performs matrix multiplication against the MxM product transpose(X) * X.
• grapp.linalg.ops_scipy.SciPyStdXOperator: The same as SciPyXOperator, except the genotype matrix is standardized by using the allele frequencies (standard deviation via binomial distribution).
• grapp.linalg.ops_scipy.SciPyStdXTXOperator: The same as SciPyXTXOperator, except the genotype matrix is standardized by using the allele frequencies (standard deviation via binomial distribution).

Additionally, there is a helpful method for eigen decomposition (grapp.linalg.eigs) and PCA (grapp.linalg.PCs).

scikit-allel Integration

The primary integration with scikit-allel is via the grapp.popgen.allele_counts() and grapp.popgen.pop_allele_counts() methods, which return output compatible with many scikit-allel functions, such as allel.sequence_diversity(), allel.watterson_theta(), allel.hudson_fst(), and many others.

See the grapp.popgen documentation for more details.

Nearest Neighbors

Experimental library for nearest neighbors search over GRG.

The grapp.nn module lets you search a dataset stored as a GRG for nearest neighbors in a variety of ways:

• Similarity is either between samples (haplotypes) or mutations (variants). I.e., you can ask to find samples that are similar to a given sample, or mutations that are similar to a given mutation.
• Similarity is defined as the Hamming distance between items. The Hamming distance is just the number of differences, i.e., Hamming(A, B) = |Muts(A)| + |Muts(B)| - 2*|Muts(A) intersect Muts(B)|.
• There are APIs that let you query using a sample/mutation already in the dataset (GRG) or more generally you can query using an external sample/mutation that is not in the dataset, though your options may be slightly more limited in the latter case.