2013-11-12 5 views
21

* Это краткое введение, конкретный вопрос выделен жирным шрифтом в последнем абзаце.Обратное расстояние Хэмминга

Я пытаюсь сгенерировать все строки с заданным расстоянием Хэмминга, чтобы эффективно решать задачи биоинформатики.

Идея состоит в том, что задана строка (например, «ACGTTGCATGTCGCATGATGCATGAGAGCT»), длина слова для поиска (то есть 4) и допустимые несоответствия при поиске этого слова в строке (например, 1), возвращают наиболее часто встречающиеся слова или «мутированные» слова.

Для того, чтобы быть ясно, слово длины 4 из заданной строки может быть это (между '[]'):

[ACGT]TGCATGTCGCATGATGCATGAGAGCT #ACGT 

этого

A[CGTT]GCATGTCGCATGATGCATGAGAGCT #CGTT 

этого

ACGTTGCATGTCGCATGATGCATGAG[AGCT] #AGCT 

Что я сделал (и его очень неэффективно, и его очень медленно, когда слова должны иметь 10 символов) генерируют все возможные слова с помощью заданное расстояние:

itertools.imap(''.join, itertools.product('ATCG', repeat=wordSize)) 

, а затем искать и сравнивать каждое слово в данной строке, если сгенерированные слова (или его мутация) появляется в цикле:

wordFromString = givenString[i:i+wordSize] 
mismatches = sum(ch1 != ch2 for ch1, ch2 in zip(wordFromString, generatedWord)) 
if mismatches <= d: 
    #count that generated word in a list for future use 
    #(only need the most repeated) 

То, что я хочу сделать это, вместо генерировать ВСЕ возможные слова, генерировать только мутации слов, которые появляются в данной строке с заданным числом несоответствий, другими словами, с учетом расстояния Хэмминга и слова, вернуть все возможные мутированные слова с этим (или менее) расстояние, а затем использовать их для поиска в данной строке.

Надеюсь, я был чист. Спасибо.

+1

вы бы рассчитывать перекрывается? Например, для строки «CCACCAC» и слова «CCAC» вы бы посчитали два появления или только один? так как «С» перекрывается. –

+0

Агностический ответ с объяснением http://stackoverflow.com/a/34523345/1090562 –

ответ

18
def mutations(word, hamming_distance, charset='ATCG'): 
    for indices in itertools.combinations(range(len(word)), hamming_distance): 
     for replacements in itertools.product(charset, repeat=hamming_distance): 
      mutation = list(word) 
      for index, replacement in zip(indices, replacements): 
       mutation[index] = replacement 
      yield "".join(mutation) 

Эта функция генерирует все мутации слова с Хэмминга расстояние меньше или равно заданному числу. Он относительно эффективен и не проверяет недопустимые слова. Однако действительные мутации могут появляться более одного раза. Используйте набор, если вы хотите, чтобы каждый элемент был уникальным.

+0

Так красиво и просто! Он работает отлично. – JackS

1
def generate_permutations_close_to(initial = "GACT",charset="GACT"): 
    for i,c in enumerate(initial): 
     for x in charset: 
      yield initial[:i] + x + inital[i+1:] 

будет генерировать перестановки с дист 1 (он также будет повторяться содержит повторы)

, чтобы получить набор все в пределах 2 ... затем вызвать это с каждым из первых решений, как исходные догадки

5

Пусть данное расстояние Хэмминга быть D и пусть ж быть «слово» подстроки. От ж, вы можете создать все слова с расстоянием ≤ D от глубины ограниченной depth-first search: (. Это будет отнюдь не быть быстрым, но она может служить в качестве вдохновения)

def dfs(w, D): 
    seen = set([w])  # keep track of what we already generated 
    stack = [(w, 0)] 

    while stack: 
     x, d = stack.pop() 
     yield x 

     if d == D: 
      continue 

     # generate all mutations at Hamming distance 1 
     for i in xrange(len(x)): 
      for c in set("ACGT") - set(x[i]) 
       y = x[:i] + c + x[i+1:] 
       if y not in seen: 
        seen.add(y) 
        stack.append((y, d + 1)) 

5

Если я правильно понимаю вашу проблему, вы хотите определить наивысший балл k-mers в геноме G. Оценка k-mer - это количество раз, которое оно появляется в геноме, а также количество раз, когда в геноме появляется любой k-mer с расстоянием Хэмминга менее m. Обратите внимание, что это предполагает, что вас интересуют только k-mers, которые появляются в вашем геноме (как указано @j_random_hacker).

Вы можете решить эту проблему в четыре основных этапа:

  1. Определить все k-меров в геноме G.
  2. Подсчитайте количество раз, когда каждый k-mer появляется в G.
  3. Для каждой пары (K1, K2) из к-меров, увеличиваем счетчик для обоих K1 и K2, если их расстояние Хэмминга меньше m.
  4. Найти max k-mer и его счет.

Вот пример кода Python:

from itertools import combinations 
from collections import Counter 

# Hamming distance function 
def hamming_distance(s, t): 
    if len(s) != len(t): 
     raise ValueError("Hamming distance is undefined for sequences of different lengths.") 
    return sum(s[i] != t[i] for i in range(len(s))) 

# Main function 
# - k is for k-mer 
# - m is max hamming distance 
def hamming_kmer(genome, k, m): 
    # Enumerate all k-mers 
    kmers = [ genome[i:i+k] for i in range(len(genome)-k + 1) ] 

    # Compute initial counts 
    initcount = Counter(kmers) 
    kmer2count = dict(initcount) 

    # Compare all pairs of k-mers 
    for kmer1, kmer2 in combinations(set(kmers), 2): 
     # Check if the hamming distance is small enough 
     if hamming_distance(kmer1, kmer2) <= m: 
      # Increase the count by the number of times the other 
      # k-mer appeared originally 
      kmer2count[kmer1] += initcount[kmer2] 
      kmer2count[kmer2] += initcount[kmer1] 

    return kmer2count 


# Count the occurrences of each mismatched k-mer 
genome = 'ACGTTGCATGTCGCATGATGCATGAGAGCT' 
kmer2count = hamming_kmer(genome, 4, 1) 

# Print the max k-mer and its count 
print max(kmer2count.items(), key=lambda (k,v): v) 
# Output => ('ATGC', 5) 
+2

Лучшее k-mer на самом деле не должно появляться в геноме. Например. предположим, что геном является «ACCA», k = 3 и порогом расстояния d = 1: тогда 3-мерный «CCC» появляется дважды (один раз начиная с первой позиции, начиная со второго, каждый с расстоянием Хэмминга 1), в то время как оба 3-мерных, которые появляются в геноме ('ACC' и' CCA') появляются только один раз. –

+0

@j_random_hacker: Хорошо, я обновил свой ответ, чтобы отметить, что оговорка. – mdml

+1

У меня другие серьезные проблемы, я боюсь. Если некоторый k-mer появляется m раз, тогда он будет считаться m * (m-1) раз в самовыравниваниях - так, например, если геном «AAAAA» и k = 2, ваш код скажет, что 2-мерный 'AA' появляется 12 раз! –

2

Вот что я думаю, что проблема, которую вы пытаетесь решить, состоит в следующем: у вас есть «геном» длины n, и вы хотите найти k-mer, который наиболее часто появляется в этом геноме, где «приблизительно появляется "появляется с расстоянием Хэмминга < = d. Этот k-mer не должен действительно отображаться точно в любом месте генома (например, для генома ACCA, k = 3 и d = 1, лучший k-mer - CCC, появляющийся дважды).

Если вы генерируете все километры расстояния Хэмминга < = d от некоторого k-mer в строке, а затем ищите каждый в строке, как вы сейчас делаете, то вы добавляете ненужный O (n) к времени поиска (если вы не будете искать их все одновременно, используя the Aho-Corasick algorithm, но это слишком много).

Вы можете сделать лучше, пройдя через геном, и в каждой позиции i, создавая набор всех k-mers, которые находятся на расстоянии < = d от k-mer, начиная с позиции i в геноме, и увеличивая счетчик для каждого.

0

Здесь есть правильные ответы, которые в значительной степени используют питон с его магическими функциями, которые делают почти все для вас. Я попытаюсь объяснить вещи математикой и алгоритмами, чтобы вы могли применить ее на любом желаемом языке.


Так у вас есть алфавит {a1, a2, ... a_a} (Кардинальность a) в вашем случае {'A', 'C', 'G', 'T'} и мощность равна 4. У вас есть строка длиной k и вы хотите создать все строки которых расстояние Хэмминга меньше или равным d.

Прежде всего, сколько из них у вас есть? Ответ не зависит от выбранной вами строки. Если вы выбрали укус, у вас будут C(d, k)(a-1)^d строки, из которых у вас есть расстояние от помех d от вашей строки.Таким образом, общее количество струн:

enter image description here

Он экспоненциально растет с точки зрения почти каждого параметра, так что вы не будете иметь какой-либо быстрый алгоритм, чтобы найти все слова.

Так как бы вы получили алгоритм, который будет генерировать все строки? Обратите внимание на то, что легко создать строку, которая не более чем на одном расстоянии от вашего wold. Вам просто нужно перебрать все символы в строке, и для каждого символа попробуйте каждую букву в алфавите. Как вы увидите, некоторые из слов будут одинаковыми.

Теперь, чтобы сгенерировать все строки, которые являются двумя расстояниями от помех от вашей строки, вы можете применить ту же функцию, которая генерирует одно слово hamming distance для каждого слова в предыдущей итерации.

Так вот псевдокод:

function generateAllHamming(str string, distance int): 
    alphabet = ['A', ...]// array of letters in your alphabet 
    results = {} // empty set that will store all the results 
    prev_strings = {str} // set that have strings from the previous iterations 
    // sets are used to get rid of duplicates 

    if distance > len(str){ distance = len(str)} // you will not get any new strings if the distance is bigger than length of the string. It will be the same all possible strings. 

    for d = 0; d < distance; d++ { 
     for one_string in prev_strings { 
      for pos = 0; pos < len(one_string); pos++ { 
       for char in alphabet { 
        new_str = substitute_char_at_pos(one_string, char, pos) 
        add new_str to set results 
       } 
      } 
     } 

     populate prev_strings with elements from results 
    } 

    return your results 
} 
Смежные вопросы