The script takes a FASTA format file in input and outputs the frequencies of each amino acid (A
, C
, T
and G
) at each site per population per iteration. The name of each sequence looks like >it0pop1ind0locus0
. In this example, this indicates that the sequence that follows is the DNA sequence for the zeroth iteration (it
), the first population (pop
), the zeroth individual (ind
) at the zeroth locus (locus
).
My goal (just to restate it) is to get the frequencies of each amino acid (A
, C
, T
and G
) at each site per population per iteration.
from __future__ import division
import os
from Bio import SeqIO
import itertools
import re
import sys
os.chdir("/path/to/directory/")
#### ReFormat the data ####
maxit = 1
maxpop = 1
maxind = 1
maxlocus = 1
for pos, i in enumerate(SeqIO.parse("0.0.1.fasta", "fasta")):
if pos == 0:
maxsite = len(i.seq)
v1 = re.findall(r'\d+', i.id)
v = []
for val in v1:
v.append(int(val))
maxit = max(v[0]+1, maxit)
maxpop = max(v[1]+1, maxpop)
maxind = max(v[2]+1, maxind)
maxlocus = max(v[3]+1, maxlocus)
all = [[[[ "x" * maxsite for locus in xrange(maxlocus)] for ind in xrange(maxind)] for pop in xrange(maxpop)] for it in xrange(maxit)]
for i in SeqIO.parse("0.0.1.fasta", "fasta"):
v1 = re.findall(r'\d+', i.id)
v = []
for val in v1:
v.append(int(val))
it = v[0]
pop = v[1]
ind = v[2]
locus = v[3]
all[it][pop][ind][locus] = str(i.seq)
#### Calculate the frequencies ####
pA = [[[2 for it in xrange(maxsite)] for pop in xrange(maxpop)] for site in xrange(maxit)]
pC = [[[2 for it in xrange(maxsite)] for pop in xrange(maxpop)] for site in xrange(maxit)]
pG = [[[2 for it in xrange(maxsite)] for pop in xrange(maxpop)] for site in xrange(maxit)]
pT = [[[2 for it in xrange(maxsite)] for pop in xrange(maxpop)] for site in xrange(maxit)]
for it in xrange(maxit):
for pop in xrange(maxpop):
for locus in xrange(maxlocus):
for site in xrange(maxsite):
x = []
for ind in xrange(maxind):
x.append(all[it][pop][ind][locus][site])
pA[it][pop][site] = x.count("A")/maxind
pT[it][pop][site] = x.count("T")/maxind
pC[it][pop][site] = x.count("C")/maxind
pG[it][pop][site] = x.count("G")/maxind
Here are the values I get:
>>> pA [[[0.16666666666666666, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.6666666666666666, 0.0], [0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0]]] >>> pC [[[0.8333333333333334, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0], [1.0, 0.0, 0.0, 0.6666666666666666, 0.0, 0.0, 0.0, 0.0, 1.0]]] >>> pG [[[0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 1.0, 0.3333333333333333, 0.0, 1.0, 0.0, 0.0, 0.0]]] >>> pT [[[0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3333333333333333, 0.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]]]
Here is the data for the example (FASTA). In this example, we have only 1 iteration (it
), 2 populations (pop
), 5 individuals (ind
), 1 locus (locus
) and 9 sites (length of the sequence).
>it0pop0ind0locus0 CGGCAGAAC >it0pop0ind1locus0 CGGCAGATC >it0pop0ind2locus0 CGGCAGAAC >it0pop0ind3locus0 CGGCAGAAC >it0pop0ind4locus0 CGGCAGATC >it0pop0ind5locus0 AGGCAGAAC >it0pop1ind0locus0 CGGCAGAAC >it0pop1ind1locus0 CGGCAGAAC >it0pop1ind2locus0 CGGGAGAAC >it0pop1ind3locus0 CGGCAGAAC >it0pop1ind4locus0 CGGGAGAAC >it0pop1ind5locus0 CGGCAGAAC