python - Inverso de la distancia de Hamming

string algorithm

Question

* Esta es una breve introducción, la pregunta específica está en negrita en el último párrafo.

Estoy tratando de generar todas las cadenas con una distancia de Hamming dada para resolver de manera eficiente una asignación bioinformática.

La idea es, dada una cadena (es decir, ''ACGTTGCATGTCGCATGATGCATGAGAGCT''), la longitud de la palabra a buscar (es decir, 4) y las discrepancias aceptables al buscar esa palabra en la cadena (es decir, 1), devolver las palabras más frecuentes o palabras ''mutadas''.

Para que quede claro, una palabra de longitud 4 de la cadena dada puede ser esto (entre ''[]''):

[ACGT]TGCATGTCGCATGATGCATGAGAGCT #ACGT

esta

A[CGTT]GCATGTCGCATGATGCATGAGAGCT #CGTT

o esto

ACGTTGCATGTCGCATGATGCATGAG[AGCT] #AGCT

Lo que hice fue (y es muy ineficiente, y es realmente lento cuando las palabras necesitan tener 10 caracteres) generar todas las palabras posibles con la distancia dada:

itertools.imap(''''.join, itertools.product(''ATCG'', repeat=wordSize))

y luego busque y compare cada palabra en la cadena dada si las palabras generadas (o su mutación) aparecen en un bucle:

wordFromString = givenString[i:i+wordSize]
mismatches = sum(ch1 != ch2 for ch1, ch2 in zip(wordFromString, generatedWord))
if mismatches <= d:
    #count that generated word in a list for future use
    #(only need the most repeated)

Lo que quiero hacer es, en lugar de generar TODAS las palabras posibles, generar solo las mutaciones de las palabras que aparecen en la cadena dada con un número dado de desajustes, en otras palabras, dada una Distancia de Hamming y una palabra, devolver todo lo posible palabras mutadas con esa (o menos) distancia , y luego úselas para buscar en la cadena dada.

Espero haber sido claro. Gracias.

Answer 1

Aquí hay respuestas correctas, que utilizan Python en gran medida con sus funciones mágicas que hacen casi todo por ti. Trataré de explicar las cosas con matemática y algoritmos, para que puedas aplicarlo a cualquier idioma que desees.

Entonces tiene un alfabeto {a1, a2, ... a_a} (la cardinalidad de a ) en su caso {''A'', ''C'', ''G'', ''T''} y la cardinalidad es 4. Usted tiene un cadena de longitud k y desea generar todas las cadenas cuya distancia de hamming es menor o igual a d .

En primer lugar, ¿cuántos de ellos tienes? La respuesta no depende de la cadena que seleccione. Si seleccionó una picadura, tendrá C(d, k)(a-1)^d cuerdas de las cuales tienen una distancia de Hamming d de su cadena. Así que el número total de cadenas es:

Aumenta exponencialmente en términos de casi todos los parámetros, por lo que no tendrá ningún tipo de algoritmo rápido para encontrar todas las palabras.

Entonces, ¿cómo derivarías un algoritmo que generará todas las cadenas? Tenga en cuenta que es fácil generar una cadena que se encuentra a una distancia máxima de Hamming de su wold. Solo necesitas iterar sobre todos los caracteres en la cadena y para cada carácter prueba cada letra en el alfabeto. Como verás, algunas de las palabras serían iguales.

Ahora, para generar todas las cadenas que están a dos distancias de hamming de su cadena, puede aplicar la misma función que genera una palabra de distancia de hamming para cada palabra en la iteración anterior.

Así que aquí hay un pseudocódigo:

function generateAllHamming(str string, distance int): 
    alphabet = [''A'', ...]// array of letters in your alphabet
    results = {} // empty set that will store all the results
    prev_strings = {str} // set that have strings from the previous iterations
    // sets are used to get rid of duplicates

    if distance > len(str){ distance = len(str)} // you will not get any new strings if the distance is bigger than length of the string. It will be the same all possible strings.

    for d = 0; d < distance; d++ {
        for one_string in prev_strings {
            for pos = 0; pos < len(one_string); pos++ {
                for char in alphabet {
                    new_str = substitute_char_at_pos(one_string, char, pos)
                    add new_str to set results 
                }
            }
        }

        populate prev_strings with elements from results
    }

    return your results
}

Answer 2

Deje que la distancia de Hamming dada sea D y sea w la subcadena de "palabra". Desde w , puede generar todas las palabras con distancia ≤ D mediante una búsqueda de profundidad limitada con profundidad limitada:

def dfs(w, D):
    seen = set([w])      # keep track of what we already generated
    stack = [(w, 0)]

    while stack:
        x, d = stack.pop()
        yield x

        if d == D:
            continue

        # generate all mutations at Hamming distance 1
        for i in xrange(len(x)):
            for c in set("ACGT") - set(x[i])
                 y = x[:i] + c + x[i+1:]
                 if y not in seen:
                     seen.add(y)
                     stack.append((y, d + 1))

(Esto no será rápido, pero puede servir como inspiración).

Answer 3

Esto es lo que creo que el problema que intentas resolver es: tienes un "genoma" de longitud n, y quieres encontrar el k-mer que aparece aproximadamente con mayor frecuencia en este genoma, donde "aproximadamente aparece" significa que aparece con distancia de Hamming <= d. No es necesario que este k-mer aparezca exactamente en cualquier lugar del genoma (por ejemplo, para el genoma ACCA , k = 3 y d = 1, el mejor k-mer es CCC , que aparece dos veces).

Si genera todos los k-mers de la distancia de Hamming <= d desde algún k-mer en la cadena y luego busca cada uno en la cadena, como parece estar haciendo actualmente, entonces está agregando una O (n) innecesaria factor al tiempo de búsqueda (a menos que los busque a todos simultáneamente usando el algoritmo Aho-Corasick , pero eso es una exageración aquí).

Puede hacerlo mejor si recorre el genoma y, en cada posición i, genera el conjunto de todos los k-mers que están a una distancia <= d desde el k-mer que comienza en la posición i en el genoma, e incrementa un contador para cada uno uno.

Answer 4

Si entiendo tu problema correctamente, quieres identificar la puntuación más alta de k-mers en un genoma G La puntuación de un k-mer es el número de veces que aparece en el genoma más el número de veces que cualquier k-mer con una distancia de Hamming menor que m también aparece en el genoma. Tenga en cuenta que esto supone que solo está interesado en los k-mers que aparecen en su genoma (como lo señala @j_random_hacker).

Puedes resolver esto en cuatro pasos básicos:

1. Identificar todos los k-mers en el genoma G
2. Cuente el número de veces que cada k-mer aparece en G
3. Para cada par ( K1 , K2 ) de k-mers, incremente la cuenta para K1 y K2 si su distancia de Hamming es menor que m .
4. Encuentra el max k-mer y su cuenta.

Aquí está el ejemplo del código de Python:

from itertools import combinations
from collections import Counter

# Hamming distance function
def hamming_distance(s, t):
    if len(s) != len(t):
        raise ValueError("Hamming distance is undefined for sequences of different lengths.")
    return sum( s[i] != t[i] for i in range(len(s)) )

# Main function
# - k is for k-mer
# - m is max hamming distance
def hamming_kmer(genome, k, m):
    # Enumerate all k-mers
    kmers = [ genome[i:i+k] for i in range(len(genome)-k + 1) ]

    # Compute initial counts
    initcount  = Counter(kmers)
    kmer2count = dict(initcount)

    # Compare all pairs of k-mers
    for kmer1, kmer2 in combinations(set(kmers), 2):
        # Check if the hamming distance is small enough
        if hamming_distance(kmer1, kmer2) <= m:
            # Increase the count by the number of times the other
            # k-mer appeared originally
            kmer2count[kmer1] += initcount[kmer2]
            kmer2count[kmer2] += initcount[kmer1]

    return kmer2count


# Count the occurrences of each mismatched k-mer
genome = ''ACGTTGCATGTCGCATGATGCATGAGAGCT''
kmer2count = hamming_kmer(genome, 4, 1)

# Print the max k-mer and its count
print max(kmer2count.items(), key=lambda (k,v ): v )
# Output => (''ATGC'', 5)

Answer 5

def generate_permutations_close_to(initial = "GACT",charset="GACT"):
    for i,c in enumerate(initial):
         for x in charset:
             yield initial[:i] + x + inital[i+1:]

generará permutaciones con un dist de 1 (también repetirá contener repeticiones)

para obtener un conjunto de todo dentro de 2 ... luego llámelo con cada una de las primeras soluciones como conjeturas iniciales

Answer 6

def mutations(word, hamming_distance, charset=''ATCG''):
    for indices in itertools.combinations(range(len(word)), hamming_distance):
        for replacements in itertools.product(charset, repeat=hamming_distance):
            mutation = list(word)
            for index, replacement in zip(indices, replacements):
                mutation[index] = replacement
            yield "".join(mutation)

Esta función genera todas las mutaciones de una palabra con una distancia de hamming menor o igual a un número dado. Es relativamente eficiente y no revisa palabras inválidas. Sin embargo, las mutaciones válidas pueden aparecer más de una vez . Usa un conjunto si quieres que cada elemento sea único.