8
\$\begingroup\$

I have this function:

import re

def check_and_clean_sequence(sequence, alphabet):
    """
    Function to check and clean up all ambiguous bases in a sequence.
    Ambigous bases are bases that are not in the sequence
    alphabet, ie. 'ACGT' for DNA sequences.
    
    Inputs:
    
        sequence - a string representing a DNA sequence.
        alphabet - set/string representing the allowed
                   characters in a sequence.
    
    Outputs:
    
        cleaned_sequence - cleaned sequence or a string.
    """    
    if set(sequence).issubset(alphabet):
        return sequence
    else:
        return cleaning_ambiguous_bases(sequence)


def cleaning_ambiguous_bases(sequence):
    """
    Function to clean up all ambiguous bases in a sequence.
    Ambiguous bases are bases that are not in the sequence
    alphabet, ie. 'ACGT' for DNA sequences.
    
    Inputs:
    
        sequence - a string representing a DNA sequence.
    
    Outputs:
    
        integer - a new clean up string representing a DNA sequence
                  without any ambiguous characteres.
    """
    # compile the regex with all ambiguous bases
    pat = re.compile(r'[NRYWXSKM]')
    # look for the ambiguous bases and replace by
    # nothing
    new_seq = re.sub(pat,  '',    sequence)
    return new_seq

The performance in jupyter notebook with timeit is around:

%timeit check_and_clean_sequence(seq, iupac_dna)
200 ms ± 436 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

This was calculated with a DNA sequence with 5539776 pb or 5539.776 kpb. The sequence can be downloaded from: E. coli genome

The function receive a big string and check if there are any characters that are not in the allowed alphabet (in the case the IUPAC DNA = 'ACGT'). If the sequence has they are cleaned up and the function returns a new sequence, otherwise return the original one.

Are there are tips or tricks to improve this time?

Toy example:

>>> import re
>>> alphabet = 'ACGT'
>>> toy_sequence = 'AGCTAGGGACATTGAGGACCACCACAXYXMAAAAAA'
>>> check_and_clean_sequence(toy_sequence, alphabet)
'AGCTAGGGACATTGAGGACCACCACAAAAAAA'
\$\endgroup\$
2

3 Answers 3

7
\$\begingroup\$

From the bioinformatics side, not the python one: Your return will be non-useful for further processing whenever an ambiguous base has been present, because it changes index locations! You'll want to fix that before worrying about how long it takes to run...

DNA and RNA are themselves a form of coding language: three base units ('codons') describe amino acids. Whenever a base is deleted, a 'frameshift' occurs: every codon after that point now starts a space sooner than it should. (In vivo, a single frameshift mutation is usually sufficient to completely shut off one - or more! - genes. In other words, this is not a trivial situation.) You cannot interpret genetic sequences correctly if you change the index location of the bases!

There are three paths to usefully 'cleaning' your genome information, depending on why you wanted to do that in the first place:

  1. If space is at a premium:

    Drawing from Kraigolas: use the code as described, but rather than replacing with the empty string, replace with a uniform "this base is ambiguous" marker, thus maintaining index locations. ('N' is the standard for this specific application.)

     pattern = re.compile(r"[^ACGT]")
     print(pattern.sub("N", toy_sequence))
    

    Another solution in this vein is to store bases as tuples of (index [type:int], base type [could handle a number of ways]) but this is probably space-inefficient

  2. If neither space nor preprocessing are an optimization concern: Don't do any deletion or substitution of these characters at all.

    Look at the documentation related to the source of your sequences; there will either be an interpretation key for their specific method or they'll be using the standard lookup table. When a subsequent step will be to interpret a sequence into an amino acid string (likely application) or to look for "similar" sequences in a database (also a likely application), knowing the difference between the code that means "either A or T" and the code that means "either C or G" will save you processing time and/or get less ambiguous results.

  3. If preprocessing is a concern:

    As in 2, refer to the appropriate lookup table for what ambiguous codes mean. Read your sequence into a [length of sequence]x4 array of binary values. Column 1 is "can be A?", column 2 is "can be T?", etc. -

    Example: AGCCNCS would become

     1  0  0  0
     0  0  1  0
     0  0  0  1
     0  0  0  1
     1  1  1  1
     0  0  0  1
     0  0  1  1
    

    As appropriate codon matches (there are almost always more than one) for specific amino acids can be easily described in the same way, it's relatively trivial to use this format for a variety of purposes, and it removes assorted other headaches with accurately interpreting the ambiguous codes.

    (For examples of how to apply them - treating the bin as if they're integers, if the dot product of each row in two sequence matrixes is at least 1, you have a 'match' of codes!)

\$\endgroup\$
1
  • \$\begingroup\$ Man you totally right I just having a hard time coding this kmers count functions efficiently. I know there are a lot of tools out there but I am learning python them I need to code to learn! Then I am write my own code! But some times some obvious details get out of control. Thanks to you I avoid a big mistake. Thank you a lot. I will go back to basics and use something like if set(kmer).issubset(set(alphabet)) to avoid ambiguous kmers in the final counts. This way I will keep all the indexes corrects. Thank you a lot! Its better a slow and correct code than a quick and wrong one. Thank you. \$\endgroup\$ Jun 29, 2021 at 20:32
13
\$\begingroup\$

Perhaps there aren't enough test cases here, because I don't see why you can't just use:

pattern = re.compile(r"[^ACGT]")
print(pattern.sub("", toy_sequence))

Which says replace all characters that are not in ACGT with blanks. Thus, you do not need to know which characters are illegal, only those which are legal. This should be very fast, doesn't require a function definition because compiling pattern really serves as your function, which you now call with pattern.sub, and is a lot easier to read.

Redundant functions

If this is a satisfactory answer, I'd like to say that

def check_and_clean_sequence(sequence, alphabet):

is not a valuable function even if my solution did not work. cleaning_ambiguous_bases will achieve the same result and it will not be any slower. Checking first will at best perform the same as just calling cleaning_ambiguous_bases because regardless, you need to check every character. However, if you check first, you will iterate through the sequence potentially twice: once to check, and then once to replace. It's faster to just walk through once.

\$\endgroup\$
3
  • \$\begingroup\$ when you get so involved in coding that sometimes this kind of redundancy escape from the eyes...Thanks a lot! \$\endgroup\$ Jun 28, 2021 at 16:54
  • 1
    \$\begingroup\$ @PauloSergioSchlogl Don't worry about it! It's happened to all of us before, that's why we're here. Glad I could help :) \$\endgroup\$
    – Kraigolas
    Jun 28, 2021 at 16:56
  • \$\begingroup\$ @PauloSergioSchlogl: I did some testing with sequences that don't need cleaning, e.g. "AGCTAGGGACATTGAGGACCACCACAAAAAAAAGCTAGGGACATTGAGGACCACCACAAAAAAA". With this string, check_and_clean_sequence() is twice as fast as cleaning_ambiguous_bases() on my system (it is still slower than the suggested pattern.sub() solution, though). The reason for this is that turning the string into a set greatly reduces the comparisons that need to be made. I'd argue, then, that we can't just assume that the check is redundant as it depends on the length and the messiness of the input. \$\endgroup\$
    – Schmuddi
    Jun 29, 2021 at 15:35
8
\$\begingroup\$

There's an inconsistency between the two functions. check_and_clean_sequence() has an alphabet parameter, but this isn't used by the inner function, which has a hard-coded list of invalid characters to remove.

The check_and_clean_sequence() function doesn't add value - all it does is add overhead of an extra pass through the string (re.sub() will just do a single pass and return the original string if there's nothing to change).

In the docstrings, there are many (different!) misspellings of the word "ambiguous". I know it's not an easy word for a non-native speaker, but it shouldn't be hard to fix those.

You could possibly use filter() to keep only the characters in alphabet, but my guess is that it's unlikely to beat re for speed. re.compile() creates a simple, fast, DFA (discrete finite automaton) that's highly optimised for matching. The best improvement you could make is to compile the regex just once, rather than every call. And then as Kraigolas suggest, replace the function call with a simple pattern.sub() (I guess you could bind the empty-string first argument, to get a function that accepts just the sequence).

\$\endgroup\$
1
  • \$\begingroup\$ Thank you for your clarification! I will check for the misspellings! 8) \$\endgroup\$ Jun 28, 2021 at 17:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.