# Efficient algorithm to calculate various population divergence statistics

Intro and description of the data

I am simulating the evolution of very long DNA sequences. The model works well, is performant and will output data in the following kind of fasta format

>it0pop0ind0locus0
ATGTTG...
>it0pop0ind1locus0
ATGTTG...
>it0pop0ind2locus0
ATGTTG...
>it0pop0ind3locus0


, where it stands for iteration (replicate of the same simulation), pop stands for population (or subpopulation of the metapopulation if you prefer), ind stands for the individual chromosome that was sampled, locus stands for locus which is defined as a very big sequence (0.1GB or 1 GB maybe).

Issues

I never had to analyze genetic data (or only during my Bachelor degree). As I know nothing about the available algorithms that exist to make this kind of analysis, I first thought I would just make my own code (in Python). It turns out that I might have a lot of data to analyse and my Python code will be way too slow. I may also experience RAM-related issues

Moreover, I want to measure a whole series of different statistics of population divergence (Fst for each site, Fst averaged over many sites following Cochramm and Weir, Gsd, absolute number of fixed sites that differ, etc..), so the algorithm would need to be quite flexible.

Question

There are a bunch of existing efficient algorithms that would eventually fit my needs. Can you please give me some recommendations and make a few comparisons of what's available out there?

An efficient algorithm that takes in input:

• a description of the naming convention in the fasta file, or a list of positions for each locus and each population in the fasta file.

• A description of how I want the statistics to be calculated

and that outputs those statistics. That'd be great!

• Have you looked at BioPython? Specifically the Genepop module (biopython.org/wiki/PopGen_Genepop) looks useful. If you've already got some Python experience, it might be the fastest way to get started. – kmm Jul 3 '15 at 1:30
• Are these sequences of equal length? – WYSIWYG Jul 3 '15 at 6:39
• I would suggest that instead of reading and outputting fasta files, you simply store the co-ordinates and the mutation locations as a numeric matrix. Use unsigned integers for storing co-ordinates. Use 32 bit (or if possible 16bit) float for storing other kinds of data such as frequency. If <1% variation in frequency won't matter to you then you can use unsigned short integers for saving them. Or if your sensitivity is 0.1% then multiply save them as per-thousand. In any case, even using regular floats should not impose a big load on memory. – WYSIWYG Jul 3 '15 at 6:53
• @WYSIWYG Yes, the solution was pretty obvious, yet I didn't think of that at first. Do you want to make an answer out of your comment? – Remi.b Jul 5 '15 at 15:35