二核苷酸計數和頻率

我試圖從文本文件中的序列中找到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

二核苷酸計數和頻率

回答

相關問題