RealSFS

* The program was called emOptim2 in earlier versions, this has now been changed to the more appropiate: 'realSFS'

This program will estimate the SFS based on a .saf file generated from the ./angsd [options] -doSaf .

Or use it to estimate a 2dsfs by supplying 2 .saf files (also generated from ./angsd [options] -doSaf).

See also SFS Estimation and 2d SFS Estimation.

Brief overview

#1d sfs
./emOptim2 afile.saf nChr [-start FNAME -P nThreads -tole tole -maxIter  -nSites -use-BFGS ]
#2dsfs
./emOptim2 pop1.saf pop2.saf nChr1 nChr2 [-start FNAME -P nThreads -tole tole -maxIter  -nSites]

The saf files are generated using ./angsd -doSaf.

nChr,nChr1, nChr2 are the number of chromosomes in the associated saf file;

-start is a file containing a log scaled estimate of the SFS that can be used as the start point for the EM optimisation.

-tole When the difference in successive likelihood values in the EM algorithm gets below this value the optimisation will stop

-P number of threads to allocate to program

-nSites Limit the optimisation to a region of this size. If nothing is supplied the program will use the entire saf file

-maxIter maximum number of iterations in the EM algorithm

-use-BFGS use bfgs algorithm instead of the EM

Program defaults to use the EM algorithm for the optimisation. See examples below for using the bfgs optimisation.

1d SFS

emOptim2 sfstest.saf 20 -P 4 >sfs.em
emOptim2 fstest.saf 20 -P 4 -use-BFGS 1 >sfs.bfgs

The emOptim2 program will read in a block of the genome (from the .saf) file, and for this region it will estimate the SFS.

The size of the block can be choosen using -nSites argument, otherwise it will try to read in the entire saf file.

If you have .saf file larger than -nSites (you can check the number of sites in the .saf.pos file), then the program will loop over the genome and output the results for each block. So each line in your Whit.saf.ml, is an SFS for a region.

2dsfs

./emOptim2 2dsfs pop1.saf pop2.saf nChr_pop1 nChr_pop2 [-start FNAME -P nThreads -tole tole -maxIter  -nSites  ]

nChr_pop1 and nChr_pop2 are the number of chromosomes in population1 and population2 respectively see 2d SFS Estimation for full example.

Output

Main results are printed to the stdout. Values are in logspace.

NB

Use as many sites as possible, for more reliable estimates.