2014-09-18 3 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)]) 

Вот результат я получаю, что неправильно, как поз 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

Смежные вопросы