• Publicidad

Cómo manejar cromosomas y posiciones

Perl aplicado a la bioinformática

Cómo manejar cromosomas y posiciones

Notapor danusol » 2011-06-27 02:31 @146

¡Hola perleros!
Tengo dos tablas de datos que quiero cruzar. En la más pequeña (llamémosla "pequeña") tengo posiciones cromosómicas en esta forma:

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
chr1    1356442
chr1    1740013
chr1    21332987
.
.
.
chr2    34298384
chr2    34299811
chr2    70007885
.
.
.
chr3    12335789
.
.
.
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Así para cada cromosoma, cientos de entradas.

Y quiero cruzar esta tabla con otra principal (llamémosla "general") enorme con casi todas las posiciones en el genoma que tiene tres columnas (cromosoma, coordenada y un valor numérico), ordenado por cromosoma y posición.

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
chr1    1       10
chr1    2       11
chr1    3       11
chr1    4       12
.
.
.
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Lo que quiero hacer es buscar las posiciones de "pequeña" en "general" y sacar a otro archivo la posición y el valor numérico para dicha posición. Es decir quiero filtrar el archivo "general" en función de las posiciones en "pequeña".

¿Cómo puedo manejar estos datos? Me imagino que el archivo "general" no lo puedo cargar en memoria por lo que éste lo tendría que abrir línea a linea. Pero entonces no sé cómo hacer la búsqueda de cada entrada de "pequeña". Como los archivos están ordenados por cromosoma, si quiero buscar el valor numérico para, por ejemplo, chr18 12224555 tendría que leer prácticamente todo "general" hasta llegar a la sección de chr18.

Espero haberlo explicado de manera entendible. Gracias por vuestra ayuda.

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

Publicidad

Re: Cómo manejar cromosomas y posiciones

Notapor explorer » 2011-06-27 04:22 @223

Yo creo que sí puedes leer "general" a memoria y convertirlo en un hash, salvo que me digas que tiene varios gigabytes de tamaño o que debes realizar esta operación varias veces al día, y tarda mucho tiempo en cada ocasión, para leerlo. En esos casos, hay que pensar algo distinto. Mientras tanto, sigue una solución normal.

El hash sería, poniendo de valor de la clave, el valor de la primera columna y el valor de la segunda columna de general. Y como valor del hash Algo así:

$cromosoma{"$primera,$segunda"} = $tercera;

Luego, solo te queda abrir "pequeña", línea a línea, y sacar el valor que buscas:

say $cromosoma{"$primera,$segunda"};
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: Cómo manejar cromosomas y posiciones

Notapor danusol » 2011-06-27 05:53 @286

Gracias explorer. Pues el "general" ocupa 700Mb, pero tampoco se va a ejecutar muchas veces el programa, la verdad.

Voy a por ello.

Gracias,

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

Re: Cómo manejar cromosomas y posiciones

Notapor danusol » 2011-06-27 10:09 @464

hola de nuevo explorer. He intentado hacer algo como lo que me has sugerido, pero creo que tengo algún error de concepto y lío en el uso de las variables, porque al ejecutar el programa me da errores de declaración de variables:

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
$ perl coberturaSelectedPositions.pl general.txt peque.txt out.txt
Global symbol "$cromosoma" requires explicit package name at coberturaSelectedPositions.pl line 34.
Global symbol "$cromosoma" requires explicit package name at coberturaSelectedPositions.pl line 36.
Global symbol "$cromosoma" requires explicit package name at coberturaSelectedPositions.pl line 45.
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


aquí debajo puedes ver mi código. Pienso que al declarar la variable $cromosoma en el ámbito general, podría usarla dentro del bucle, pero parece que no la he declarado bien?

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl -w
  2. use common::sense;
  3. use Data::Dumper::Simple;
  4.  
  5. ### Leemos el archivo grande, contiene posiciones y valor
  6. my %cromosomas;
  7.  
  8. open (IN1, $ARGV[0]) || die "File not found\n"; #archivo "general" contiene todos los datos
  9.  
  10. while (<IN1>) {
  11.     chomp;
  12.  
  13.     my ($cromosoma, $posicion, $valor) = split;
  14.  
  15.     push @{ $cromosomas{ $cromosoma } }, [ $posicion , $valor];   # estructura de hash de arrays de arrays
  16.                                                                 # por cada $cromosoma, creamos un array donde guardamos,
  17.                                                                 # en cada elemento, otro array, con dos elementos:
  18.                                                                 # la $posición y el $valor
  19. }
  20.  
  21. close (IN1);
  22.  
  23. say Dumper %cromosomas;                                        
  24.  
  25. ### Leer archivo de interés, contiene las posiciones de interés
  26. open (IN2, $ARGV[1]) || die "File not found\n";
  27. open OUT, ">", $ARGV[2] or die $!; ##crear un archivo de resultados, donde se guardara el gene_symbol y su Id
  28.  
  29. while (<IN2>) {
  30.     chomp;
  31.    
  32.     my ($chr, $coord) = split; #declaro mis columnas en el nuevo archivo
  33.  
  34.     if ( $cromosomas{ $cromosoma } eq $chr ) {                           # si el cromosoma lo tenemos
  35.    
  36.         my @array = @{ $cromosomas{ $cromosoma } };              # sacamos todas sus posiciones y valores
  37.  
  38.         for my $array_ref (@array) {                             # por cada una de estas posiciones y valores
  39.  
  40.             my($posicion, $valor) = @$array_ref;                   # sacamos la posición y el valor
  41.  
  42.             if ($posicion == $coord) {    #y comparo coordenadas en ambos archivos    
  43.                                                                                        
  44.  
  45.             print OUT "$cromosoma\t$coord\t$valor";                    # escribir resultado a archivo salida
  46.  
  47.             last;                               # saltar a la siguiente posición
  48.           }
  49.       }
  50.   }
  51.  
  52. }
  53.  
  54. close (IN2);
  55. close (OUT);
  56.  
  57.  
Coloreado en 0.005 segundos, usando GeSHi 1.0.8.4
Última edición por danusol el 2011-06-28 01:57 @123, editado 1 vez en total
danusol
Perlero nuevo
Perlero nuevo
 
Mensajes: 46
Registrado: 2010-04-22 07:08 @339

Re: Cómo manejar cromosomas y posiciones

Notapor explorer » 2011-06-27 10:12 @466

Perl dice que no sabe de dónde viene $cromosoma. Es decir, que no ha sido declarada antes de ser usada por primera vez.

Debes declararla con un my() antes de usarla.

Y, sobre todo, darle una valor cada vez.

Por otra parte, me parece que la línea 34 está mal... hay un operador de asignación, dentro de un if().
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: Cómo manejar cromosomas y posiciones

Notapor danusol » 2011-06-28 02:03 @127

Gracias explorer, he modificado el if() de la línea 34 para que la comparación sea de strings, pero sigo sin entender por qué tengo que declarar $cromosoma, si la he declarado globálmente en la línea 13, ¿o no?. Quiero usar la variable $cromosoma que viene del archivo "general" (ARGV[0]) y compararla con la variable $chr declarada dentro del while (<IN2>). Si vuelvo a declararla con un my() dentro del while, ¿no se comportaría como una variable local?
danusol
Perlero nuevo
Perlero nuevo
 
Mensajes: 46
Registrado: 2010-04-22 07:08 @339

Re: Cómo manejar cromosomas y posiciones

Notapor explorer » 2011-06-28 04:06 @212

En la línea 13 estás declarando (y definiendo) la variable $cromosoma como local dentro del primer while() (fíjate en las llaves, que son las que delimitan el contexto en que es conocida la variable).

En cambio, en la línea 34, el programa se encuentra con una variable desconocida, $cromosoma (porque era local en el while anterior, no global), y además, no tiene ningún valor (no se lo has puesto).

Prueba a poner esta línea:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. my $cromosoma = $chr;
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


que, aunque no es lo que te he respondido en mi respuesta, es así como estás construyéndolo en la primera estructura.
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: Cómo manejar cromosomas y posiciones

Notapor explorer » 2011-06-28 04:18 @221

Esta es mi versión:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use common::sense;
  3. use autodie;
  4.  
  5. # use Data::Dumper::Simple;
  6.  
  7. @ARGV == 3 or die "Uso: $0 <fichero general> <fichero> <fichero salida>\n";
  8.  
  9. my %cromosomas;
  10.  
  11. ### Leemos el archivo grande, contiene posiciones y valor
  12. open my $IN1, "<", $ARGV[0];
  13.  
  14. while (<$IN1>) {
  15.     chomp;
  16.  
  17.     my ($cromosoma, $posicion, $valor) = split;
  18.  
  19.     $cromosomas{ "$cromosoma,$posicion" } = $valor;     # un hash
  20. }
  21.  
  22. close $IN1;
  23.  
  24. #say Dumper %cromosomas;                                        
  25.  
  26. ### Leer archivo de interés, contiene las posiciones de interés
  27. open my $IN2, "<", $ARGV[1];
  28. open my $OUT, ">", $ARGV[2];    # archivo de resultados, donde se guardara el gene_symbol y su Id
  29.  
  30. while (<$IN2>) {
  31.     chomp;
  32.  
  33.     my ($cromosoma, $posicion) = split;
  34.  
  35.     my $clave = "$cromosoma,$posicion";
  36.     my $valor = $cromosomas{ $clave };
  37.  
  38.     if ( $valor ) {
  39.         say $OUT join "\t", $cromosoma, $posicion, $valor;
  40.     }
  41. }
  42.  
  43. close $IN2;
  44. close $OUT;
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4
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


Volver a Bioinformática

¿Quién está conectado?

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