二核苷酸计数和频率

Question

我试图从文本文件中的序列中找到dinuc计数和频率，但我的代码只输出单核苷酸计数。在

e = "ecoli.txt"

ecnt = {}

with open(e) as seq:
    for line in seq:
        for word in line.split():
            for i in range(len(seqr)):
                dinuc = (seqr[i] + seqr[i:i+2])
                for dinuc in seqr:
                    if dinuc in ecnt:
                        ecnt[dinuc] += 1
                    else:
                        ecnt[dinuc] = 1  

for x,y in ecnt.items():
    print(x, y)

示例输入：“aaatttcgtcgtgccc”

样本输出： AA:2个 电话：3 温度：2 重心：2 GT:2台 气相色谱：1 抄送：2

现在，我只得到一个核苷酸作为我的输出：

55083600摄氏度 甲60342100 电话：88192300 国92834000

对于重复的核苷酸，即“AAA”，计数必须返回连续“AA”的所有可能组合，因此输出应该是2而不是1。不管二核苷酸以什么顺序排列，我只需要所有的组合，并让代码返回重复核苷酸的正确计数。我问助教，她说我唯一的问题是让“for”循环把二核苷酸加到字典里，我想我的范围可能错了，也可能没有错。这个文件非常大，所以序列被分成几行。在

提前谢谢你！！！在

Answer 1

使用defaultdict的绝佳机会：

from collections import defaultdict

file_name = "ecoli.txt"

dinucleotide_counts = defaultdict(int)

sequence = ""

with open(file_name) as file:
    for line in file:
        sequence += line.strip()

for i in range(len(sequence) - 1):
    dinucleotide_counts[sequence[i:i + 2]] += 1

for key, value in sorted(dinucleotide_counts.items()):
    print(key, value)

Answer 2

我查看了您的代码，发现了一些您可能需要查看的内容。在

为了测试我的解决方案，因为我没有生态.txt，我用随机核苷酸生成了自己的一个，其功能如下：

import random
def write_random_sequence():
    out_file = open("ecoli.txt", "w")
    num_nts = 500
    nts_per_line = 80
    nts = []
    for i in range(num_nts):
        nt = random.choice(["A", "T", "C", "G"])
        nts.append(nt)
    lines = [nts[i:i+nts_per_line] for i in range(0, len(nts), nts_per_line)]
    for line in lines:
        out_file.write("".join(line) + "\n")
    out_file.close()
write_random_sequence()

请注意，这个文件有一个由500个核苷酸组成的序列，每个序列分成80个核苷酸行。为了计算第一个核苷酸在一个行的末尾，第二个核苷酸在下一行的开头，我们需要将所有这些单独的行合并成一个字符串，没有空格。我们先这么做：

^{pr2}$

试着打印出“seq”，注意它应该是一个包含所有核苷酸的巨大字符串。接下来，我们需要找到序列串中的二核苷酸。我们可以用切片来做，我看你试过了。因此，对于字符串中的每个位置，我们同时查看当前的核苷酸和它后面的核苷酸。在

for i in range(len(seq)-1):#note the -1
    dinuc = seq[i:i+2]

然后我们可以计算核苷酸，并把它们存储在字典“ecnt”中，就像你以前那样。最终代码如下：

ecnt = {}
seq = ""
with open("ecoli.txt", "r") as seq_data:
    for line in seq_data:
        seq += line.strip()
for i in range(len(seq)-1):
    dinuc = seq[i:i+2]
    if dinuc in ecnt:
        ecnt[dinuc] += 1
    else:
        ecnt[dinuc] = 1
print ecnt