Diversity reduction due to linked deleterious mutation

This is bgshr, a package for computing the expected reduction in diversity due to linked deleterious mutations.

Installation

Dependencies are listed in requirements.txt. bgshr can be installed directly from github using pip. In this sample we create a Python virtual environment and install the package:

python -m venv ~/path/to/venv
source ~/path/to/venv/bin/activate
pip install git+https://github.com/apragsdale/bgshr.git

Alternatively, we can clone and build the package locally:

git clone https://github.com/apragsdale/bgshr.git

Then we cd into the cloned directory and perform a local installation.

cd bgshr
pip install .

Dependencies

If bgshr was cloned locally, we can install dependencies listed in requirements.txt using:

pip install -r requirements.txt

Minimal dependencies are the ubiquitous numpy, scipy and pandas. It is highly useful to also have ipython installed. Further, jupyter is necessary to run the example notebooks in examples/.

Usage

To calculate diversity reduction (B values) at a site, we use a lookup table of two-locus predictions for relative reduction in diversity, and combine the effects of many loci multiplicatively. A first-order correction for interference can be applied by adjusting parameters based on the local rescaling of effective population sizes. We generate lookup tables using the moments++ model available here: https://github.com/gvbarroso/momentspp/tree/main.

In addition to the lookup table, we need one or more arrays of constrained genomic elements that experience selected mutations, an array indicating estimates of the per-base deleterious mutation rate, and a recombination map or per-base recombination rate. Along with these data we require parameters for the relevant gamma-distributed DFEs, which may include a point mass of neutral mutations.

Example usage is shown in human_chr22_cds_regulatory_model.ipynb.

Current features include:

• Discretized gamma DFEs
• Handling of mutation maps
• Allows multiple element types (with their own DFEs)
• Likelihood function for inferring drift-effective Ne

Command-line interface

Predicting B-values with a known demographic model

You can predict B-values and diversity from the command line using bgshr predict_B. This tool takes model predictions stored in a lookup table and several pieces of genomic data described below, and produces a table of predictions averaged in regular genomic intervals (saved as .csv). By default, fields saved to output are:

• chrom: Chromosome number (obtained from the mask file).
• chromStart, chromEnd: Starts and ends of windows.
• num_sites: Number of accessible sites, as specified by a genetic mask.
• exp_pi: Expected diversity at accessible sites, averaged in windows.
• B: Predicted diversity reduction due to BGS at the centers of windows (but see --resolution for a complication).

Additional fields can be added by giving the flag --rich:

• exp_del_pi: Expected diversity at constrained sites (equals 0 in windows with no constrained sites).
• avg_mut: Average mutation rate.
• del_mut: Average mutation rate at constrained sites.
• avg_rec: Average recombination rate across the interval of the window.

Required arguments for this command are:

-t, --lookup_tbl: A lookup table generated by moments++, holding two-locus model predictions, saved as a .csv file (may be gzipped). See the documentation at https://github.com/gvbarroso/momentspp/tree/main for more details.

-b, --bed: One or more BED files recording intervals where sites are selectively constrained. If more than one file is given, DFE parameters for each set of intervals should be given in the same order. If a site is present in more than one BED file, order determines priority and the site is assigned to the first constraint class where it occurs.

--shapes: Shape parameters of the DFE(s), given in the same order as BED files. All DFEs are gamma or gamma-neutral distributions; gamma-neutral distributions place a fraction of probability mass (p_neu) at selective neutrality.

--scales: Scale parameters of the DFE(s), in physical units (s).

--out: Path to output .csv file. The output file is a table in the BEDGRAPH format, with predicted B, expected diversity averaged across windows, and average recombination and mutation rates.

Additional optional arguments:

--help: Prints descriptions of each argument.

--p_neus: Neutral fraction parameters for each DFE. To use a mix of gamma and gamma-neutral DFEs, specify 0 for the pure gamma distributions, e.g., --p_neus 0.3 0 0 determines that the DFE of the first class of elements is gamma-neutral, while the second and third are pure gamma distributions.

--rmap: Path to a recombination map file. This is expected to be tab- separated text file (.tsv, .txt) with columns specifying the recombination rate on irregular genomic intervals. The file may be gzipped.

--rmap_pos_col: Specifies the column in --rmap from which positions that define the starts and ends of genomic intervals are loaded. By default, this is Position(bp).

--rmap_rate_col: The column in --rmap from which recombination rates should be loaded (by default, Rate(cM/Mb)). It is expected that rates are given in centiMorgans/Mb.

-r, --rec_rate: Optional uniform recombination rate, which can be used in place of an estimated map. If used, you must also give the length of the chromosome to be modeled with -L.

-L, --L: The length of the chromosome; required if a uniform recombination (-r) or mutation (-u) map is to be used.

--umap": Path to a file specifying the mutation map. This can be a binary .npy file holding rates at each base pair (in which missing data should be nan, or a comma or tab-separated text file (.csv, .tsv), which may be gzipped. A mutation map text file should have this format;

chrom	start	end	rate
22	0	1000	nan

where the starts and ends of intervals must be given in the (1-indexed) second and third columns (the names of the columns do not matter).

--umap_rate_col: Specifies the column from which mutation rates should be loaded in --umap is a text file (defaults to rate).

-u, --mut_rate: Optional uniform mutation rate.

-c, --n_corrs: Number of interference corrections to run.

--save_corrs: If given, the B landscape at each interference correction is retained and saved in the output file, with column names B_0, B_1, ..., B.

--spacing: Distance between focal sites (sites where B is computed), in base pairs (default 1000).

--window_size: Window size to use for aggregating constrained elements together (default 1000).

--resolution: Window size of output table. By default, this will equal --spacing, but a larger or smaller value may be used if desired. If --resolution does not equal --spacing, cubic interpolation is used to calculate B at the center of each window in the output table.

--n_cores: Number of cores across which to parallelize computation (with the default value 1, bgshr loops across blocks of focal sites).

--chunk_size: Number of focal sites for which computations should be performed at once on each processor core (default 100). Specifying many cores and a large chunk size will result in a large memory footprint.

--B_unlinked: Optional unlinked diversity reduction due to deleterious mutations on other chromosomes.

--cbgs_start: Optional minimum selection coefficient for classic BGS (CBGS) predictions. Note that we define the fitness of the heterozygote as 1+s, so s<0. Defaults to the smallest selection coefficient in the input lookup table, following any scaling.

--n_s_cbgs: Number of selection coefficients to append to the lookup table with CBGS (default 20).

-v, --verbose: If given, prints updates throughout prediction.

--rich: If given, emits extra information in the output table.

--Ne: Effective population size to use for prediction. The parameters of the input lookup table are scaled to match this parameter. If we use a lookup table computed for an equilibrium population, the scaling is exact to numerical error. For non-equilibrium tables, scaling is carried out with respect to the ancestral population size. In this case, a twofold scaling (--Ne is twice the ancestral size of the lookup table demography) can be interpreted as a doubling of all time and effective size parameters in the model.

--mask: BED mask file to use for filtering sites recorded in the output table. If no mask is given, the coverage of the provided mutation map is used to filter sites represented in the output table. Hence, if a uniform mutation map is used and no mask is provided, all sites from 1 to L will be represented in output.

Fitting the drift-effective population size

If we wish to fit the effective population size parameter for an equilibrium population, we can do so with bgshr fit_Ne, which optimizes the parameter by maximizing the likelihood of an empirical diversity dataset given Ne. Note that because prediction must be carried out several times to fit Ne, this procedure is computationally expensive and may take a long time. This tool takes all the arguments that bgshr predict_B does, with the exception that

--fit_Ne: Specifies the initial guess for optimization. Unlike with the prediction script, this argument is required for fitting.

Additional required arguments are:

--ndns: Path to a .npy binary file with base-pair-level counts of same-by-state and different-by-state allele copies (nd and ns) in its zeroth and first columns.

--max_iter: Maximum number of function calls during fitting. To achieve convergence, values near 30 are desirable.

--log_out: Path for the output log file. This file records the fitted value, log-likelihood, and warning flag output by the scipy function used for optimization. It also retains the runtime, Ne parameter and likelihood at each function call.

The fit_Ne tool outputs a table of predictions with the same data columns as does predict_B (note that the --rich argument is not presently implemented for the fitting script: the additional fields are always saved).