Genotype Likelihoods

Genotype likelihoods from alignments

<classdiagram> // [input|bam files;SOAP files{bg:orange}]->[sequence data]

[sequence data]->[genotype likelihoods|samtools;GATK;soapSNP;kim et.al]
</classdiagram>

-GL [int]

If your input is sequencing file you can estimate genotype likelhoods from the mapped reads. Four different methods are available.

Samtools

-GL 1

This methods has a random component. In same tools there is a stocastic component so to get the exact same results as samtools use nThreads=1. However, the method is still the same with multiple threads but some sites will have small differences compared to the samtools output bacause of the stocastic component.

options

-minQ [int]

default 13. The minimum allowed base quality score.

-minMapQ [int]

default 0; The minimum allowed mapping quality score.

example

./angsd -bam bam.filelist -GL 1 -out outfile

GATK

-GL 2

options

-minQ [int]

default 13. The minimum allowed base quality score.

-minMapQ [int]

default 0; The minimum allowed mapping quality score.

example

./angsd -bam bam.filelist -GL 2 -out outfile

soapSNP

-GL 3 When estimating GL with soapSNP we need to generate a calibration matrix. This is done automaticly if these doesn't exist. These are located in angsd_tmpdir/basenameNUM.count,angsd_tmpdir/basenameNUM.qual

options

-minQ [int]

default 13. The minimum allowed base quality score.

-minMapQ [int]

default 0; The minimum allowed mapping quality score.

-maxq [int]

default 51; The maximum allowed base quality score.

-L [int]

default 150; The maximum read length (choosing one that is too large is not a problem)

example

./angsd -bam bam.filelist -GL 3 -out outfile -minQ 0 -ref hg19.fa 

This first loop doesn't estimate anything else than the calibration matrix. So now we can do the analysis we want

./angsd -bam bam.filelist -GL 3 -out outfile -minQ 0 -ref hg19.fa

Kim et al.

-GL 4 Citation Citation

options

-error [filename]

A file with the estimated type specific error rates (see Error_estimation).

example

./angsd -bam bam.filelist -GL 4 -out outfile -error error.file 

output genotype likelihoods

-doGlf [int]

Output the log genotype likelihoods to a file

binary

Glf file in binary doubles. All 10 genotype likelihoods are printed to a file. For each printed site there are 10*N doubles where N is the number of individuals. The order of the 10 genotypes are alphabetical AA AC AG AT CC CG CT GG GT TT.

Beagle format

Beagle haplotype imputation and be performed directly on genotype likelhoods. To generate beagle input file use

-doGlf 2

In order to make this file the major and minor allele has the be inferred (-doMajorMinor). It is also a good idea to only use the polymorphic sites.

Example

In this example our input files are bam files. We use the samtools genotype likelihood methods. We use 10 threads. We infer the major and minor allele from the likelihoods and estimate the allele frequencies. We test for polymorphic sites and only outbut the ones with are likelhood ratio test statistic of minimum 24 (ca. p-value<1e-6).

./angsd -GL 1 -out genolike -nThreads 10 -doGlf 2 -doMajorMinor 1 -minLRT 24 -doMaf 2 -doSNP 2 -bam bam.filelist

output

The above command generates the file genolike.beagle.gz that can be use as input for the beagle software

marker  allele1 allele2 Ind0    Ind0    Ind0    Ind1    Ind1    Ind1    Ind2    Ind2    Ind2    Ind3    Ind3    Ind3 
1_14000023      1       0       0.941177        0.058822        0.000001        0.799685        0.199918        0.000397        0.666316        0.333155        0.000529 
1_14000072      2       3       0.709983        0.177493        0.112525        0.941178        0.058822        0.000000        0.665554        0.332774        0.001672
1_14000113      0       2       0.855993        0.106996        0.037010        0.333333        0.333333        0.333333        0.799971        0.199989        0.000040 
1_14000202      2       0       0.835380        0.104420        0.060201        0.799685        0.199918        0.000397        0.333333        0.333333        0.333333
...

Note that the above values sum to one per sites for each individuals. This is just a normalization of the genotype likelihoods in order to avoid underflow problems in the beagle software it does not mean that they are genotype probabilities.