How to find inverted repeated pattern in a FASTA sequence?

Question

Suppose my long sequence looks like:

5’-AGGGTTTCCC**TGACCT**TCACTGC**AGGTCA**TGCA-3

The two italics subsequences (here within the two stars) in this long sequence are together called as inverted repeat pattern. The length and the combination of the four letters such as A,T,G,C in those two subsequences will be varying. But there is a relation between these two subsequence. Notice that, when you consider the first subsequence then its complementary subsequence is ACTGGA (according to A combines with T and G combine with C) and when you invert this complementary subsequence (i.e. last letter comes first) then it matches with the second subsequence.

There are large number of such patterns present in a FASTA sequence (contains 10 million ATGC letters) and I want to find such patterns and their start and end positions.

Answer 1

I'm new to both Python and bioinformatics, but I'm working my way through the rosalind.info web site to learn some of both. You do this with a suffix tree. A suffix tree (see http://en.wikipedia.org/wiki/Suffix_tree) is the magical data structure that makes all things possible in bioinformatics. You quickly locate common substrings in multiple long sequences. Suffix trees only require linear time, so length 10,000,000 is feasible.

So first find the reverse complement of the sequence. Then put both into the suffix tree, and find the common substrings between them (of some minimum length).

The code below uses this suffix tree implementation: http://www.daimi.au.dk/~mailund/suffix_tree.html. It's written in C with Python bindings. It won't handle a large number of sequences, but two is no problem. However I can't say whether this will handle length 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)

This produces output

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]

It finds each match twice, once from each end. And there's a reverse palindrome in there, too!

Answer 2

Here's a brute force implementation, although it's probably not very useful on super long sequences.

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

Let's find "inverted repeat patterns" of length 6 (and their start positions and lengths):

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

This doesn't check if the two patterns overlap, though. For instance, 'CTGCAG' matches itself.

Answer 3

I have an idea for finding the inverted palindrome sequence of a long sequence. Consider a section of a DNA sequence of the whole sequence and generate its complement. Then reversing the section of this complement sequence. Then perform the dynamical alignment between these two sections and calculate the cost of it(one from actual sequence and other from reverse complement sequence). The cost will be giving an idea which alignment is best. Now if the cost of the best alignment >=threshold cost then select that section and find the common region. The two common regions of this particular section will be one inverted repeat unit. Once find the unit then repeat it for the next section other wise increase the length of the sections continuously. Can anyone implement this algorithm. May be it will be a useful solution.

Answer 4

I took a shot at it using list comprehensions. I'm still new to python, but have been a C# dev for the last 5 or so years. This produces the output you desire, although it's definitely not optimized to be able to efficiently handle a 10 million character string.

Note: because we convert the list into a set in frequent_words, in order to remove duplicate entries, the results are not ordered.

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']]]]