• Publicidad

Localizar cadenas de ADN con comodines

Perl aplicado a la bioinformática

Localizar cadenas de ADN con comodines

Notapor danusol » 2010-12-03 05:36 @275

Hola perleros,

tengo un archivo de texto con secuencia de ADN en el que tengo que buscar varias cadenas dadas como argumento al programa, por ejemplo CAGTT y ACGTA. Estas cadenas las guardo en un array y en un bucle recorro el array buscando en el texto todas las apariciones de cada una.

De la siguiente manera abreviada: (nota: he usado la librería Tie para poder leer el archivo de texto como un array aunque requiere mucha memoria porque los archivos son de cientos de megas)

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
tie my @array, 'Tie::File', $ARGV[0] or die;
my $tags = <STDIN>;
my @tags = split (',',$tags);
for (@array)
        {
                if (/^$tags[$t]/i)
               

 
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4


Si ahora quiero permitir que la primera base de cada cadena pueda ser una indeterminación N, es decir que me busque la presencia de CAGTT o NAGTT, ¿cómo lo podría hacer?

Editado: de hecho, me gustaría que esa indeterminación pudiese tener lugar en cualquiera de las posiciones pero sólo en una, es decir, permito NAGTT o CANTT o CAGTN o ..., pero no CNNTT.

Mil gracias.
danusol
Perlero nuevo
Perlero nuevo
 
Mensajes: 46
Registrado: 2010-04-22 07:08 @339

Publicidad

Re: Localizar cadenas de ADN con comodines

Notapor explorer » 2010-12-03 08:23 @391

Yo no lo habría hecho así... Tie se habría leído el fichero de texto completo a memoria, y si solo me interesa localizar las cadenas, línea a línea, sin importar en qué línea están, entonces es mejor hacer un bucle por todo el fichero, línea a línea.

El proceso -que yo seguiría- es:
* abrir el fichero de secuencias y lo leería línea a línea
* por cada línea, miraría si están las cadenas que estoy buscando. Si están, pinto toda la secuencia.

No nos dices ni muestras qué quieres hacer cuando encuentras esas cadenas, así que he supuesto que solo queremos pintarlas.

Para el caso de los comodines, e incluso para el caso de buscar más de una cadena a la vez, pues Perl viene fenómeno, pues podemos construir una expresión regular desde cero, con todas esas alternativas.

Algo así (no probado):
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;
  4. use diagnostics;
  5.  
  6. my $fichero = $ARGV[0] or die "Uso: $0 <fichero de secuencias>\n";
  7.  
  8. print "Introduzca la cadena o cadenas a buscar, separadas por comas: ";
  9. my $cadenas = <>;                   # leemos las cadenas
  10. $cadenas =~ s/\s+//g;               # quitamos todos los caracteres en blanco y el final de línea
  11. $cadenas = uc $cadenas;             # todo a mayúsculas, por si acaso
  12.  
  13. for (split /,/, $cadenas) {         # para cada cadena
  14.     my $enes = tr/N/N/;             # contamos el número de 'N' que tiene
  15.     die "ERROR: No se puede indicar más de un comodín en cada cadena. Cadena errónea: [$_]\n"
  16.         if $enes > 1;
  17. }
  18.  
  19. $cadenas =~ tr/,N/|./;              # convertimos las comas en alternativas y las 'N' en comodines
  20.  
  21. my $regex = qr/(?:$cadenas)/;       # construimos la expresión regular
  22.  
  23. open my $fh, q[<], $fichero  or  die "ERROR: No puedo leer el fichero $fichero: $!\n";
  24. while (<$fh>) {
  25.     print if $regex;                # pintamos la línea si contiene alguna de las cadenas
  26. }
  27. close $fh;
  28.  
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


Lo que me despista es tu '^' en tu expresión regular... estás reduciendo la búsqueda de la cadena a solo en el comienzo de la secuencia.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Localizar cadenas de ADN con comodines

Notapor danusol » 2010-12-03 08:46 @406

Hola explorer muchas gracias una vez más por tu ayuda. Me extiendo más.
El archivo a parsear tiene bloques de cuatro líneas de la siguiente manera, a cada secuencia de DNA con su encabezado le corresponde una secuencia de la misma longitud de caracteres ASCII con su encabezado, y hay muchos bloques:

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
@SEQ1
CAGTTTCGGGTTTTATGCC
+SEQ1
+++-)+++++abtttbcbm
@SEQ2
ACGTATCGGGTCCTATTCC
+SEQ2
+n+-n++b++abtntbccm
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


y tengo que buscar estas cadenas de cinco bases, en total 10 cadenas, permitiendo que haya una única indeterminación en cualquier posición. Cuando encuentre la primera cadena con o sin indeterminación quiero que me elimine esa cadena de la línea 2 del bloque y los cinco caracteres correspondientes de la línea 4 del bloque. Quedaría de la siguiente manera en el ejemplo buscando CAGTT y ACGTA:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
@SEQ1
TCGGGTTTTATGCC
+SEQ1
+++++abtttbcbm
@SEQ2
TCGGGTCCTATTCC
+SEQ2
++b++abtntbccm
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Y quiero que me guarde todos los bloques específicos de cada cadena en un archivo diferente, de tal manera que si tengo diez cadenas originales, acabe con diez archivos de salida en los que se han eliminado esas cadenas y sus correspondientes caracteres ASCII.

El hecho de tener que meter las cadenas con sus posibles opciones de indeterminación como argumento daría lugar a tener que meter cincuenta cadenas en el argumento. Esto aunque tedioso se podría llegar a hacer, pero entonces acabaría con 50 archivos y luego tendría que unir los cinco correspondientes a cada cadena.

Espero haberme explicado mejor.

Gracias por tu paciencia.
danusol
Perlero nuevo
Perlero nuevo
 
Mensajes: 46
Registrado: 2010-04-22 07:08 @339

Re: Localizar cadenas de ADN con comodines

Notapor explorer » 2010-12-03 11:55 @538

Pero, ¿las cadenas de búsqueda de 5 bases solo las buscamos al principio de las secuencias? ¿No pasa nada porque aparezcan en mitad de las secuencias?
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Localizar cadenas de ADN con comodines

Notapor danusol » 2010-12-03 12:11 @549

Hola de nuevo, experto.

No, solo pueden estar al principio de la secuencia, por eso pongo el símbolo ^.

Se me ha ocurrido, para no tener que introducir todas las opciones, crear dentro de cada bucle que busque cada cadena, primero generar cinco variables con las posibilidades de indeterminación, usando las funciones substr() y tr///, pero substr() no me selecciona la letra que quiero.

Como ejemplo pongo cambiar la primera y la segunda base por N para una de las cadenas, y quedaría algo así como:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
my $cadena_v1 = $cadena;
$cadena_v1=~tr/substr $cadena_v1,0,1/N/; #cambiar la primera base por N, pero no me funciona
my $cadena_v2 = $cadena;
$cadena_v2=~tr/substr $cadena_v1,1,1/N/; #cambiar la segunda base por N, pero no me funciona


if (/^$cadena/i) || (/^$cadena_v1/i) || (/^cadena_v2/i)
....

 
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


Gracias
danusol
Perlero nuevo
Perlero nuevo
 
Mensajes: 46
Registrado: 2010-04-22 07:08 @339

Re: Localizar cadenas de ADN con comodines

Notapor explorer » 2010-12-03 12:31 @563

O sea... recapitulando...

Tenemos un fichero con un montón de bloques, cada bloque son 4 líneas, de las que solo nos interesan la segunda y la cuarta.

Le pedimos al usuario que nos dé una o más cadenas de búsqueda de 5 bases cada una.

Por cada cadena, generamos seis combinaciones, que son la misma cadena y cinco más, de resultas de una indeterminación en cada base.

Cada combinación la buscamos en los comienzos de las secuencias (segundas líneas). Si están, recortamos las comienzos de las líneas 2 y las líneas 4, y grabamos el bloque resultado en un fichero aparte (del que no sabemos cómo nombrarlo y tampoco qué hacer con los bloques que no concuerdan con ninguna cadena de búsqueda).

¿Es eso?

Este programa da un resultado así, excepto la parte de creación de los ficheros.
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use common::sense;
  3. use autodie;
  4.  
  5. ## Pedimos las cadenas a buscar
  6. print "Introduzca la cadena o cadenas a buscar, separadas por comas: ";
  7. my $cadenas = <>;                      # leemos las cadenas
  8. #my $cadenas = 'CAGTT,ACGTA';
  9.  
  10. $cadenas =~ s/\s+//g;                   # quitamos todos los caracteres en blanco y el final de línea
  11. $cadenas = uc $cadenas;                 # todo a mayúsculas, por si acaso
  12.  
  13. my @cadenas = split /,/, $cadenas;
  14.  
  15. ## Generamos las combinaciones por cada cadena
  16. for my $cadena (@cadenas) {
  17.  
  18.     my $nueva_cadena = $cadena;         # combinación sin incertidumbre
  19.    
  20.     for my $i (0 .. 4) {                # resto de combinaciones de incertidumbres
  21.         substr(my $c = $cadena, $i, 1) = '.';
  22.         $nueva_cadena .= "|$c";
  23.     }
  24.    
  25.     $cadena = $nueva_cadena;            # guardamos en el array
  26. }
  27.  
  28. $cadenas = join q[|], @cadenas;
  29.  
  30. my $regex = qr/^($cadenas)/;            # aquí construiremos la expresión regular  
  31.  
  32. #say $regex;
  33.  
  34.  
  35. ## Lectura del fichero de bloques
  36. my ($titulo, $seq, $subtitulo, $texto);
  37.  
  38. open my $fh, q[<], 'kk.txt';
  39.  
  40. while (my $linea = <$fh>) {
  41.  
  42.     die if $linea !~ /^@/;      # Oops, algo fue mal
  43.  
  44.     $titulo    = $linea;        # leemos el bloque
  45.     $seq       = <$fh>;
  46.     $subtitulo = <$fh>;
  47.     $texto     = <$fh>;
  48.                                 # quitar finales de líneas
  49.     chomp($titulo, $seq, $subtitulo, $texto);
  50.  
  51. #    say "[$titulo, $seq, $subtitulo, $texto]";
  52.  
  53.     if ($seq =~ $regex) {       # buscamos en $seq las cadenas
  54. #        say "Encontrado $1 en $seq";
  55.  
  56.         say substr $seq,   5;           # Recorte
  57.         say substr $texto, 5;
  58.     }
  59. }
  60.  
  61. close   $fh;
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4
La salida es
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
TCGGGTTTTATGCC
+++++abtttbcbm
TCGGGTCCTATTCC
++b++abtntbccm
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


La expresión regular que formamos a partir de las cadenas de búsqueda tiene un aspecto como este:
(?-xism:^(CAGTT|.AGTT|C.GTT|CA.TT|CAG.T|CAGT.|ACGTA|.CGTA|A.GTA|AC.TA|ACG.A|ACGT.))

que, como ves, es una alternancia de todas las combinaciones que queremos probar.

De esta manera, nos ahorramos bucles y pruebas, dando un resultado muy rápido.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Localizar cadenas de ADN con comodines

Notapor danusol » 2010-12-09 02:54 @162

Muchísimas gracias Exlorer, y perdón por este parón en la discusión.

El resultado sí que tiene que incluir las cuatro líneas de cada bloque, no sólo la línea 2 y 4 (las que has llamado $título y $subtítulo), ya que el bloque completo está representado por las cuatro juntas.

La solución al problema de las indeterminaciones es muy buena; sí quiero todos los resultados a la vez, pero haciendo la búsqueda de todas las cadenas a la vez sería complicado que me separase cada cadena a un archivo.

Me explico. Si busco las cadenas CAGTT y ACGTA que he metido como argumento, quiero que me guarde
- por un lado los bloques de cuatro líneas correspondientes a CAGTT, .AGTT, C.GTT, CA.TT, CAG.T, y CAGT. en los que se han eliminado de la línea 2 la cadena y de la cuatro 4 los primeros caracteres correspondientes a la longitud de la cadena (cinco en este caso) en un archivo con nombre CAGTT_out para poder saber de dónde viene.
- por otro lado lo mismo para los correspondientes a ACGTA, .CGTA, A.GTA, AC.TA, ACG.A, y ACGT., en otro archivo distinto con nombre ACGTA_out.

Siempre tengo el mismo problema cuando intento leer archivos en los que tengo que encontrar una referencia y si existe, tengo que volver unas líneas atrás y coger lo anterior y lo siguiente a la referencia que busco. Si no estoy equivocado, $_ te guarda en memoria la linea que está leyendo pero ¿cómo puedo usar la línea anterior y las dos siguientes? ¿Creando una variable para cada línea? En este caso en el que estoy trabajando sé que cada bloque empieza por '@' y que los bloques son de cuatro líneas, pero ¿y si un bloque tiene longitud variable?

Gracias por la ayuda y los comentarios.
danusol
Perlero nuevo
Perlero nuevo
 
Mensajes: 46
Registrado: 2010-04-22 07:08 @339

Re: Localizar cadenas de ADN con comodines

Notapor explorer » 2010-12-09 13:33 @606

Bueno, entonces hay que hacer un doble bucle: uno por cada entrada del usuario, y otro interno, que abra el fichero, busque los patrones en él y cree los ficheros de salida.

En cuanto a la forma de leer registros, si sabes cuál es su separador, aunque tengan distinto tamaño, puedes usar la variable especial $/. Si estás seguro que la '@' no aparece en ningún otro sitio más que en el comienzo de cada registro, puedes hacer la lectura, en lugar de línea a línea, de bloque en bloque.

Más información en perldoc perlvar
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España


Volver a Bioinformática

¿Quién está conectado?

Usuarios navegando por este Foro: No hay usuarios registrados visitando el Foro y 3 invitados