• Publicidad

Script que funcione como trimmer de secuencias de nts

Perl aplicado a la bioinformática

Script que funcione como trimmer de secuencias de nts

Notapor 3mgcantarero » 2017-11-12 04:52 @245

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
3mgcantarero
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2017-11-12 04:41 @237

Publicidad

Re: Script que funcione como trimmer de secuencias de nts

Notapor explorer » 2017-11-12 14:21 @639

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
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14480
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