8.2 sfsStats

Estimating mutation load across the genome

mutationLoad estimates the proportion of sites in the genome which are homozygous for the preferred allele (Pi_rr), heterozygous for the preferred and any derived allele (Pi_ra), homozygous for any derived allele, with both derived alleles being the same base (Pi_aa) or heterozygous for two different derived alleles (Pi_ab). Recessive load is defined as Pi_aa + Pi_ab and dominant load is defined as Pi_ra + Pi_aa + Pi_ab.

The task accepts either a glf or a bam file as input and requires a file defining the preferred allele at every position where mutation load is estimated from. This file can take on two different formats. i) It can be in the format of an alleles file with custom preferred allele at every position, usually generated by assuming that the ancestral allele of outgroup data (the ancestral state) is the preferred allele, e.g. using the task ancestralAlleles. ii) It can be in the format of a bed file. In that case, the task requires the reference genome to be provided as well and will automatically take the reference base at each position as the preferred allele.

8.2.1 Input

Required inputs :

--bam Input_bam_file.bam	Input bam file (can be replaced with a glf file)
--glf Input_glf_file.glf.gz	Input glf file (can be replaced with a bam file)
--alleles Input_ancestral_alleles.alleles	Input alleles file defining the position to be used and the preferred allele at every position (can be replaced with a bed + reference genome)
--regions Input_regions.bed	Input bed file defining the positions to be used. Using this option also requires a reference genome. (can be replaced with an alleles file)
--fasta Input_refrence_genome_file.fasta	Reference genome. Required if a bed file is provided (can be replaced with an alleles file)

Example alleles file:

Chr Pos Allele1
chr1 1000 A
chr1 1001 T
chr2 3000 C

See Filter parameters to apply specific filters for bases, reads and parsing window setting. Engine parameters that are common to all tasks can be found here.

8.2.2 Output

*_mutationLoad.txt

Genomewide mutation load output file showing proportions of sites in mutation load as described above.

8.2.3 Usage Example

#! /bin/bash

# Set atlas path
atlas=$(dirname "$0")/../build/atlas

# Simulate a BAM File
$atlas simulate --chrLength 100000

# Create a bed file with every 7th position
printf "chr1\t0\t1\n" > bed.bed
for i in {7..10000..7}; do
    printf "chr1\t$i\t$((i+1))\n" >> bed.bed
done

# Calculate mutation load
$atlas mutationLoad --bam ATLAS_simulations.bam --fasta ATLAS_simulations.fasta --regions bed.bed