6

I am writing a function that is supposed to go through a .fasta file of DNA sequences and create a dictionary of nucleotide (nt) and dinucleotide (dnt) frequencies for each sequence in the file. I am then storing each dictionary in a list called "frequency". This is the piece of code that is acting strange: 

for fasta in seq_file:
    freq = {}
    dna = str(fasta.seq)
    for base1 in ['A', 'T', 'G', 'C']:
        onefreq = float(dna.count(base1)) / len(dna)
        freq[base1] = onefreq
        for base2 in ['A', 'T', 'G', 'C']:
            dinucleotide = base1 + base2
            twofreq = float(dna.count(dinucleotide)) / (len(dna) - 1) 
            freq[dinucleotide] = twofreq
    frequency.append(freq)

(I am using biopython, by the way, so that I do not have to commit the entire fasta file to memory. Also this is for ssDNA, so I don't have to account for anti-sense dnt)

The frequencies that are recorded for the single nt add to 1.0, but the frequencies for the dnt do not add to 1.0. Which is od since the method of calculating the two types of frequencies are identical in my eyes.

I left the diagnostic print-statements and the "check" variables in:

for fasta in seq_file:
    freq = {}
    dna = str(fasta.seq)
    check = 0.0
    check2 = 0.0
    for base1 in ['A', 'T', 'G', 'C']:
        onefreq = float(dna.count(base1)) / len(dna)
        freq[base1] = onefreq
        check2 += onefreq
        for base2 in ['A', 'T', 'G', 'C']:
            dinucleotide = base1 + base2
            twofreq = float(dna.count(dinucleotide)) / (len(dna) - 1) 
            check += twofreq
            print(twofreq)
            freq[dinucleotide] = twofreq
    print("\n")
    print(check, check2)
    print(len(dna))
    print("\n")
    frequency.append(freq)

to get this output: (for only one sequence)

0.0894168466523 
0.0760259179266
0.0946004319654
0.0561555075594
0.0431965442765
0.0423326133909
0.0747300215983
0.0488120950324
0.0976241900648
0.0483801295896
0.0539956803456
0.0423326133909
0.0863930885529
0.0419006479482
0.0190064794816
0.031101511879


(0.9460043196544274, 1.0)
2316

Here we can see the frequency of each of the 16 different dnt possible, the sum of all dnt frequencies (0.946) and the sum of all nt frequencies (1.0) and the length of the sequence.

Why does the dnt frequency not add up to 1.0?

Thanks for your help. I'm very new to python, and this is my first question, so I hope that this submissions is acceptable.

Bantha
  • 63
  • 4
  • Are you sure you're using Python 2? It looks like you're using the `print()` function, not the `print` statement. – TigerhawkT3 May 27 '15 at 16:36
  • I'm using python 2. I'm just used to writing print() functions from a little adventure I once had in Uni with C. And since it has never given me trouble for it, I never changed to the print statement. – Bantha May 27 '15 at 16:43
  • In that case, I'd recommend using `from __future__ import print_function` so that (for example) your `print(check, check2)` prints `check` and then `check2` separated by a space instead of printing the tuple `(check, check2)`. – TigerhawkT3 May 27 '15 at 16:56
  • your problem ..... try with "AAAAA" ;-) – Jose Ricardo Bustos M. May 27 '15 at 17:00

3 Answers3

3

your problem, try with following fasta:

>test
AAAAAA

"AAAAAA".count("AA")

you get:

3

It should be

5

reason

from documentation: count return the number of (non-overlapping) occurrences of substring sub in string s[start:end]

solution using Counter and chunk function

from Bio import SeqIO

def chunks(l, n):
  for i in xrange(0, len(l)-(n-1)):
    yield l[i:i+n]

from collections import Counter

frequency = []
input_file = "test.fasta"
for fasta in SeqIO.parse(open(input_file), "fasta"):
  dna = str(fasta.seq)
  freq = Counter(dna)   #get counter of single bases
  freq.update(Counter(chunks(dna,2))) #update with counter of dinucleotide
  frequency.append(freq)

for "AAAAAA" you get:

Counter({'A': 6, 'AA': 5})
Jose Ricardo Bustos M.
  • 8,016
  • 6
  • 40
  • 62
2

str.count() doesnt count overlapping motif it find.

Exemple:

If you have 'AAAA' in your sequence and you look for the dinucleotide 'AA', you expect than 'AAAA'.count('AA') return you 3, but it will return 2. So:

print float('AAAA'.count('AA')) / (len('AAAA') - 1)
0.666666

instead of 1

You can just change the line where you count the frequency by:

twofreq = len([i for i in range(len(dna)-1) if dna[i:i+2] == dinucleotide]) / float((len(dna) - 1))
Data_addict
  • 300
  • 1
  • 9
2

You're scanning the string far more than you need to - 20 times, in fact. This may not matter for small test sequences, but it will be noticeable as they get larger. I would recommend a different approach, which resolves the issue with the overlap as a side effect:

nucleotides = [ 'A', 'T', 'G', 'C' ]
dinucleotides = [ x+y for x in nucleotides for y in nucleotides ]
counts = { x : 0 for x in nucleotides + dinucleotides }

# count the first nucleotide, which has no previous one
n_nucl = 1
prevn = dna[0]
counts[prevn] += 1

# count the rest, along with the pairs made with each previous one
for nucl in dna[1:]:
    counts[nucl] += 1
    counts[prevn + nucl] += 1
    n_nucl += 1
    prevn = nucl

total = 0.0
for nucl in nucleotides:
    pct = counts[nucl] / float(n_nucl)
    total += pct
    print "{} : {} {}%".format(nucl, counts[nucl], pct)
print "Total : {}%".format(total) 

total = 0.0
for dnucl in dinucleotides:
    pct = counts[dnucl] / float(n_nucl - 1)
    total += pct
    print "{} : {} {}%".format(dnucl, counts[dnucl], pct)
print "Total : {}%".format(total)

This approach only scans through the string once, although it is admittedly more code...

twalberg
  • 59,951
  • 11
  • 89
  • 84