I'm doing some coding with DNA sequences and I'm interested in a function to find sequential repeats (which could represent where primers could 'slip' AKA do bad stuff).
An example of what I'm interested in would be as follows:
longest_repeat('ATTTTCCATGATGATG')
which would output the repeat length and coordinates, in this case 9 long and 7:15. The function should have picked up the ATGATGATG at the end and since it is longer than the TTTT repeat and the TGATGA repeat, it would only report the ATGATGATG. In the case of ties, I'd like if it could report all the ties, or at least one of them.
It would also be nice to set a threshold to only report these sequential repeats if they're over a specific length.
I have some experience in python, but this specific question has me stumped, since if I code it inefficiently and put in a 50 character long string it could take forever. I appreciate all the help!
Here is a solution:
def longest_repeat(seq, threshold):
results = []
longest = threshold
# starting position
for i in range(len(seq)):
# pattern period
for p in range(1, (len(seq)-i)//2+1):
# skip unecessary combinations
if results != [] and results[-1][0] == i and results[-1][3] % p == 0: continue
# max possible number of repetitions
repetitions = len(seq)//p
# position within the pattern's period
for k in range(p):
# get the max repetitions the k-th character in the period can support
m = 1
while i+k+m*p < len(seq) and seq[i+k] == seq[i+k+m*p]:
m += 1
repetitions = min(m, repetitions)
# check if we're already below the best result so far
if repetitions*p < longest: break
# save the result if it's good
if repetitions > 1 and repetitions*p >= longest:
# overwrite lesser results
if repetitions*p > longest: results = []
# store the current one (with ample information)
results += [(i, seq[i:i+p], repetitions, repetitions*p)]
longest = max(longest, repetitions*p)
return results
The logic is that you run through each starting position in the sequence (i), you check every sensible pattern period (p) and for that combination you check if they result in a substring at least as good as the best one so far (or the threshold, if no result has been found yet).
The result is a list of tuples of the form (starting index, period string, repetitions, total length). Running your example
threshold = 5
seq = 'ATTTCCATGATGATG'
t = time.time()
results = longest_repeat(seq, threshold)
print("execution time :", time.time()-t)
for t in results:
print(t)
we get
exec : 0.00010848045349121094
(6, 'ATG', 3, 9)
From there, it is trivial to get the full matched string (simply do period_string * repetitions)
For a random input of 700 characters, the execution time is ~6.8 seconds, compared to ~20.2 seconds using @IoaTzimas's answer.