Página 1 de 1

Script que funcione como trimmer de secuencias de nts

NotaPublicado: 2017-11-12 04:52 @245
por 3mgcantarero
Hola.

Tengo que realizar un script que me permita leer un archivo fasta usando biotools (proporcionado por el usuario en la consola al invocar el programa) y que elimine las N del extremo 5' las del 3', que recorte los dos extremos en un número de nucleótidos indicados por el usuario en la consola al invocar el programa y que me devuelva un archivo con las secuencias recortadas e imprima por pantalla la cabecera del fasta más el tamaño de origen más el tamaño resultante más la diferencia.

He empezado con ello pero estoy atascadísima ya que al intentar eliminar las N con un bucle if(), parece que no me leyera el archivo.

¿Alguien me podría echar una mano?

Esto es lo que llevo hecho hasta ahora y el archivo de salida trimmedfile.fa sale vacío. Si quito el primer bucle if() después del while(), me imprime en el archivo de salida, pero no elimina ninguna N. ¿Alguien me puede ayudar?

Muchas gracias
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2.  
  3. # Introduzco los módulos de Perl que voy a necesitar
  4.  
  5. use strict;
  6. use warnings;
  7. use Getopt::Long;
  8. use Bio::SeqIO;
  9. open( OUT, ">trimmedfile.fa" );
  10.  
  11. # Declaro los indicadores que voy a obtener con GetOpt al invocar el programa
  12.  
  13. my $file;
  14. my $trim5;
  15. my $trim3;
  16. my $help;
  17.  
  18. # Asocio los parámetros proporcionados por la consola con las variables declaradas
  19.  
  20. GetOptions(
  21.     "file=s" => \$file,
  22.  
  23.     #   "trim5"=>\$trim5,
  24.     #   "trim3"=>\$trim3,
  25.     #   "help"=>\$help
  26. );
  27.  
  28. # Leo el fichero FASTA con BioPerl
  29.  
  30. if ($file) {
  31.     my $seqio = Bio::SeqIO->new( -format => "fasta", -file => "sequences.fa" );
  32.     while ( my $seq = $seqio->next_seq ) {
  33.  
  34.         # Busco líneas que comiencen con N
  35.         if ( $seq =~ /^N/ ) {
  36.             $seq =~ s/N//;                     # Pretendo quitar las N. Pendiente ver cómo itero cada línea hasta eliminar todas
  37.             print OUT $seq->display_id() . "\n";    # Imprimir la cabecera
  38.             print OUT $seq->seq . "\n";        # Imprimir la secuencia sin una N al principio
  39.         }
  40.     }
  41. }
  42. else {
  43.     # Como no sé ha incluido ningún parámetro o se han incluido parámetros incorrectos ponemos otra ayuda
  44.     print "\nWrong parameters!!!:\n";
  45.     print "\nHelp:\n\nperl $0 -file -trim5 -trim3\n";
  46.  
  47. }
  48. close OUT;
  49.  
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4

Re: Script que funcione como trimmer de secuencias de nts

NotaPublicado: 2017-11-12 14:21 @639
por explorer
Bienvenida a los foros de Perl en Español, 3mgcantarero.

Lo que devuelve el método next_seq() no es la secuencia de nucleótidos, sino un objeto Bio::Seq.

Para acceder a la secuencia de nucleótidos debes usar el método seq(), igual a como lo haces en la línea 38:

my $secuencia = $seq->seq();

Atención aquí: no sabemos si las letras llegan en mayúsculas o minúsculas (o mezclado), así que puedes hacer varias cosas:

1) pasar la secuencia a mayúsculas: $secuencia = uc $secuencia; y luego ya le puedes pasar la expresión regular

2) no hacer nada, y modificar la expresión regular para que no tenga en cuenta el tamaño de caja de las letras: agregas '/i': $secuencia =~ s/N//i;

Quedaría algo así (no probado):
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. if (-e $file) {
  2.     my $seqio = Bio::SeqIO->new( -format => "fasta", -file => $file );
  3.  
  4.     while ( my $seq = $seqio->next_seq ) {
  5.  
  6.         my $secuencia = $seq->seq();
  7.  
  8.         # Si la secuencia comienza con N
  9.         if ($secuencia =~ m/^N/i ) {
  10.  
  11.             $secuencia =~ s/N//gi;                      # Quitar todas las N
  12.  
  13.             print OUT $seq->display_id() . "\n";        # Imprimir la cabecera
  14.             print OUT $secuencia         . "\n";        # Imprimir la secuencia
  15.         }
  16.     }
  17. }
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4