Найти все (k, d) -мотифы в наборе строк ДНК

Question

В основном, вопрос заключается в том, чтобы попросить выяснить все возможные мотивы (k-mers long) с не более чем d-несоответствиями среди коллекции строк ДНК. Я могу написать код ниже, чтобы найти все мотивы (k, d) для одной строковой ДНК. Я не знаю, как изменить свой код, когда выйдет многострочная строка ДНК.

Пример ввода:

к = 3, г = 1

ATTTGGC

TGCCTTA

CGGTATC

GAAAATT

Пример вывода:

ATA

ATT

GTT

ТТТ

import collections 

    kmer = 5; 
    in_genome = "GGGGCTTCACAGCGCCCCTACAATACAATAGCCCTCGAATACCTACTTGCCACTATGTTCGGCGTCATTACATACGACCCGCATGCTCGGCAGTATGTCTCTACTCAGGATCCCTCAATATTACTTACGCCAATATGTCTAAGGTTTAGA"; 
    in_mistake = 1; 
    out_result = []; 
    mismatch_list = [] 

    def hamming_distance(s1, s2): 
     # Return the Hamming distance between equal-length sequences 
     if len(s1) != len(s2): 
      raise ValueError("Undefined for sequences of unequal length") 
     else: 
      return sum(ch1 != ch2 for ch1, ch2 in zip(s1, s2)) 

    for i in xrange(len(in_genome)-kmer + 1): 
     v = in_genome[i:i + kmer] 
     out_result.append(v) 

    for t_kmer in set(out_result): 
     for s_kmer in out_result: 
      if hamming_distance(t_kmer, s_kmer) <= in_mistake: 
       mismatch_list.append(t_kmer) 

    mismatch_count = collections.Counter(mismatch_list) 

    print mismatch_count 

Answer 1

Проблема заключается в том, чтобы переключить код с помощью внутренних переменных на чтение ввода из файла. Вы не можете просто конкатенировать нити ДНК файла и запускать их по-прежнему, поскольку это изменяет результат, где встречаются концы нитей. Вы также должны обрабатывать первую строку ввода по-разному от других, поскольку она содержит параметры программы, а остальные - необработанные данные:

import re 
import sys 
import collections 

mismatch_list = [] 

def hamming_distance(s1, s2): 
    """ Returns the Hamming distance between equal-length sequences """ 
    if len(s1) != len(s2): 
     raise ValueError("Undefined for sequences of unequal length") 
    return sum(ch1 != ch2 for ch1, ch2 in zip(s1, s2)) 

with open(sys.argv[1]) as file: 

    kmer = None 
    in_mistake = None 

    parameters = file.readline().rstrip() # first line of file has parameters 

    matchobj = re.search(r"k\s*=\s*(\d+)", parameters) 
    if matchobj: 
     kmer = int(matchobj.group(1)) 

    matchobj = re.search(r"d\s*=\s*(\d+)", parameters) 
    if matchobj: 
     in_mistake = int(matchobj.group(1)) 

    assert kmer is not None and in_mistake is not None, "file parameters misread" 

    for sequence in file: # subsequent lines of file are DNA strands 
     sequence = sequence.rstrip() 
     if not sequence: 
      continue # ignore blank lines 

     result = [] 

     for i in range(len(sequence) - kmer + 1): 
      v = sequence[i:i + kmer] 
      result.append(v) 

     for t_kmer in set(result): 
      for s_kmer in result: 
       if hamming_distance(t_kmer, s_kmer) <= in_mistake: 
        mismatch_list.append(t_kmer) 

mismatch_count = collections.Counter(mismatch_list) 

print(mismatch_count)