Hola, soy bastante nuevo en Python. Tengo el siguiente problema: Quiero escribir un script que, dada una secuencia (de adn) con ambigüedades, escriba todas las secuencias posibles (si hay menos de 100, si hay más de 100 secuencias posibles, se imprime un mensaje de error apropiado) Para ambigüedades de nucleótidos del ADN: http://www.bioinformatics.org/sms/iupac.html

Ejemplo: para la secuencia “AYGH”, la salida del script sería “ACGA”, “ACGC”, “ACGT”, “ATGA”, “ATGC” y “ATGT”. A, C, G y T son los nucleótidos predeterminados. TODOS los demás pueden tener valores diferentes (ver enlace).

Entonces escribí esto:

def possible_sequences (seq):
    poss_seq = ''
    for i in seq:
        if i=='A'or i=='C'or i=='G'or i=='T': 
            poss_seq += i 
        else: 
            if i== 'R':  
                poss_seq += 'A' # OR 'G', how should i implement this? 
            elif i == 'Y': 
                poss_seq += 'C' # OR T 
            elif i == 'S': 
                poss_seq += 'G' # OR C
            elif i == 'W': 
                poss_seq += 'A' # OR T 
            elif i == 'K': 
                poss_seq += 'G' # OR T
            elif i == 'M': 
                poss_seq += 'A' # OR C
            elif i == 'B': 
                poss_seq += 'C' # OR G OR T 
            elif i == 'D': 
                poss_seq += 'A' # OR G OR T 
            elif i == 'H': 
                poss_seq += 'A' # OR C OR T 
            elif i == 'V': 
                poss_seq += 'A' # OR C OR G 
            elif i == 'N': 
                poss_seq += 'A' # OR C OR G OR T 
            elif i == '-' or i == '.': 
                poss_seq += ' '
    return poss_seq

Cuando pruebo mi función: possible_sequences ('ATRY-C') obtuve:

'ATAC C'

Pero debería haber conseguido:

'ATAC C'
'ATAT C' 
'ATGC C'
'ATGT C'

¿Alguien puede ayudarme? Entiendo que tengo que recapitular y escribir un segundo poss_seq cuando hay una ambigüedad presente, pero no sé cómo ...

0
H.F.S C. 30 nov. 2016 a las 14:16
Solo está haciendo un bucle una vez. Deberá recorrer la secuencia y luego, cuando encuentre una letra para cambiar (es decir, no A, C, G o T), repita los posibles reemplazos. Tendrá que anidar estos bucles para obtener todas las permutaciones.
 – 
samb8s
30 nov. 2016 a las 14:21
Solo un consejo; puede hacer if i in 'RWMDHVN': poss_seq += 'A' en lugar de escribir o escribir varias declaraciones if-elif
 – 
Ted Klein Bergman
30 nov. 2016 a las 14:22
Sí, eso es lo que quería decir en la última frase. Lo entiendo pero no sé cómo implementarlo.
 – 
H.F.S C.
30 nov. 2016 a las 14:22

1 respuesta

La mejor respuesta

Puede utilizar itertools.product para generar las posibilidades :

from itertools import product

# List possible nucleotides for each possible item in sequence
MAP = {
    'A': 'A',
    'C': 'C',
    'G': 'G',
    'T': 'T',
    'R': 'AG',
    'Y': 'CT',
    'S': 'GC',
    'W': 'AT',
    'K': 'GT',
    'M': 'AC',
    'B': 'CGT',
    'D': 'AGT',
    'H': 'ACT',
    'V': 'ACG',
    'N': 'ACGT',
    '-': ' ',
    '.': ' '
}

def possible_sequences(seq):
    return (''.join(c) for c in product(*(MAP[c] for c in seq)))

print(list(possible_sequences('AYGH')))
print(list(possible_sequences('ATRY-C')))

Salida:

['ACGA', 'ACGC', 'ACGT', 'ATGA', 'ATGC', 'ATGT']
['ATAC C', 'ATAT C', 'ATGC C', 'ATGT C']

En lo anterior, primero iteramos sobre los elementos en la secuencia dada y obtenemos la lista de posibles nucleótidos para cada elemento:

possibilities = [MAP[c] for c in 'ATRY-C']
print(possibilities)

# ['A', 'T', 'AG', 'CT', ' ', 'C']

Luego, el iterable se desempaqueta como argumentos dados a product que devolverán el producto cartesiano:

products = list(product(*['A', 'T', 'AG', 'CT', ' ', 'C']))
print(products)

# [('A', 'T', 'A', 'C', ' ', 'C'), ('A', 'T', 'A', 'T', ' ', 'C'), 
#  ('A', 'T', 'G', 'C', ' ', 'C'), ('A', 'T', 'G', 'T', ' ', 'C')]

Finalmente, cada uno de los productos se convierte en una cadena con join:

print(list(''.join(p) for p in products))

# ['ATAC C', 'ATAT C', 'ATGC C', 'ATGT C']

Tenga en cuenta que possible_sequences devuelve un generador en lugar de construir todas las secuencias posibles a la vez para que pueda detener fácilmente la iteración cuando lo desee en lugar de tener que esperar a que se generen todas las secuencias.

7
niemmi 30 nov. 2016 a las 15:10
Todavía no he jugado bien con itertools, pero parece increíblemente poderoso. Esta parece la mejor manera de hacer esto
 – 
samb8s
30 nov. 2016 a las 14:36
Me ganaste. :) ¿Quieres que edite mi versión completa de MAP en tu respuesta? FWIW, probablemente solo haría MAP[c] en lugar de MAP.get(c, ' ') ya que probablemente queremos detectar el error de clave si obtenemos datos incorrectos.
 – 
PM 2Ring
30 nov. 2016 a las 14:49
Gracias por la oferta pero puedo hacerlo. La única razón para usar get fue que el segundo ejemplo tenía - pero en los resultados faltaba. Si es necesario detectar errores, entonces el operador de índice es el camino a seguir.
 – 
niemmi
30 nov. 2016 a las 14:52
1
Tenga en cuenta que la función possible_sequences de niemmi es un generador, por lo que si no necesita almacenar todas las combinaciones posibles en una lista, puede hacer for s in possible_sequences('AYGH'): para iterar sobre las combinaciones una por una. Obviamente, si la secuencia de entrada es grande con muchas letras que no están en 'ACGT', puede obtener muchas combinaciones. :)
 – 
PM 2Ring
30 nov. 2016 a las 14:53
1
Está bien. FWIW, puse '-': ' ', '.': ' ', en mi versión de MAP.
 – 
PM 2Ring
30 nov. 2016 a las 14:57