Foro - Perl en Español

Publicado: **2011-09-06 03:44 @197**

Hola,

necesito ayuda con una tarea que en principio parece sencilla, pero que dado el tamaño de los ficheros de datos que se manejan se me está volviendo compleja.

El caso es que tengo dos ficheros, uno de ellos con 22.500 líneas y en cada línea varias columnas de las que me interesan las dos primeras (cromosoma y posición). En otro fichero, de 180.000 líneas, tengo 3 columnas (cromosoma, inicio y fin). La cosa es: necesito por cada posición del primer fichero ver si está o no en alguno de los rangos del segundo.

He conseguido una solución "a lo bestia" que consiste en cargar los dos ficheros en sendos arrays y recorrer el primero. Por cada posición del primero me recorro el segundo para ver si está en algún rango (con alguna medida para no recorrerme el fichero entero, claro).

Ejemplo fichero de valores:

Sintáxis: [ Descargar ] [ Ocultar ]

Using text Syntax Highlighting

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  
chr10   76210   .       A       G       88.8    .       DP=2;AF1=1;CI95=0.5,1;DP4=0,0,1,1;MQ=60;FQ=-33  GT:PL:GQ        1/1:120,6,0:10
chr10   94238   .       A       T       28      .       DP=3;AF1=1;CI95=0.5,1;DP4=0,0,0,3;MQ=45;FQ=-36  GT:PL:GQ        1/1:60,9,0:16
chr10   131577  .       A       T       4.77    .       DP=14;AF1=0.4999;CI95=0.5,0.5;DP4=6,6,1,1;MQ=60;FQ=6.99;PV4=1,4.9e-07,1,0.053   GT:PL:GQ        0/1:33,0,255:33
chr10   146528  .       G       A       19.8    .       DP=2;AF1=1;CI95=0.5,1;DP4=0,0,1,1;MQ=60;FQ=-33  GT:PL:GQ        1/1:51,6,0:10
 
Coloreado en 0.000 segundos,  usando GeSHi 1.0.8.4

Ejemplo fichero de rangos:

Sintáxis: [ Descargar ] [ Ocultar ]

Using text Syntax Highlighting

track:exome and plus captures           
chr1    29775   30931
chr1    68569   70529
chr1    227733  228854
chr1    227971  229211
chr1    367147  369108
chr1    470471  471830
chr1    620584  622545
chr1    740665  741785
 
Coloreado en 0.000 segundos,  usando GeSHi 1.0.8.4

Mi "solución":

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

#! /usr/bin/perl
 
my @ranges;
my @values;
 
open( RANGES, $ARGV[0] ) or die "Can't open file: $!";
while ( my $line = <RANGES> ) {
    chomp $line;
    my @val = split( "\t", $line );
    push( @ranges, $val[0] . "_" . $val[1] . "_" . $val[2] );
}
close(RANGES);
 
open( VALUES, $ARGV[1] ) or die "Can't open the file $!";
while ( my $line = <VALUES> ) {
    chomp $line;
    my @val = split( "\t", $line );
    if ( $val[0] !~ "#" ) { $leer = 1; }
    if ( $leer == 1 ) {
        push( @values, $val[0] . "_" . $val[1] );
    }
}
close(VALUES);
 
foreach $Val (@values) {
    my @aux = split( "_", $Val );
    foreach $range (@ranges) {
        my @aux2 = split( "_", $range );
        if ( $aux[1] >= $aux2[1] && $aux[1] <= $aux2[2] ) {
            chomp $aux2[2];
            print "$aux[1] esta en el rango [$aux2[1]-$aux2[2]] \n";
            last;
        }
        elsif ( $aux[1] < $aux2[1] ) {
            print "$aux[1] no esta en ningun rango\n";
            last;
        }
    }
}
 
Coloreado en 0.003 segundos,  usando GeSHi 1.0.8.4

Viendo el problema, ¿existe alguna solución más optimizada que la de recorrerse ambos arrays a piñón?

¡Muchas gracias!

Publicado: **2011-09-06 05:33 @273**

Bienvenido a los foros de Perl en español, Dixan_Megaperl.

En este hilo comentamos otra solución.

Primero, leemos el fichero con los cromosomas y sus posiciones. Guardamos la información en una estructura de hash de array de array (por cada cromosoma tenemos un array, que guarda arrays con la posición e ID.

Segundo, hacemos una reordenación numérica (de menor a mayor) de las posiciones de los cromosomas.

Tercero, abrimos el fichero con los rangos. Por cada línea, primero comprobamos si tenemos información de ese cromosoma, y si es así, hacemos un bucle por todas las posiciones de ese cromosoma. Como están ordenadas numéricamente, es fácil saltar por las posiciones que son menores al inicio del rango. Una vez que encontramos una posición mayor que el inicio del rango, comprobamos si está antes del final del rango, y si es así, pintamos la información de la coincidencia, y terminamos el bucle.

Y es casi lo mismo a lo que tienes hecho, salvo que tu no miras por la coincidencia del nombre del cromosoma, sino que estás mirando todas las posiciones para todos los rangos.

Tu solución quizás se pueda acelerar, si eliminas los split(), y cambias la estructura de datos, de una dimensión a dos (array de array) (habría que lanzar un Benchmark, para comprobarlo).

Publicado: **2011-09-06 07:40 @361**

¡Muchas gracias!

Con la información del enlace he conseguido mejorar el tiempo de ejecución del programa. He bajado de 6 minutos a 1 y poco.

Foro - Perl en Español

Búsqueda de valores en rangos

Búsqueda de valores en rangos

Re: Búsqueda de valores en rangos

Re: Búsqueda de valores en rangos