Как найти перевернутый повторяющийся шаблон в последовательности FASTA?

Я новичок как в Python, так и в биоинформатике, но я просматриваю веб-сайт rosalind.info, чтобы изучить и то, и другое. Вы делаете это с суффиксным деревом. Дерево суффиксов (см. http://en.wikipedia.org/wiki/Suffix_tree) волшебная структура данных, которая делает возможным все в биоинформатике. Вы быстро находите общие подстроки в нескольких длинных последовательностях. Деревьям суффиксов требуется только линейное время, поэтому возможна длина 10 000 000.

Итак, сначала найдите обратное дополнение последовательности. Затем поместите их в дерево суффиксов и найдите между ними общие подстроки (некоторой минимальной длины).

В приведенном ниже коде используется такая реализация дерева суффиксов: http://www.daimi.au.dk/~mailund/suffix_tree.html. Он написан на C с привязками к Python. Он не справится с большим количеством последовательностей, но две — не проблема. Однако я не могу сказать, будет ли это обрабатывать длину 10 000 000.

from suffix_tree import GeneralisedSuffixTree

baseComplement = { 'A' : 'T', 'C' : 'G', 'G' : 'C', 'T' : 'A' }

def revc(seq):
    return "".join([baseComplement[base] for base in seq[::-1]])

data = "AGGGTTTCCCTGACCTTCACTGCAGGTCATGCA"
# revc  TGCATGACCTGCAGTGAAGGTCAGGGAAACCCT
#       012345678901234567890123456789012
#                 1         2         3
minlength = 6

n = len(data)
tree = GeneralisedSuffixTree([data, revc(data)])
for shared in tree.sharedSubstrings(minlength):
    #print shared
    _, start, stop = shared[0]
    seq = data[start:stop]
    _, rstart, rstop = shared[1]
    rseq = data[n-rstop:n-rstart]
    print "Match: {0} at [{1}:{2}] and {3} at [{4}:{5}]".format(seq, start, stop, rseq, n-rstop, n-rstart)

Это производит вывод

Match: AGGTCA at [23:29] and TGACCT at [10:16]
Match: TGACCT at [10:16] and AGGTCA at [23:29]
Match: CTGCAG at [19:25] and CTGCAG at [19:25]

Он находит каждое совпадение дважды, по одному разу с каждого конца. А еще там есть перевернутый палиндром!

Вот реализация грубой силы, хотя она, вероятно, не очень полезна для очень длинных последовательностей.

def substrings(s, lmin, lmax):
    for i in range(len(s)):
        for l in range(lmin, lmax+1):
            subst = s[i:i+l]
            if len(subst) == l:
                yield i, l, subst

def ivp(s, lmin, lmax):
    mapping = {'A': 'T', 'G': 'C', 'T': 'A', 'C': 'G'}
    for i, l, sub in substrings(s, lmin, lmax):
        try:
            from string import maketrans
        except ImportError: # we're on Python 3
            condition = sub.translate(
                       {ord(k): v for k, v in mapping.items()})[::-1] in s
        else: # Python 2
            condition = sub.translate(maketrans('ATGC', 'TACG'))[::-1] in s
        if condition:
            yield i, l, sub

Давайте найдем «перевернутые повторяющиеся шаблоны» длины 6 (и их начальные позиции и длины):

>>> list(ivp('AGGGTTTCCCTGACCTTCACTGCAGGTCATGCA', 6, 6))
[(10, 6, 'TGACCT'), (19, 6, 'CTGCAG'), (23, 6, 'AGGTCA')]

Однако это не проверяет, перекрываются ли два шаблона. Например, 'CTGCAG' соответствует самому себе.

У меня есть идея найти перевернутую последовательность палиндромов длинной последовательности. Рассмотрите участок последовательности ДНК всей последовательности и сгенерируйте ее дополнение. Затем реверсирование участка этой последовательности комплемента. Затем выполните динамическое выравнивание между этими двумя разделами и рассчитайте его стоимость (один из фактической последовательности, а другой из последовательности с обратным дополнением). Стоимость будет давать представление о том, какое выравнивание лучше. Теперь, если стоимость наилучшего выравнивания >= пороговой стоимости, выберите этот раздел и найдите общий регион. Две общие области этого конкретного участка будут одной инвертированной повторяющейся единицей. Как только вы найдете блок, повторите его для следующего раздела, в противном случае непрерывно увеличивайте длину разделов. Кто-нибудь может реализовать этот алгоритм. Может быть, это будет полезное решение.

Я попробовал это, используя понимание списков. Я все еще новичок в python, но последние 5 лет разрабатываю C#. Это дает желаемый результат, хотя он определенно не оптимизирован для эффективной обработки строки из 10 миллионов символов.

Примечание: поскольку мы преобразовываем список в набор в частых_словах, чтобы удалить повторяющиеся записи, результаты не упорядочены.

def pattern_matching(text, pattern):
    """ Returns the start and end positions of each instance of the pattern  """
    return [[x, str(len(pattern) + x)] for x in range(len(text) - len(pattern) + 1) if
            pattern in text[x:len(pattern) + x]]


def frequent_words(text, k):
    """ Finds the most common k-mers of k """
    counts = [len(pattern_matching(text, text[x:x + k])) for x in range(len(text) - k)]
    return set([text[x:x + k] for x in range(len(text) - k) if counts[x] == max(counts)])


def reverse_complement(pattern):
    """ Formed by taking the complement of each nucleotide in Pattern """
    complements = {'A': 'T', 'C': 'G', 'T': 'A', 'G': 'C'}
    return ''.join([complements.get(c, c) for c in pattern][::-1])


def find_invert_repeats(text, pattern_size):
    """ Returns the overlap for frequent words and its reverse """
    freqs = frequent_words(text, pattern_size)
    rev_freqs = frequent_words(reverse_complement(text), pattern_size)
    return [[x, pattern_matching(text, x)] for x in set(freqs).intersection(rev_freqs)]

print(find_invert_repeats("AGGGTTTCCCTGACCTTCACTGCAGGTCATGCA", 6))

[['TGACCT', [[10, '16']]], ['AGGTCA', [[23, '29']]], ['CTGCAG', [[19, '25']]]]