python - 频率加起来不为一

Question

我正在编写一个函数，它应该通过 DNA 序列的 .fasta 文件并为文件中的每个序列创建核苷酸 (nt) 和二核苷酸 (dnt) 频率的字典。然后我将每个字典存储在一个名为“频率”的列表中。这是一段奇怪的代码:

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)

(顺便说一句，我正在使用 biopython，这样我就不必将整个 fasta 文件提交到内存中。这也是针对 ssDNA，所以我不必考虑反义 dnt)

记录的单个 nt 的频率加起来为 1.0，但 dnt 的频率加起来不为 1.0。这很奇怪，因为在我看来，计算这两种频率的方法是相同的。

我将诊断打印语句和“检查”变量留在了:

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)

得到这个输出:(只针对一个序列)

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

在这里我们可以看到 16 种不同 dnt 可能的频率、所有 dnt 频率的总和 (0.946) 和所有 nt 频率的总和 (1.0) 以及序列的长度。

为什么dnt频率加起来不等于1.0？

感谢您的帮助。我是 python 的新手，这是我的第一个问题，所以我希望这个提交是可以接受的。

Answer 1

您的问题，请尝试以下 fasta:

>testAAAAAA

"AAAAAA".count("AA")

你得到:

3

应该是

5

原因

来自文档:count 返回子字符串 sub 在字符串 s[start:end] 中(非重叠)出现的次数

解决方案 使用 Counter和 block 函数

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)

对于“AAAAAA”你得到:

Counter({'A': 6, 'AA': 5})