• Publicidad

Cruzar dos listas

Perl aplicado a la bioinformática

Cruzar dos listas

Notapor danusol » 2011-02-21 09:46 @449

Hola perleros,

Tengo dos archivos de datos. El primero me indica el cromosoma, el nombre del gen (gene_symbol) y la posición central del gen en el cromosoma. Tengo 2000 filas. El segundo me indica el cromosoma, la posición inicial y final del gen y el identificador (ID) del gen. Tengo 35000 filas. Lo que necesito es obtener en un archivo resultante, el gene_symbol y su ID propio.

Para ello he pensado leer el archivo menor y guardar gene_symbol cromosoma y la posición central del gen en un array cada uno; luego leer el segundo archivo y pasar cada miembro de los arrays e interrogar el segundo archivo con condiciones secuenciales para que si la primera no se cumple pase a buscar en la siguiente fila, del estilo
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #en el archivo pequeño
  2. open(my $fh, '<', "lista.txt") or die $!;
  3. my @cromosoma,@symbol,@pos;
  4. #leer el archivo por líneas y separar cada campo para introducirlo en el array respectivo
  5. while (<>) {
  6.  chomp;
  7.  $cromosoma,$pos,$symbol = split(/\t/);
  8. }
  9. close $fh;
  10.  
  11. # en la base de datos
  12. open(my $fh, '<', "BD.txt") or die $!;
  13. while (<>) {
  14.  chomp;
  15.  my ($Xsoma,$start,$end,$ID) = split(/\t/);
  16.  my $n=0; #contador
  17.  foreach (@pos){
  18.    n++;
  19.   if ($cromosoma eq $Xsoma[$n]){   #si coinciden los cromosomas
  20.      if ($pos > $start & $pos <$end){ #si la posición está entre el incio y el fin
  21.        my @result_symbol = push $symbol[$n]; #guardo el gene_symbol correspondiente
  22.        my @result_ID = push $ID; #guardo el ID correspondiente
  23.        #¿sería mejor ir imprimiendo los resultados en el archivo de salida?
  24.      }
  25.    }
  26.  }
  27. }
  28.  
  29. close $fh;
  30.  
  31. # abrir el archivo de salida y escribir los dos arrays en dos columnas,
  32. # o igual mejor, abrir el archivo de salida antes e ir escribiendo los resultados a medida que salen, porque ¿sería más eficiente?
Coloreado en 0.004 segundos, usando GeSHi 1.0.8.4


Pero me surge la duda de, aparte de si está bien codificado, si seguirá buscando el mismo $pos dentro de la base de datos hasta que la lea entera una vez encontrado un valor que cumple ambas condiciones o si lo encuentra una vez y ya pasa al siguiente $pos.

Gracias por vuestra ayuda,

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

Publicidad

Re: Cruzar dos listas

Notapor pvaldes » 2011-02-21 10:16 @469

if ($pos > $start & $pos <$end)

De entrada probablemente quieras incluir un >= aquí
pvaldes
Perlero nuevo
Perlero nuevo
 
Mensajes: 129
Registrado: 2011-01-22 12:56 @580

Re: Cruzar dos listas

Notapor danusol » 2011-02-21 10:33 @481

Gracias pvaldes. Ya estuve pensando en eso, pero no creo que sea necesario porque $pos ha de estar incluido en el intervalo ($start,$end). Se supone que, aunque no sea exacto, $pos es el centro del gen que ha de estar en ese intervalo.
danusol
Perlero nuevo
Perlero nuevo
 
Mensajes: 46
Registrado: 2010-04-22 07:08 @339

Re: Cruzar dos listas

Notapor explorer » 2011-02-21 10:44 @489

¿Podrías publicar un ejemplo breve de los ficheros de entrada y el fichero de salida?

¿Dónde se inicializa @pos?

Faltan unas cuantas líneas...
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Cruzar dos listas

Notapor danusol » 2011-02-21 12:06 @545

Hola Explorer,

aquí va un trocito de la lista
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
Chr1    19160151        4CL1
Chr4    13221122        ABI1
Chr5    25894043        ABR1
Chr1    4511450 ACA.l
Chr3    21213915        ACA11
Chr1    11411539        ACBP6
Chr5    6590160 ACL5
Chr2    8476733 ACO1
Chr4    13546819        ACO2
Chr1    23083211        ACO2
Chr2    1138898 ACR5
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


y un trocito de la base de datos
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
AT3G18710       6434083 6435561 Chr3
AT4G25880       13154708        13159382        Chr4
AT1G71695       26964249        26966688        Chr1
AT5G41480       16595927        16598636        Chr5
AT5G15008       4856975 4857139 Chr5
AT3G18310       6284356 6287144 Chr3
AT1G62380       23082170        23084253        Chr1
AT4G26970       13542918        13548629        Chr4
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


La base de datos es el archivo más grande que incluye a todos los de la lista. No sé si me coincidirán más pero el ejemplo de "lista" incluye dos genes que se llaman ACO2 (hacia el final) pero en distintos cromosomas, y el programa tiene que ver que uno coincide con el último de la base de datos porque su $pos está comprendido entre start-end de AT4G26970 y el otro con el penúltimo (AT1G62380) de la base de datos. Por tanto el resultado tendría que ser

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
ACO2 AT1G62380
ACO2 AT4G26970
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Por otro lado, pensaba que @pos estaba inicializado en la línea 3 del script, pero igual lo tendría que modificar a

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. my @cromosoma,@symbol,@pos = ();
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4
danusol
Perlero nuevo
Perlero nuevo
 
Mensajes: 46
Registrado: 2010-04-22 07:08 @339

Re: Cruzar dos listas

Notapor explorer » 2011-02-21 12:20 @556

He dicho inicializar... es decir: la línea 17 hace un bucle por los elementos de @pos, pero en ninguna parte se han guardado esos valores dentro de @pos.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Cruzar dos listas

Notapor danusol » 2011-02-21 12:33 @565

cierto, lo acabo de encontrar y justo lo iba a indicar,

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. while (<>) {
  2.    
  3.        chomp;
  4.    
  5.        my $cromosoma,$pos,$symbol = split(/\t/);
  6.            @pos = push $pos;
  7.            @cromosoma = push $cromosoma;
  8.            @symbol = push $symbol;
  9.    
  10.       }
  11.    
  12.       close $fh;
  13.  
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4
danusol
Perlero nuevo
Perlero nuevo
 
Mensajes: 46
Registrado: 2010-04-22 07:08 @339

Re: Cruzar dos listas

Notapor explorer » 2011-02-21 13:08 @589

¿ @pos = push $pos; ? ¿Esto te funciona?

Me temo que no has ejecutado el programa... :)
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Cruzar dos listas

Notapor danusol » 2011-02-21 13:35 @607

Desgraciadamente, lo he escrito antes de probarlo, y ahora estoy revisando varias cosas. Lo vuelvo a adjuntar un poco mejorado ¡pero se me queda colgado en el primer while! normal, está mal empleado el filehandle, lo modifico y funciona.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.  open(fh, '<', "lista.txt") or die $!;
  2.    
  3.       my (@cromosoma,@symbol,@pos) = ();
  4.                 my ($cromosoma,$pos,$symbol);
  5.      #leer el archivo por líneas y separar cada campo para introducirlo en el array respectivo
  6.    
  7.       while (<fh>) {
  8.  
  9.        chomp;
  10.        ($cromosoma,$pos,$symbol) = split('\t');
  11.            push (@pos,$pos);
  12.            push (@cromosoma,$cromosoma);
  13.            push (@symbol,$symbol);
  14.    print "hola\n";#imprime
  15.       }
  16.    
  17.       close (fh);
  18.  
  19.        
  20.  
  21.       # en la base de datos
  22.  
  23.       open(fh2, '<', "TrozoDB.txt") or die $!;
  24.  
  25.       while (<fh2>) {
  26.  
  27.        chomp;
  28.  
  29.        my ($ID,$start,$end,$Xsoma) = split('\t');
  30.  
  31.        my $n=0; #contador
  32.  
  33.        foreach (@pos){
  34.  
  35.            
  36.         if ($cromosoma[$n] eq $Xsoma){   #si coinciden los cromosomas
  37.  
  38.            if (($pos[$n] > $start) && ($pos <$end)){ #si la posición está entre el incio y el fin
  39.  
  40.              #my @result_symbol = push $symbol[$n]; #guardo el gene_symbol correspondiente
  41.                                 print $symbol[$n],"\n";
  42.  
  43.              #my @result_ID = push $ID; #guardo el ID correspondiente
  44.                                 print $ID,"\n";
  45.              #¿sería mejor ir imprimiendo los resultados en el archivo de salida?
  46.                                 $n++;
  47.            }
  48.  
  49.          }
  50.  
  51.        }
  52.  
  53.       }
  54.  
  55.        
  56.  
  57.       close (fh2);
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4
danusol
Perlero nuevo
Perlero nuevo
 
Mensajes: 46
Registrado: 2010-04-22 07:08 @339

Re: Cruzar dos listas

Notapor pvaldes » 2011-02-21 16:23 @724

danusol escribiste:> incluye dos genes que se llaman ACO2 (hacia el final) pero en distintos cromosomas


Los genes duplicados, sí, iba a preguntarte por ellos... tienes que definir previamente qué vas a hacer con eso... veo tres posibles enfoques:

1.- Simplemente quieres una lista única de todos los genes y sus ID y al mismo nombre corresponde la misma ID siempre, en cuyo caso tu sistema actual para separarlos no es óptimo porque corre muchos más ciclos de lo necesario, estás haciendo 70 millones de comparaciones para encontrar 2000 nombres y 2000 números.

Meter una comprobación previa debería ahorrar un considerable tiempo de cálculo.

if gen_ID YA existe en tu arreglo ... saltar al siguiente elemento en la comparación
else push gen_ID

2.- O bien quieres una lista de todos los genes por cromosoma, serán "diferentes" genes si están en diferentes cromosomas aunque tengan la misma información y tendrán ID distintas, en cuyo caso tienes que encontrar el modo de diferenciarlos, (y de entrada sería práctico denominarlos como cromosoma_gen-symbol, es decir aplicar un join() de los elementos 1+3 después de hacer el split() y trabajar sobre ello)

3.- O lo que buscas en realidad es sacar un hash con key = gen_symbol y como valor una lista de valores ID con un número indeterminado de elementos, que sería una variante de la segunda vía.
pvaldes
Perlero nuevo
Perlero nuevo
 
Mensajes: 129
Registrado: 2011-01-22 12:56 @580

Siguiente

Volver a Bioinformática

¿Quién está conectado?

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

cron