Foro - Perl en Español

por **Giskard** » 2014-02-24 05:43 @279

Muy buenas.

He empezado con Perl hace relativamente poco, y estoy tratando de hacer el primer script de algo que sea útil y que no sea solo un ejemplo chorra.

Lo que intento conseguir con esto, es que este programa ejecute él solo varios BLAST a partir de un archivo de secuencias en FASTA, compruebe que resultados son comunes entre varias secuencias y te digas cuáles son, número de acceso, descripción y con qué longitud aparean cada una de las secuencias.

Para ello he usado BioPerl y he canibalizado un par de programas que nos pasó un profesor que nos estaba explicando esto (de hecho aún están los comentarios originales por ahí xD) para que lea las secuencias de FASTA y lance los BLAST, ya que no me llevo muy bien aún con BioPerl y con los objetos; y luego la parte de las comparaciones ya las he hecho yo.

Pues bien, cuando lo ejecuto me da el error:

Not a GLOB reference at comparacion.pl line 37.

He buscado por qué puede ser y no acabo de entender lo de los GLOB, y como es una parte de las que no he hecho yo, puede que, o bien la esté cagando y me falte algo, o me haya colado o lo que sea, por lo que si me podéis echar una mano en cómo solucionarlo os lo agradecería infinitamente. Aquí va el código, no seáis muy duros con él

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

use strict;
use Bio::Tools::Run::RemoteBlast;
 
#  use strict;
my $ficherosalida = "loquesea.txt";
my $prog          = 'blastp';
my $db            = 'swissprot';
my $e_val         = '1e-10';
my @secuencia_name;
my @secuencia_length;
my @secuencia_accesion;
my @secuencia_description;
my @params = (
    '-prog'       => $prog,
    '-data'       => $db,
    '-expect'     => $e_val,
    '-readmethod' => 'SearchIO'
);
 
my $factory = Bio::Tools::Run::RemoteBlast->new(@params);
 
#change a paramter
#  $Bio::Tools::Run::RemoteBlast::HEADER{'ENTREZ_QUERY'} = 'Homo sapiens [ORGN]';
 
#remove a parameter
# delete $Bio::Tools::Run::RemoteBlast::HEADER{'FILTER'};
 
my $v = 1;
my %hits_name;
my %hits_length;
my %hits_accesion;
my %hits_description;
 
#$v is just to turn on and off the messages
 
my $str = Bio::SeqIO->new( -file => 'sequences.fasta', '-format' => 'fasta' );
my @seq_array_name = <$str>;
for ( my $i = 0; $i < int(@seq_array_name); $i++ ) {
    my @array_name;
    my @array_length;
    my @array_accesion;
    my @array_description;
 
    while ( my $input = $str->next_seq() ) {
 
        #Blast a sequence against a database:
 
        #Alternatively, you could  pass in a file with many
        #sequences rather than loop through sequence one at a time
        #Remove the loop starting 'while (my $input = $str->next_seq())'
        #and swap the two lines below for an example of that.
        my $r = $factory->submit_blast($input);
 
        #my $r = $factory->submit_blast('amino.fa');
 
        print STDERR "waiting..." if ( $v > 0 );
        while ( my @rids = $factory->each_rid ) {
            foreach my $rid (@rids) {
                my $rc = $factory->retrieve_blast($rid);
                if ( !ref($rc) ) {
                    if ( $rc < 0 ) {
                        $factory->remove_rid($rid);
                    }
                    print STDERR "." if ( $v > 0 );
                    sleep 5;
                }
                else {
                    my $result = $rc->next_result();
 
                    #save the output
                    my $filename = $result->query_name() . "\.out";
                    $factory->save_output($filename);
                    $factory->remove_rid($rid);
                    $secuencia_name[$i]        = $result->query_name();
                    $secuencia_accesion[$i]    = $result->query_accesion();
                    $secuencia_length[$i]      = $result->query_length;
                    $secuencia_description[$i] = $result->query_description();
                    while ( my $hit = $result->next_hit ) {
                        next unless ( $v > 0 );
                        push( @array_name,        $hit->name );
                        push( @array_length,      $hit->length );
                        push( @array_accesion,    $hit->accesion );
                        push( @array_description, $hit->description );
                        print "\nhit name is ", $hit->name, "\n";
                    }
                }
            }
        }
 
    }
    $hits_name{$i}        = join( ";", @array_name );
    $hits_length{$i}      = join( ";", @array_length );
    $hits_accesion{$i}    = join( ";", @array_accesion );
    $hits_description{$i} = join( ";", @array_description );
}
 
my $tamano = scalar( keys %hits_name );
my @final;
for ( my $i = 0; $i < $tamano; $i++ ) {
    my @comparar = [];
    my @resultados1 = split /;/ => $hits_name{$i};
    for ( my $j = $i++; $j = $tamano; $j++ ) {
        my @resultados2 = split /;/ => $hits_name{$j};
        @comparar = comparar( $i, $j, @resultados1, @resultados2 );
 
        if ( $j == $tamano - 1 ) {
            my @igualados = igualdades(@comparar);
            @final = ( @final, redaccion(@igualados) );
        }
 
    }
}
guardar_datos( $ficherosalida, @final );
 
#################################
###########FUNCIONES#############
#################################
 
#Compara si los hits son el mismo y guarda la posición de los mismos dentro de la línea de datos del hash de guardado
sub comparar {
    my ( $cadena1, $cadena2, @result1, @result2 ) = shift;
    my @comparacion = [];
    my $tamres1     = int @result1;
    my $tamres2     = int @result2;
    for ( my $i = 0; $i < $tamres1; $i++ ) {
        for ( my $j = 0; $j < $tamres2; $j++ ) {
            if ( $result1[$i] eq $result2[$j] ) {
                my $a = int @comparacion;
                $comparacion[$a] = "$cadena1" . ";" . "$i" . ";" . "$cadena2" . ";" . "$j";
            }
 
        }
    }
    return @comparacion;
}
 
#Ve si hay más de un hit igual para algún gen y los ordena en un único resultado
sub igualdades {
    my (@resultados) = shift;
    my $tam = ( int @resultados ) - 1;
 
    for ( my $i = 0; $i < $tam; $i++ ) {
        my @res1 = split /;/ => $resultados[$i];
 
        for ( my $j = $i + 1; $j < $tam; $j++ ) {
 
            if ( $resultados[$j] != 0 ) {
                my @res2 = split /;/ => $resultados[$j];
 
                if ( $res1[1] == $res2[1] ) {
                    $resultados[$i] = $resultados[$i] . ";" . "$res2[2]" . ";" . "$res2[3]";
                    $resultados[$j] = 0;
                }
            }
        }
    }
    return @resultados;
}
 
#Con los datos de los hits comunes redacta una salida de datos para cada hit
sub redaccion {
    my @datos     = shift;
    my @redaccion = [];
    my $h         = 0;
    foreach my $resultado (@datos) {
        my @informacion = split /;/ => $resultado;
        my $tam = int @informacion;
        my @queries;
        my $hit;
        my $acceso;
        my $descripcion;
        my @longitudes;
        my $j = 0;
 
        for ( my $i = 0; $i < $tam; $i = $i + 2 ) {
            my $secuencia  = $informacion[$i];
            my $tama       = $informacion[ $i + 1 ];
            my $longitud   = $hits_length{$secuencia};
            my @todas_long = split /;/ => $longitud;
            $queries[$j]    = $secuencia_description[$secuencia];
            $longitudes[$j] = $todas_long[$tama];
            $j++;
        }
        my $sec1            = $informacion[0];
        my $pos1            = $informacion[1];
        my $hit_comun       = $hits_name{$sec1};
        my $hit_acceso      = $hits_accesion{$sec1};
        my $hit_descripcion = $hits_description{$sec1};
        my @hitcom          = split /;/ => $hit_comun;
        my @hitacc          = split /;/ => $hit_acceso;
        my @hitdes          = split /;/ => $hit_descripcion;
        my $nombrehit       = $hitcom[$pos1];
        my $nombreacc       = $hitacc[$pos1];
        my $nombredes       = $hitdes[$pos1];
        my $numerosecu      = int @queries;
        my $secredactado    = "";
        my $lonredactado    = "";
 
        for ( my $k = 0; $k < $numerosecu; $k++ ) {
            $secredactado = $secredactado . " y " . $queries[$k];
            $lonredactado = $lonredactado . " y " . $longitudes[$k];
        }
        $redaccion[$h]
            = "Las secuencias $secredactado tienen en comun el gen $nombrehit\n"
            . "cuya descripcion es $nombredes , su numero de acceso $nombreacc \n;"
            . "y cada uno de ellos aparea las siguientes longitudes respectivamente: $lonredactado ";
        $h++;
    }
    return @redaccion;
}
 
#Guarda los datos en un fichero de salida
sub guardar_datos {
    my ( $b, @salida ) = shift;
    ( open Fichero, ">$b" );
    my $N = int @salida;
    for ( my $i = 0; $i < $N; $i++ ) {
        print Fichero $salida[$i] . "\n";
    }
    close Fichero;
}
Coloreado en 0.013 segundos,  usando GeSHi 1.0.8.4

Cierto, no contaba con que al ponerlo aquí se descuadrase con respecto a mi ordenador: sería la línea 37.

Por lo que entiendo, lo que hace es un bucle para que se repita tantas veces como secuencias haya leído. Por ello le he puesto un bucle como el while de la línea 44 para recoger "a mano" las secuencias que hay, y luego meter ese número en el for, aunque no sé si estaría bien.

por **explorer** » 2014-02-24 11:48 @533

El error indica que en la línea 37 estás intentando usar $str como un gestor de archivo para leer el contenido del fasta, pero... no lo es.

Según la documentación, para obtener un gestor de archivo a través de Bio::SeqIO, se debe usar el método newFh(), con lo que la línea 36 debería quedar así:

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

my $str = Bio::SeqIO->newFh( -file => 'sequences.fasta', '-format' => 'fasta' );
Coloreado en 0.001 segundos,  usando GeSHi 1.0.8.4

por **Giskard** » 2014-02-25 05:36 @275

Pero por lo que puedo ver, new() y newFh() no es lo mismo, ¿no? De hecho, por lo que he leído aqui el objeto que crean los dos no es el mismo. De hecho, habiéndolo cambiado, ahora el error me da en la línea 44, de que no se puede usar next_seq en esa línea:

Can't locate object method "next_seq" via package "IO::File" at comparacionx.pl line 41, <GEN0> line 2.

Y cuando he probado ese comando habiendo usado new no me daba el problema.

por **explorer** » 2014-02-25 07:06 @337

Vale... se puede arreglar eso, pero... yo veo un problema.

En la línea 37 se leen todas las secuencias, ¿no? Entonces, ¿qué es lo que se pretende en la línea 44? Leer la siguiente secuencia, ¿de dónde? Si suponemos que las hemos leído todas, entonces están almacenadas en @seq_array_name. Y no se puede leer más.

¿O estoy equivocado?

¿Cual es el propósito del doble bucle de las líneas 38 y 44?

En cuanto a lo de los errores...

Para poder ejecutar la línea 37 -¡ojo!, escrita de esa manera-, es necesario el uso de newFh(). En cambio, para poder usar la línea 44, es necesario usar new().

Se puede solventar con el método fh():

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

my $str  = Bio::SeqIO->new(
            -file    => 'sequences.fasta',
           '-format' => 'fasta',
);                                           # Crear un flujo de objetos de secuencias FASTA
my $fh   = $str->fh();                       # Obtención del gestor de flujo asociado a $str
...;
my $seq  = <$fh>;                            # Forma de acceder al flujo, en forma de lectura de archivo, con $fh
...;
my $next = $str->next_seq();                 # Forma de acceder al flujo, en forma orientado a objetos, con $str
Coloreado en 0.001 segundos,  usando GeSHi 1.0.8.4

Yo creo que sí se pueden mezclar los dos métodos de lectura, pero sigo sin ver claro que queden secuencias posibles por leer, si se leen todas de golpe y se meten en un array.

por **Giskard** » 2014-02-26 06:33 @315

Joder, muchísimas gracias.

Lo del bucle cuando lo puse ya me pareció raro, pero era de esas partes que comentaba que estaban en otros programas que he reciclado y no lo toqué, aunque claro, sí que debería haberlo quitado. Lo he cambiado como me has dicho, tal que así:

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

  my $str  = Bio::SeqIO->new(
            -file    => 'sequences.fasta',
           '-format' => 'fasta',
  );                                           # Crear un flujo de objetos de secuencias FASTA
  my $fh   = $str->fh();                       # Obtención del gestor de flujo asociado a $str
  my $seq  = <$fh>;                            # Forma de acceder al flujo, en forma de lectura de archivo, con $fh
  my $numero_secuencias = 1;
  my @array_name;
  my @array_length;
  my @array_accession;
  my @array_description;
  my $i = 0;
    while (my $input = $str->next_seq()){
      #Blast a sequence against a database:
  
      #Alternatively, you could  pass in a file with many
      #sequences rather than loop through sequence one at a time
      #Remove the loop starting 'while (my $input = $str->next_seq())'
      #and swap the two lines below for an example of that.
      my $r = $factory->submit_blast($input);
      #my $r = $factory->submit_blast('amino.fa');
Coloreado en 0.001 segundos,  usando GeSHi 1.0.8.4

Y toda esa parte funciona. Ahora me estoy peleando con la parte de abajo, pero eso ya es más fácil xD

En serio, me has salvado, porque con esto voy a poder acortar gran parte de mi trabajo fin de máster una barbaridad, que si no me tocaba hacer 1000 BLAST o así, a mano, y mirar yo los resultados de uno en uno xD (sí, mi tutor es un poco bruto...).

por **Giskard** » 2014-02-27 04:54 @245

¡Puffff!, siento ser tan plasta, de verdad.

Me he puesto a comprobar los resultados del BLAST con dos secuencias, porque cada vez que hago alguna prueba se me eterniza esto, así que quería copiar y pegar los resultados, para luego modificar las funciones y cosas de más abajo sin tener que lanzar el BLAST. El caso es que me he dado cuenta de que no lanza, o al menos no obtiene los resultados de las dos secuencias, sino solo de la segunda, la cual además no se queda en la posición 2 (bueno, 1 si empezamos a contar de 0) sino que se queda el 1º (vamos, en la posición 0). Vuelvo a poner el código como lo tengo después de las modificaciones y lo que saca al final:

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

use strict;
use Bio::Tools::Run::RemoteBlast;
#  use strict;
  my $ficherosalida = "loquesea.txt";
  my $prog = 'blastp';
  my $db   = 'swissprot';
  my $e_val= '1e-10';
  my @secuencia_name;
  my @secuencia_length;
  my @secuencia_accession;
  my @secuencia_description;
  my @params = ( '-prog' => $prog,
         '-data' => $db,
         '-expect' => $e_val,
         '-readmethod' => 'SearchIO' );
 
  my $factory = Bio::Tools::Run::RemoteBlast->new(@params);
 
  #change a paramter
#  $Bio::Tools::Run::RemoteBlast::HEADER{'ENTREZ_QUERY'} = 'Homo sapiens [ORGN]';
 
  #remove a parameter
 # delete $Bio::Tools::Run::RemoteBlast::HEADER{'FILTER'};
 
  my $v = 1;
  my %hits_name;
  my %hits_length;
  my %hits_accession;
  my %hits_description;
  #$v is just to turn on and off the messages
 
  my $str  = Bio::SeqIO->new(
            -file    => 'sequences.fasta',
           '-format' => 'fasta',
  );                                           # Crear un flujo de objetos de secuencias FASTA
  my $fh   = $str->fh();                       # Obtención del gestor de flujo asociado a $str
  my $seq  = <$fh>;                            # Forma de acceder al flujo, en forma de lectura de archivo, con $fh
  my $numero_secuencias = 1;
  my @array_name;
  my @array_length;
  my @array_accession;
  my @array_description;
  my $i = 0;
    while (my $input = $str->next_seq()){
      #Blast a sequence against a database:
  
      #Alternatively, you could  pass in a file with many
      #sequences rather than loop through sequence one at a time
      #Remove the loop starting 'while (my $input = $str->next_seq())'
      #and swap the two lines below for an example of that.
      my $r = $factory->submit_blast($input);
      #my $r = $factory->submit_blast('amino.fa');
  
      print STDERR "waiting..." if( $v > 0 );
      while ( my @rids = $factory->each_rid ) {
        foreach my $rid ( @rids ) {
          my $rc = $factory->retrieve_blast($rid);
          if( !ref($rc) ) {
            if( $rc < 0 ) {
              $factory->remove_rid($rid);
            }
            print STDERR "." if ( $v > 0 );
            sleep 5;
          } else {
            my $result = $rc->next_result();
            #save the output
            my $filename = $result->query_name()."\.out";
            $factory->save_output($filename);
            $factory->remove_rid($rid);
            $secuencia_name[$i] = $result->query_name();
            $secuencia_accession[$i] = $result->query_accession();
            $secuencia_length[$i] = $result->query_length;
            $secuencia_description[$i] = $result->query_description();
            while ( my $hit = $result->next_hit ) {
              next unless ( $v > 0);
              push(@array_name,$hit->name);
              push(@array_length,$hit->length);
              push(@array_accession,$hit->accession);
              push(@array_description,$hit->description);
              print "\nhit name is ", $hit->name, "\n";
            }
          }
        }
        
      }
    $hits_name{$i}= join(";",@array_name);
    print "para la clave $i da $hits_name{$i}";
    $hits_length{$i}= join(";",@array_length);
    $hits_accession{$i}= join(";",@array_accession);
    $hits_description{$i}= join(";",@array_description);
    $numero_secuencias = $numero_secuencias + 1;
    $i++;
    }
    
  
  
  
  my @final;
  print "los hits name 0 son", $hits_name{0},"\n";
  print "los hits name 1 son", $hits_name{1}, "\n";
  print "los hits length 0 son", $hits_length{0}, "\n";
  print "los hits length 1 son", $hits_length{1}, "\n";
  print "los hits accession 0 son", $hits_accession{0}, "\n";
  print "los hits accession 1 son", $hits_accession{1}, "\n";
  print "los hits description 0 son", $hits_description{0}, "\n";
  print "los hits description 1 son", $hits_description{1}, "\n";
  print "los query name 0 son", $secuencia_name[0], "\n";
  print "los query name 1 son", $secuencia_name[1], "\n";
Coloreado en 0.003 segundos,  usando GeSHi 1.0.8.4

Foro - Perl en Español

Problemas al leer FASTA y hacer BLAST

Problemas al leer FASTA y hacer BLAST

Publicidad

Re: Problemas al leer FASTA y hacer BLAST

Re: Problemas al leer FASTA y hacer BLAST

Re: Problemas al leer FASTA y hacer BLAST

Re: Problemas al leer FASTA y hacer BLAST

Re: Problemas al leer FASTA y hacer BLAST

Re: Problemas al leer FASTA y hacer BLAST

Re: Problemas al leer FASTA y hacer BLAST

¿Quién está conectado?