2014-09-18 80 views
0

業餘Python編碼器試圖在這裏學習,所以我想問我的腳本怎麼回事。我真的不知道哪裏出了問題。 (考慮第24行或difference = "%s-%s [%d]" %(seq1[i], seq2[i], i))。序列比較功能沒有按預期工作

該函數用於載入序列,重命名它們(此位正在工作),然後將每個序列與參考序列(此處僅爲文件中的第一個序列)進行比較,如果序列中的某個字母不匹配參考打印的差異和位置。但是,你可以看到下面這不工作

下面是一個實物模型輸入文件 - http://pastebin.com/AH2zxdBn

import re 
from Bio.Alphabet import generic_dna, generic_protein 
from Bio import SeqIO 

def compare_seqs(seq1, seq2): 

    similar = 0 
    diff = 0 

    diff_positions = [] 

    for i in range(0, len(seq1)): 
    if (seq1[ i ] != seq2[ i ]): 
     difference = "%s-%s [%d]" %(seq1[i], seq2[i], i) 
     diff_positions.append(difference) 
# else: 
#  similar += 1 


    return ",".join(diff_positions) 


new_seq = [] 

reference_sequence = "" 
reference_name  = "" 

outfile = open("test_out.csv", 'w') 

for record in SeqIO.parse(open('test.fa', 'ru'), 'fasta', generic_protein): 
    record_id = re.sub(r'\d+_(\d+_\d\#\d+)_\d+', r'\1', record.id) 


    if (not reference_sequence): 
     reference_sequence = record.seq 
     reference_name  = record_id 
     #continue 
    print "\t".join([reference_name, record_id, compare_seqs(reference_sequence, record.seq)]) 

這是我得到的輸出,這是不正確的,POS機454 7065_8#4其實= P

7065_8#1 7065_8#1  
7065_8#1 7065_8#2  
7065_8#1 7065_8#3  
7065_8#1 7065_8#4 E-G [245] 
7065_8#1 7065_8#5 

回答

2

解決此問題的最佳方法是將其分解爲更小的部分並驗證每個部分。

這裏有一個最小差異實現:

def compare_sequences(seq1, seq2): 
    for index, (a, b) in enumerate(zip(seq1, seq2)): 
     if a != b: 
      yield index, a, b 

這是工作:

print list(compare_sequences('abcdef', 'abddef')) 

這給了我

[(2, 'c', 'd')] 

您可以使用它作爲一個簡單的證明,它的工作。我推薦做的是將輸入放入你的函數中,並驗證它是否按預期工作。

也許有一個問題與輸入有空白或換行符,你不指望它扔掉一切?

+0

感謝你們,經過了一番,它似乎是輸入文件而不是編碼的問題。 – user3234810 2014-09-18 17:46:37