8

Взял курс по биоинформатике на coursera. Есть задачка, которую я не могу решить. Смысл в том что у нас есть две переменные:

  1. перечень букв в одно длинное слово, в котором есть повторяющиеся сочетания букв (напр.ACGTTGCATGTCGCATGATGCATGAGAGCT и CATG GCAT)
  2. Цифра которая задает количество повторяющихся букв(в примере выше задана цифра
  3. Надо узнать какое из сочетаний самое повторяющееся. Я нашел решение и не могу понять одну строку:

    Count = CountDict(Text, k)
    

    У меня выкидывает ошибку, а ее смысл как я понимаю в создании словаря.

    # Input:  A string Text and an integer k
    # Output: A list containing all most frequent k-mers in Text
    def FrequentWords(Text, k):
        FrequentPatterns = []
        Count = CountDict(Text, k)
        m = max(Count.values())
        for i in Count:
            if Count[i] == m:
                FrequentPatterns.append(Text[i:i+k]
        FrequentPatternsNoDuplicates = remove_duplicates(FrequentPatterns)
        return FrequentPatternsNoDuplicates
    
3
  • 4
    Перечитал вопрос 3 раза, но так и не понял суть вопроса... "Цифра которая задает количество повторяющихся букв(в примере выше задана цифра" - какая цифра? в каком примере? 25 июл 2018 в 20:47
  • я наверное неправильно выразился. Text = ACGTTGCATGTCGCATGATGCATGAGAGCT k = 4 CATG GCAT - это самые часто употребляемые комбинации 4 символов в Text.
    – Tony Stark
    25 июл 2018 в 21:56
  • связанный вопрос Algorithm to find the most common substrings in a string
    – jfs
    26 июл 2018 в 5:18

5 ответов 5

8

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

А одним из эффективных решений, особенно для больших k, в данном случае будет использование суффиксного массива в связке с определением наибольшего общего префикса двух соседних суффиксов (LCP). Построение суффиксного массива несложными средствами имеет сложность O(NlogN) или O(Nlog^2(N)). Есть и более сложные в реализации алгоритмы за линейное время.

Edit:
В коде закомментирована краткая реализацию на Python c плохой асимптотикой.
Заменил на реализацию алгоритма Manber-Myers отсюда (должна быть асимптотика O(NlogN), но, как я понял, реализация этого не обеспечивает для произвольных входных данных)

LCP строится за линейное время - я привел эффективный и простой алгоритм Kасаи.

После построения LCP нужно из него вычленить наиболее длинные серии, все элементы которых не меньше заданной длины k (думаю, на Python это делается одной строчкой). Время линейное.

Длина наибольшей серии соответствует (точнее - на единицу меньше) количеству повторов самой частой подстроки нужной длины. Для получения самой этой подстроки нужно взять подстроку, начиная с индекса из соответствующего элемента суффиксного массива.

#медленная реализация
#def get_suffix_array(s):
#    return sorted(range(len(s)), key=lambda i: s[i:])

from collections import defaultdict

def sort_bucket(s, bucket, order):
    d = defaultdict(list)
    for i in bucket:
        key = s[i + order // 2:i + order]
        d[key].append(i)
    result = []
    for k, v in sorted(d.items()):
        if len(v) > 1:
            result += sort_bucket(s, v, 2 * order)
        else:
            result.append(v[0])
    return result


def suffix_array_ManberMyers(s):
    return sort_bucket(s, range(len(s)), 1)

def lcp_kasai(s, suffarr):
    n = len(suffarr)
    k = 0
    lcp = [0] * n
    rank = [0] * n
    for i in range(n):
        rank[suffarr[i]] = i

    for  i in range(n):
        if (k>0):
            k -= 1
        if(rank[i]==n-1):
             k = 0
             continue
        j = sa[rank[i]+1]
        while((i+k<n) & (j+k<n) & (s[i+k]==s[j+k])):
            k += 1
        lcp[rank[i]] = k
    return lcp

sa = suffix_array_ManberMyers("ACGTTGCATGTCGCATGATGCATGAGAGCT$")
print(sa)
lc = lcp_kasai("ACGTTGCATGTCGCATGATGCATGAGAGCT$", sa)
print(lc)

вывод суффиксного массива и lcp:

[30, 0, 24, 26, 21, 14, 17, 7, 20, 13, 6, 11, 1, 28, 23, 25, 16, 19, 12, 5, 
 27, 9, 2, 29, 10, 22, 15, 18, 4, 8, 3]
[0, 1, 2, 1, 4, 3, 3, 0, 5, 4, 1, 2, 1, 0, 3, 2, 1, 6, 5, 2, 1, 2, 0, 1, 1,
 3, 2, 6, 2, 1, 0]

в lcp мы видим два куска длиной два (это означает, что две подстроки встречаются трижды) со значением >=4: 5, 4 и 6, 5, начинающиеся на 8 и 17 позициях. Эти позиции в суффиксном массиве содержат индексы в исходной строке 20 и 19, чему соответствуют подстроки CATG и GCAT

7
  • у вас квадратичная реализация get_suffix_array()¶ При сравнениях с реализацией в лоб, полезно посмотреть при каких больших n, алгоритм с лучшей асимптотикой становится лучше в фактических измерениях.
    – jfs
    27 июл 2018 в 18:05
  • на Питоне 2 можно было бы легко сделать алгоритм O(n log n) вместо вашего O(n * n) с помощью buffer функции, пример
    – jfs
    27 июл 2018 в 18:11
  • Да, я об этом и написал.
    – MBo
    27 июл 2018 в 18:26
  • Если у вас одна из частей квадратична, то не имеет значения, что другая линейная.
    – jfs
    27 июл 2018 в 18:28
  • Я в курсе. Добавлю реализацию алгоритма Manber-Myers
    – MBo
    27 июл 2018 в 18:46
4

Простой алгоритм в лоб: подсчитать количество вхождений для каждой [перекрывающейся] подстроки длиной k, а затем вывести подстроку с наибольшим количеством повторений:

#!/usr/bin/env python
"""Find most common substring of length k."""
from collections import Counter

text = "ACGTTGCATGTCGCATGATGCATGAGAGCT"
k = 4

# count occurrences of all k-mers: words of length k in the text
words = Counter(text[i:i+k] for i in range(len(text) - k + 1))  # O(n*k)
print(words.most_common(1)[0][0])  #-> GCAT

<script src="https://cdn.rawgit.com/brython-dev/brython/3.4.0/www/src/brython.js"></script><script src="https://cdn.rawgit.com/brython-dev/brython/3.4.0/www/src/brython_stdlib.js"></script><body onload="brython()"><script type="text/python">
from collections import Counter
from browser import document

@document["mybutton"].bind("click")
def on_click(event):
    text = document['text'].value
    k = int(document['k'].value)
    words = Counter(text[i:i+k] for i in range(len(text) - k + 1))
    print(words.most_common()[0][0])

on_click('dummy on start')
</script><label for="k">k&nbsp;=</label>&nbsp;<input id="k" value="4"> <label for="text">text&nbsp;=</label>&nbsp;<input id="text" value="ACGTTGCATGTCGCATGATGCATGAGAGCT"> <button id="mybutton">Найти k-мер</button></body>

1
  • Мне понравилась простота цикла для получения подпоследовательностей, вы не против, если я ее для моего примера в ответе использую?
    – gil9red
    26 июл 2018 в 6:08
2
from collections import defaultdict, Counter

text = 'ACGTTGCATGTCGCATGATGCATGAGAGCT'
key_len = 4

accumulator = Counter(text[i: i + key_len] for i in range(len(text) - key_len + 1))

print(accumulator)  # {'A': 7, 'C': 6, 'G': 9, 'T': 7,

max_items = defaultdict(list)

for k, v in accumulator.items():
    max_items[v].append(k)

print(max_items)  # ... , 2: ['TGCA', 'ATGA'], 3: ['GCAT', 'CATG']})

# Находим ключ с максимальным значением
max_key = max(max_items)
print(max_key, max_items[max_key])  # 3 ['GCAT', 'CATG']
4
  • спасибо, теперь понятно. Я еще не изучал collections и не знал Counter и defaultdict.
    – Tony Stark
    26 июл 2018 в 8:50
  • 1
    @AntonZubochenko Counter подсчитывает элементы коллекции, а defaultdict устанавливает значением по умолчанию при обращении по ключу, т.е. в момент вызова max_items[v] по ключу v был создан и возвращен список, если ключа не было, иначе просто возвращен список
    – gil9red
    26 июл 2018 в 9:13
  • Пока не разобрался в чем дело, но работает некорректно. Если тестить на вот этом геноме: CAGTGGCAGATGACATTTTGCTGGTCGACTGGTTACAACAACGCCTGGGGCTTTTGAGCAACGA GACTTTTCAATGTTGCACCGTTTGCTGCATGATATTGAAAACAATATCACCAAATAAATAACGC CTTAGTAAGTAGCTTTT То ответ неверный. Должен быть ТТТТ, очевидно последние 4 символа не считаются. Есть подозрение, что должно быть: accumulator = Counter(text[i: i + key_len] for i in range(len(text) - key_len+1)) т.к. нумерация идет в пайтоне с 0.
    – Tony Stark
    26 июл 2018 в 12:41
  • 1
    @AntonZubochenko Посмотрю в чем дело, в ответе jfs тоже, получается такая-же проблема? Попробуйте: accumulator = Counter(text[i: i + key_len] for i in range(len(text) - key_len + 1)). Обновил ответ, вроде бы теперь правильно
    – gil9red
    26 июл 2018 в 12:55
0

Попробуйте нечто вроде:

def CountDict(Text, k):
    # Создаём полный список всех подстрок длиной k
    if k > len(Text): return None
    if k == len(Text): return {Text:1}

    list_substr = list()   # Список подстрок длиной k
    for j in range(len(Text) - k + 1):
        list_substr.append( Text[j:j+k])

    # Список превращаем в словарь, с подсчётом частот
    RetMap = dict()
    for j in range(len(list_substr)):
        key_val = list_substr[j]
        RetMap[key_val] = RetMap.get(key_val, 0) + 1

    return RetMap

dnk_map = CountDict('ACGTTGCATGTCGCATGATGCATGAGAGCT', 4)
print(dnk_map)
-1

еще один вариант (сложность алгоритма затрудняюсь оценить):

from itertools import product
text = "ACGTTGCATGTCGCATGATGCATGAGAGCT"

cnt = [[i,text.count(''.join(i))] for i in product('ACGT',repeat=4)]
m = max(cnt,key=lambda x: x[1])[1]
patterns = [''.join(i) for i,_ in filter(lambda x: x[1]==m,cnt)]

print(patterns) # ['CATG', 'GCAT']
1
  • 1
    добро пожаловать на stack overflow на русском! если у вас появился новый вопрос, задайте его, пожалуйста, с помощью кнопки «задать вопрос» в правом верхнем углу страницы; если нужно указать контекст, дайте ссылку на этот вопрос. а этот «ответ» пожалуйста удалите. 19 авг 2021 в 11:48

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.