You should try MOODS: it's a suite of algorithms for matching position weight matrices (PWM) against DNA sequences. I use it on a daily basis, it's a very good piece of software written in C++ with interface for python (a simple import MOODS, and up you go!). The must difficult step is to convert your query sequences in PWM. Here is the function I use:

def primer2pwm(primer):
    """
    Write a primer sequence as a position weight matrix.
    """
    # Create 4 lists of length equal to primer's length.
    matrix = [[0] * len(primer) for i in range(4)]
    # List of correspondance IUPAC.
    IUPAC = {
        "A" : ["A"], 
        "C" : ["C"], 
        "G" : ["G"],        
        "T" : ["T"],        
        "U" : ["U"],        
        "R" : ["G", "A"], 
        "Y" : ["T", "C"], 
        "K" : ["G", "T"], 
        "M" : ["A", "C"], 
        "S" : ["G", "C"], 
        "W" : ["A", "T"], 
        "B" : ["C", "G", "T"], 
        "D" : ["A", "G", "T"], 
        "H" : ["A", "C", "T"], 
        "V" : ["A", "C", "G"], 
        "N" : ["A", "C", "G", "T"]
    }
    # Position of nucleotides in the PWM.
    dico = {"A" : 0,  "C" : 1, "G" : 2, "T" : 3}
    # Read each IUPAC letter in the primer.
    for index, letter in enumerate(primer):
        for nuc in IUPAC.get(letter):
            i = dico.get(nuc)
            matrix[i][index] = 1
    return matrix

Perhaps it would also be helpful to use the "matching core" of cutadapt?

It's used to strip adapters from sequencing data, which sounds like a problem very much (in spirit) to what you are trying to solve.

There is a fast Levenshtein edit distance module for Python here:

http://code.google.com/p/pylevenshtein/

Have a look at Pyhton's difflib. You can use it for approximate pattern matching.

An example taken from there:

s = SequenceMatcher(None, "abcd", "bcde")
>>> s.ratio()
0.75

The first argument to SequenceMatcher allows you to ignore certain characters, which might be handy for ambiguity characters like N. If you want to search on the reverse complement as well, then you will have to create it (e.g. using Biopython's reverse_complement())

Matching against billions of sequences will be slow though. Two additional functions (quick_ratio() and real_quick_ratio()) might help a bit, but if calling an external program is an option for you, then have a look at e.g. Mummer or similar programs.

Andreas