二核苷酸计数和频率

我试图从文本文件中的序列中找到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)

样品输入： “AAATTTCGTCGTTGCCC”

示例输出： AA：2 TT：3 TC：2 CG：2 GT：2 GC：1 CC：2

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

C 83550600 A 60342100 牛逼88192300 摹92834000

对于重复即“AAA”的核苷酸，计数必须返回的连续的“AA”所有可能的组合，所以输出应该是2，而不是1。它不事关什么样的顺序列出了二核苷酸，我只需要所有组合，并且让代码返回重复核苷酸的正确计数。我问我的助教，她说我唯一的问题是让我的'for'循环将二核苷酸添加到我的字典中，并且我认为我的范围可能是错误的也可能不错。该文件是一个非常大的文件，所以序列被分成几行。

非常感谢你提前！

来源

2016-09-15 Sabs hogan

显示样品输入的短节和相应的期望的输出。 – John1024

什么是seqr？它没有在你发布的代码段中定义 –

你的代码在很多方面都被破坏了。什么是'seqr'。为什么你在这里用空格分隔行'for line.split（）：'，是不是它应该是DNA序列呢？你不会删除换行符号。 –

我看了一下你的代码，发现了一些你可能想要看的东西。

为了测试我的解决方案，因为我没有ecoli.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个核苷酸的行。为了计算在第一行末尾有第一个核苷酸和下一行开头第二个核苷酸的二核苷酸，我们需要将所有这些单独的行合并成一个单独的字符串，而不是空格。让我们先做：

seq = "" 
with open("ecoli.txt", "r") as seq_data: 
    for line in seq_data: 
     seq += line.strip()

试着打印出“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

来源

2016-09-15 22:17:32 CodeSurgeon

使用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)

来源

2016-09-16 06:44:59 cdlane

二核苷酸计数和频率

回答

相关问题