Faster replace method

user_12670608 0
Viewed 113

I have written a function, which takes aligned DNA sequences and replaces "-" with "Z" if the "-" is both preceded and followed by other bases. The goal is to make those bases a "Z" so I can differentiate unsequenced regions of a genome from insertions/deletions. Here's the function:

  def find_insertion_deletion(sequence):
    pattern = r'[A-Z]+-+(?=[A-Z]+)'
    new_sequence = re.sub(pattern, lambda x: x.group().replace('-', 'Z'), sequence)
    return new_sequence

This is obviously using regex to find the pattern and then make the replacements. He's an example of how it works.

sequence = '-----------AGCATCGACGTCTAGTAC---CGTACGTA--CGTACGTAGCTA-GCTAGCTAGCTGATCGATGCTAGCA---------------'
new_sequence = find_insertion_deletion(sequence)

output:

new_sequence = '-----------AGCATCGACGTCTAGTACZZZCGTACGTAZZCGTACGTAGCTAZGCTAGCTAGCTGATCGATGCTAGCA---------------'

It works exactly how I want it to, but it's very slow when running this on many sequences in an alignment. Is there a way I can speed this up by a significant amount? I assumed regex would be the fastest way, but maybe there's another way I'm not aware of..

Thanks!

python string
1 Answers

I'll assume that your example is representative, i.e., that you simply want to replace all hyphens except those at the start and at the end. This uses basic string functions, which are much faster:

def find_insertion_deletion(sequence):
    stripped = sequence.strip('-')
    if not stripped:
        return sequence
    start = sequence.index(stripped[0])
    end = len(sequence) - start - len(stripped)
    return '-' * start + stripped.replace('-', 'Z') + '-' * end

Or:

def find_insertion_deletion(sequence):
    lstripped = sequence.lstrip('-')
    start = len(sequence) - len(lstripped)
    stripped = lstripped.rstrip('-')
    end = len(lstripped) - len(stripped)
    return '-' * start + stripped.replace('-', 'Z') + '-' * end

Benchmarks

With your example sequence:

   80.5 us  original
   20.5 us  Wiktor_1
   14.4 us  Wiktor_2
    3.6 us  Kelly_1
    3.3 us  Kelly_2

With a longer sequence (sequence *= 1000):

 5931.9 us  original
20896.0 us  Wiktor_1
 7498.8 us  Wiktor_2
  150.5 us  Kelly_1
  160.9 us  Kelly_2

Code:

from timeit import repeat
import re
import regex

def original(sequence):
    pattern = r'[A-Z]+-+(?=[A-Z]+)'
    new_sequence = re.sub(pattern, lambda x: x.group().replace('-', 'Z'), sequence)
    return new_sequence

def Wiktor_1(sequence):
    return regex.sub(r'(?:\G(?!\A)|[A-Z](?=-+[A-Z]))\K-', 'Z', sequence)

def Wiktor_2(sequence):
    return re.sub(r'\b-+\b', lambda x: x.group().replace('-', 'Z'), sequence)

def Kelly_1(sequence):
    stripped = sequence.strip('-')
    if not stripped:
        return sequence
    start = sequence.index(stripped[0])
    end = len(sequence) - start - len(stripped)
    return '-' * start + stripped.replace('-', 'Z') + '-' * end

def Kelly_2(sequence):
    lstripped = sequence.lstrip('-')
    start = len(sequence) - len(lstripped)
    stripped = lstripped.rstrip('-')
    end = len(lstripped) - len(stripped)
    return '-' * start + stripped.replace('-', 'Z') + '-' * end

funcs = original, Wiktor_1, Wiktor_2, Kelly_1, Kelly_2

sequence = '-----------AGCATCGACGTCTAGTAC---CGTACGTA--CGTACGTAGCTA-GCTAGCTAGCTGATCGATGCTAGCA---------------'
sequence *= 1   # or 1000 with number = 10**2
number = 10**5

expect = original(sequence)
for func in funcs:
    print(func(sequence) == expect, func.__name__)

for _ in range(3):
    print()
    for func in funcs:
        t = min(repeat(lambda: func(sequence), number=number)) / number
        print('%7.1f us ' % (t * 1e6), func.__name__)
user_126710570
Related