Foro - Perl en Español

La comunidad de programadores en Perl de habla hispana.

Obviar


  • Publicidad

Filtrado genes en archivos FASTA

Perl aplicado a la bioinformática
Publicar una respuesta

Filtrado genes en archivos FASTA

Notapor jalapea » 2015-03-04 10:16 @469

Un buen día.

Estoy tratando de filtrar los genes de un cromosoma que se encuentran dentro del genoma en un archivo FASTA.

Me podrían dar luces de cómo iniciar el proceso.
jalapea
Perlero nuevo
Perlero nuevo
 
Mensajes: 7
Registrado: 2015-03-03 22:42 @987
Arriba

Publicidad

Re: Filtrado genes en archivos FASTA

Notapor explorer » 2015-03-04 14:39 @652

Bienvenido a los foros de Perl en Español, jalapea.

Si nos das alguna pista, mejor.

De momento, mira este foro de Bioinformática, porque hemos hablado bastantes veces de los archivos FASTA.

Usa el sistema de búsqueda. O danos esa pista.
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
Arriba

Filtrado de secuencias de un archivo FASTA

Notapor jalapea » 2015-03-05 11:40 @528

¿Me pueden colaborar en cómo hacer para realizar un filtrado de algunas secuencias?

Estoy tratando de realizar un script, que el filtrado lo realice teniendo en cuenta el siguiente parámetro: :GRCh38:7:, es decir, que solo me filtre esas secuencias y guardarlas en un nuevo archivo FASTA.

Les agradezco la colaboración.

He pensado en realizar un hash con arreglos para hacer ese filtrado.

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
>ENST00000415118 havana_ig_gene:known chromosome:GRCh38:14:22438547:22438554:1 gene:ENSG00000223997 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene
GAAATAGT
>ENST00000448914 havana_ig_gene:known chromosome:GRCh38:12:22449113:22449125:1 gene:ENSG00000228985 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene
ACTGGGGGATACG
>ENST00000431870 havana_ig_gene:known chromosome:GRCh38:7:105894508:105894523:-1 gene:ENSG00000227800 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
TGACTACGGTGACTAC
>ENST00000414852 havana_ig_gene:known chromosome:GRCh38:14:105913222:105913237:-1 gene:ENSG00000233655 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
TGACTACAGTAACTAC
>ENST00000390578 havana_ig_gene:known chromosome:GRCh38:7:105897957:105897987:-1 gene:ENSG00000211918 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
AGGATATTGTAGTGGTGGTAGCTGCTACTCC
>ENST00000390571 havana_ig_gene:known chromosome:GRCh38:14:105886031:105886061:-1 gene:ENSG00000211911 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
GTATTACTATGATAGTAGTGGTTATTACTAC
>ENST00000390577 havana_ig_gene:known chromosome:GRCh38:12:105895634:105895670:-1 gene:ENSG00000211917 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
GTATTATGATTACGTTTGGGGGAGTTATCGTTATACC
>ENST00000390575 havana_ig_gene:known chromosome:GRCh38:14:105893542:105893561:-1 gene:ENSG00000211915 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
GTGGATACAGCTATGGTTAC
>ENST00000452198 havana_ig_gene:known chromosome:GRCh38:14:105881539:105881556:-1 gene:ENSG00000225825 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
GGGTATAGCAGCGGCTAC
>ENST00000604446 havana_ig_gene:known chromosome:GRCh38:15:21010494:21010516:-1 gene:ENSG00000270824 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
GTGGATATAGTGTCTACGATTAC
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
Última edición por explorer el 2015-03-05 13:02 @585, editado 1 vez en total
Razón: parametro => parámetro; colaboracion => colaboración;
jalapea
Perlero nuevo
Perlero nuevo
 
Mensajes: 7
Registrado: 2015-03-03 22:42 @987
Arriba

Re: Filtrado genes en archivos FASTA

Notapor explorer » 2015-03-05 19:37 @859

¿Tienes hecho algo de código? Es un ejercicio muy sencillo.

Incluso se puede resolver desde la línea de comandos:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. $ perl -ne 'print if /:GRCh38:7:/' archivo.fasta > nuevo_archivo.fasta
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: 14480
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España
Arriba

Re: Filtrado genes en archivos FASTA

Notapor jalapea » 2015-03-05 20:22 @890

explorer, te agradezco mucho.
jalapea
Perlero nuevo
Perlero nuevo
 
Mensajes: 7
Registrado: 2015-03-03 22:42 @987
Arriba

Re: Filtrado genes en archivos FASTA

Notapor jalapea » 2015-03-27 11:15 @510

Un buen día, explorer.

He creado un script para obtener unas secuencias de un archivo fasta.

Con el código anterior que me diste solo me saca el encabezado pero no me saca la secuencia.

Los datos son estos:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
>ENST00000415118 havana_ig_gene:known chromosome:GRCh38:14:22438547:22438554:1 gene:ENSG00000223997 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene
GAAATAGT
>ENST00000448914 havana_ig_gene:known chromosome:GRCh38:12:22449113:22449125:1 gene:ENSG00000228985 gene_biotype:TR_D_gene transcript_biotype:TR_D_gene
ACTGGGGGATACG
>ENST00000431870 havana_ig_gene:known chromosome:GRCh38:7:105894508:105894523:-1 gene:ENSG00000227800 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
TGACTACGGTGACTAC
>ENST00000414852 havana_ig_gene:known chromosome:GRCh38:14:105913222:105913237:-1 gene:ENSG00000233655 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
TGACTACAGTAACTAC
>ENST00000390578 havana_ig_gene:known chromosome:GRCh38:7:105897957:105897987:-1 gene:ENSG00000211918 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
AGGATATTGTAGTGGTGGTAGCTGCTACTCC
>ENST00000390571 havana_ig_gene:known chromosome:GRCh38:14:105886031:105886061:-1 gene:ENSG00000211911 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
GTATTACTATGATAGTAGTGGTTATTACTAC
>ENST00000390577 havana_ig_gene:known chromosome:GRCh38:12:105895634:105895670:-1 gene:ENSG00000211917 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
GTATTATGATTACGTTTGGGGGAGTTATCGTTATACC
>ENST00000390575 havana_ig_gene:known chromosome:GRCh38:14:105893542:105893561:-1 gene:ENSG00000211915 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
GTGGATACAGCTATGGTTAC
>ENST00000452198 havana_ig_gene:known chromosome:GRCh38:14:105881539:105881556:-1 gene:ENSG00000225825 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
GGGTATAGCAGCGGCTAC
>ENST00000604446 havana_ig_gene:known chromosome:GRCh38:15:21010494:21010516:-1 gene:ENSG00000270824 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
GTGGATATAGTGTCTACGATTAC
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4

El resultado que espero en mi archivo de salida es:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
>ENST00000431870 havana_ig_gene:known chromosome:GRCh38:7:105894508:105894523:-1 gene:ENSG00000227800 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
TGACTACGGTGACTAC
>ENST00000390578 havana_ig_gene:known chromosome:GRCh38:7:105897957:105897987:-1 gene:ENSG00000211918 gene_biotype:IG_D_gene transcript_biotype:IG_D_gene
AGGATATTGTAGTGGTGGTAGCTGCTACTCC
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


El script es:

use strict;
use Data::Dumper;

my @list_of_lines;
my @list_out;
my $r=0;
open my $INPUT, "< $ARGV[0]";
open my $OUT, " > $ARGV[0].out";

### Puts every line of the archive in an array called @list_of_lines
while (<($INPUT)>) {
chomp;
push @list_of_lines, $_;
}
CONTINUE:
for (my $i = $r ; $i <= $#list_of_lines; $i++ ) {
### Read until the header is found
if ($list_of_lines[$i] =~ /^>:GRCh38:7:/) {
### $r counts the number of lines before the next header
for ($r = $i; $r <= $#list_of_lines; $r++) {
push @list_out, $list_of_lines[$r];
if ($list_of_lines[$r+1] =~ s/\r\n/n/) {
goto CONTINUE;
}
}
}

}

for (my $d <= $#list_out) {
if ($list_out[$d] =~ /^>/) {
$list_out[$d] =~ s/\s.+//;
}
}
warn Dumper\@list_out;
print $OUT join ("\n",@list_out);

close $INPUT;
close $OUT;
exit;

pero no me funciona. Le agradezco su colaboración.
Última edición por explorer el 2015-03-27 11:52 @536, editado 1 vez en total
Razón: codigo => código; colaboracion => colaboración;
jalapea
Perlero nuevo
Perlero nuevo
 
Mensajes: 7
Registrado: 2015-03-03 22:42 @987
Arriba

Re: Filtrado genes en archivos FASTA

Notapor explorer » 2015-03-27 14:23 @641

Lo tenías casi hecho:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;
  4.  
  5. use Data::Dumper;
  6.  
  7. my @list_of_lines;
  8. my @list_out;
  9. my $r = 0;
  10.  
  11. open my $INPUT ,        "< $ARGV[0]";
  12. open my $OUT,           "> $ARGV[0].out";
  13.  
  14. ### Puts every line of the archive in an array called @list_of_lines
  15. while (<$INPUT>) {                                              # quitar los paréntesis. Con ellos, no estamos leyendo las líneas
  16.                                                                 # asociadas al gestor de archivo $INPUT, sino el propio gestor $INPUT
  17.     chomp;
  18.     push @list_of_lines, $_;
  19. }
  20. warn Dumper \@list_of_lines;
  21.  
  22. CONTINUE:
  23. for ( my $i = $r; $i <= $#list_of_lines; $i++ ) {
  24.     ### Read until the header is found
  25.     if ( $list_of_lines[$i] =~ /:GRCh38:7:/ ) {                 # la expresión regular no coincidía con ninguna línea actual del archivo
  26.         ### $r counts the number of lines before the next header
  27.         for ( $r = $i; $r <= $#list_of_lines; $r++ ) {
  28.             push @list_out, $list_of_lines[$r];
  29.             if ( $list_of_lines[ $r + 1 ] =~ /^>/ ) {           # se ha cambiado una sustitución por la detección de inicio de cabecera
  30.                 goto CONTINUE;
  31.             }
  32.         }
  33.     }
  34.  
  35. }
  36.  
  37. #for my $d (0 .. $#list_out ) {                                 # estas lineas cortan las cabeceras, algo que no queremos
  38. #    if ( $list_out[$d] =~ /^>/ ) {
  39. #        $list_out[$d] =~ s/\s.+//;
  40. #    }
  41. #}
  42.  
  43. warn Dumper \@list_out;
  44.  
  45. print $OUT join( "\n", @list_out), "\n";                        # faltaba un "\n" al final
  46.  
  47. close $INPUT;
  48. close $OUT;
  49. exit;
  50.  
  51. __END__
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4

Ahora bien, como se trata de un caso sencillo en donde tenemos que tomar una decisión por cada línea que leemos, podemos simplificar el procesamiento a... justamente eso: tomar la decisión de pintar o no la línea si hemos leído una cabecera interesante:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use v5.14;
  3. use autodie;
  4.  
  5. ## Leer archivo de entrada
  6. open my $INPUT , '<',  $ARGV[0];
  7. my @lineas = <$INPUT>;
  8. close $INPUT;
  9.  
  10. ## Procesar las líneas
  11. open my $OUTPUT, '>', "$ARGV[0].out";
  12.  
  13. my $lo_hemos_encontrado = 0;
  14.  
  15. for (my $i = 0; $i < @lineas; $i++ ) {          # para todas las líneas
  16.  
  17.     if ($lineas[$i] =~ /^>/) {                  # si hemos encontrado una cabecera
  18.         if ($lineas[$i] =~ /:GRCh38:7:/) {      # y es la que estamos buscando...
  19.             $lo_hemos_encontrado = 1;           # ponemos la bandera de que $lo_hemos_encontrado
  20.         }
  21.         else {                                  # si no, es que se trata de una secuencia que no nos interesa
  22.             $lo_hemos_encontrado = 0;
  23.         }
  24.     }
  25.  
  26.     if ($lo_hemos_encontrado) {                 # mientras $lo_hemos_encontrado
  27.         print $OUTPUT $lineas[$i];              # vamos guardando el resultado
  28.     }
  29. }
  30.  
  31. close $OUTPUT;
  32.  
  33. __END__
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

Pero hay una forma más cómoda de recorrer los valores de un array:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use v5.14;
  3. use autodie;
  4.  
  5. ## Leer archivo de entrada
  6. open my $INPUT , '<',  $ARGV[0];
  7. my @lineas = <$INPUT>;
  8. close $INPUT;
  9.  
  10. ## Procesar las líneas
  11. open my $OUTPUT, '>', "$ARGV[0].out";
  12.  
  13. my $lo_hemos_encontrado = 0;
  14.  
  15. for my $linea (@lineas) {                       # para todas las líneas
  16.  
  17.     if ($linea =~ /^>/) {                       # si la $linea es una cabecera
  18.         if ($linea =~ /:GRCh38:7:/) {           # y es la que estamos buscando...
  19.             $lo_hemos_encontrado = 1;           # ponemos la bandera de que $lo_hemos_encontrado
  20.         }
  21.         else {                                  # si no, es que se trata de una secuencia que no nos interesa
  22.             $lo_hemos_encontrado = 0;
  23.         }
  24.     }
  25.  
  26.     if ($lo_hemos_encontrado) {                 # si $lo_hemos_encontrado
  27.         print $OUTPUT $linea;                   # vamos guardando el resultado
  28.     }
  29. }
  30.  
  31. close $OUTPUT;
  32.  
  33. __END__
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

Incluso podemos ir procesando cada línea, por lo que ni siquiera es necesario almacenarlas en un array. Solo hay que leerlas de una en una:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use v5.14;
  3. use utf8;
  4. use autodie;
  5.  
  6. ## Procesar todas las líneas, una a una
  7. open my $INPUT , '<',  $ARGV[0];
  8. open my $OUTPUT, '>', "$ARGV[0].out";
  9.  
  10. my $lo_hemos_encontrado = 0;
  11.  
  12. while (my $línea = <$INPUT>) {
  13.  
  14.     if ($línea =~ /^>/) {                      # si hemos encontrado una cabecera
  15.         if ($línea =~ /:GRCh38:7:/) {          # y es la que estamos buscando...
  16.             $lo_hemos_encontrado = 1;           # ponemos la bandera de que $lo_hemos_encontrado
  17.         }
  18.         else {                                  # si no, es que se trata de una secuencia que no nos interesa
  19.             $lo_hemos_encontrado = 0;
  20.         }
  21.     }
  22.  
  23.     if ($lo_hemos_encontrado) {                 # mientras $lo_hemos_encontrado
  24.         print $OUTPUT $línea;                  # vamos guardando el resultado
  25.     }
  26. }
  27.  
  28. close $INPUT;
  29. close $OUTPUT;
  30.  
  31. __END__
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

Podemos pensar que un archivo fasta es un conjunto de registros, separados con "\n>". Entonces podemos hacer que Perl lea esos registros de forma completa (no líneas, sino registros).
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use v5.14;
  3. use autodie;
  4.  
  5. ## Procesar todas las líneas, una a una
  6. open my $INPUT , '<',  $ARGV[0];
  7. open my $OUTPUT, '>', "$ARGV[0].out";
  8.  
  9. $/ = "\n>";                             # separador entre registros
  10.  
  11. while (my $registro = <$INPUT>) {
  12.     if ($registro =~ /:GRCh38:7:/) {    # si es el que nos interesa
  13.         $registro =~ s/^>|\n>$//g;      # quitar el marcador '>'
  14.         say $OUTPUT ">$registro";       # resultado
  15.     }
  16. }
  17.  
  18. close $INPUT;
  19. close $OUTPUT;
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

Pero... como Perl es así de potente para el manejo de archivos, también se puede resolver con una línea, desde la consola:
Sintáxis: [ Descargar ] [ Ocultar ]
Using bash Syntax Highlighting
  1. perl -ne 'print and print scalar <> if /:GRCh38:7:/' archivo.fasta > nuevo_archivo.fasta
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4

Funciona así:
  • la opción -n hace que lea el archivo de entrada (archivo.fasta), línea a línea
  • por cada línea, comprobamos si contiene la cadena que nos interesa, con la expresión regular
  • si la condición se cumple, pintamos esa línea (print) y la siguiente (print scalar <>)
  • el operador diamante (<>) significa leer desde la entrada estándar. Con la función scalar(), decimos que solo queremos leer una sola línea. Y print la saca por la salida estándar
  • la salida estándar está redirigida hacia el archivo de salida (nuevo_archivo.fasta)
Claro que esto funciona si la secuencia solo ocupa una línea ;)
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
Arriba

Re: Filtrado genes en archivos FASTA

Notapor rednet » 2016-09-01 09:36 @441

Hola, explorer.

¿Cómo podría hacer lo mismo que hace este script, pero en vez de ser una sola palabra, como en este caso ':GRCh38:7:', pueda leer palabras de búsqueda desde una lista contenida en otro archivo?

Saludos.
rednet
Perlero nuevo
Perlero nuevo
 
Mensajes: 14
Registrado: 2007-11-15 09:15 @427
Arriba

Re: Filtrado genes en archivos FASTA

Notapor explorer » 2016-09-01 12:04 @544

Yo lo que haría sería primero leer todas las palabras de ese archivo.

Luego, en el núcleo principal, donde vamos mirando línea a línea, lo metemos dentro de un bucle for() que vaya recorriendo las palabras, algo así:

for my $palabra (@palabras) {

y luego ya podemos mirar a ver si esa $palabra está en el $registro:

if ($registro =~ /$palabra/) {

y si es así, hacemos la sustitución o lo que sea, dentro del $registro.
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
Arriba

Re: Filtrado genes en archivos FASTA

Notapor rednet » 2016-09-01 15:15 @677

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.     #!/usr/bin/perl
  2.  
  3.     ## Procesar todas las líneas, una a una
  4.     open (F,"<lista_de_secuencias.txt") or die "can't open lista\n";
  5.     my @palabras=<F>;
  6.     close F;
  7.     #print "@palabras||\n";
  8.     open (INPUT , "<$ARGV[0]");
  9.     open (OUTPUT, ">$ARGV[0].out");
  10.      
  11.     $/ = "\n>"; # separador entre registros
  12.     chomp $palabras;
  13.     #print $palabras;
  14.     while(my $registro = <INPUT>) {
  15.         chomp $registro;
  16.         #print "$registro||\n";
  17.         for my $palabra (@palabras) {
  18.             chomp $palabra;
  19.             #print "$registro||$palabra|!!!|\n";
  20.             if ($registro =~ /$palabra/) {    # si es el que nos interesa
  21.                                              # quitar el marcador '>'
  22.                 say OUTPUT ">$registro\n";       # resultado
  23.             }
  24.             last;
  25.         }
  26.     }
  27.     close INPUT;
  28.     close OUTPUT;
  29.     exit;
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

Tengo algunos problemas en script ya que solo toma el primer nombre de la lista, y el resto no. Además no detecta algunas secuencias. Pienso que es por la línea de separador de registro.
rednet
Perlero nuevo
Perlero nuevo
 
Mensajes: 14
Registrado: 2007-11-15 09:15 @427
Arriba
Siguiente
Publicar una respuesta

Volver a Bioinformática

¿Quién está conectado?

Usuarios navegando por este Foro: No hay usuarios registrados visitando el Foro y 1 invitado

Powered by phpBB® Forum Software © phpBB Group
Traducción al Español Mexicano por nextgen en colaboración con phpbb-es.com
phpBB SEO
Advertisements by Advertisement Management
cron