Foro - Perl en Español

La comunidad de programadores en Perl de habla hispana.

Obviar


  • Publicidad

Programa para comparar dos ID y recuperar FASTA

Perl aplicado a la bioinformática
Publicar una respuesta

Programa para comparar dos ID y recuperar FASTA

Notapor alemanmd » 2012-01-31 11:27 @518

Hola a todos, tengo un programa que hice para extraer las secuencias FASTA a partir de una lista, lo que necesito es que el programa use los id de la lista 1 y los compare contra una sección del encabezado de mi archivo FASTA. Si el id de la lista 1 existe entonces imprimir el nombre del archivo FASTA y las secuencias. Además en un archivo adicional imprime el id de la lista 1 con el correspondiente GI. Por último, me gustaría que un archivo adicional me imprimiera los id que no fueron encontrados en el fichero FASTA para poder hacer las búsquedas en otro archivo.

El programa hace correctamente las comparaciones e imprime lo que estoy buscando e inclusive me cuenta cuántos id encontró, pero no puedo hacer que recupere las secuencias pues solo me imprime la cabecera.

Espero me puedan dar un consejo de por dónde puedo irme. Enseguida les dejo mi código y un ejemplo de los archivos de entrada así como un ejemplo de cómo debería ser el archivo de salida.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. ########################################################################################################
  2. # Este programa fue creado para extraer las secuencias FASTA resultado de una comparación de diferentes#
  3. # identificadores; si el id de la lista existe en el fichero FASTA, imprime el nombre con su secuencia #
  4. # correspondiente, al mismo tiempo cree un archivo adicional donde imprime el Gi, id y JGI comparados  #
  5. ########################################################################################################
  6.  
  7. #!/usr/bin/perl
  8. use strict;
  9. use warnings;
  10.  
  11. my $file  = $ARGV[0] or die "Uso: $0 <lista de ID>\n";
  12. my $file2 = $ARGV[1] or die "Uso: $1 <fichero FASTA>\n";
  13.  
  14. my $DIR = 'effectors_NCBI';
  15.  
  16. if ( !-d $DIR ) {                      # si no existe ese directorio
  17.     mkdir($DIR);                       # lo creamos
  18. }
  19.  
  20. open( LISTA, "$file" ) || die("No puedo abrir $file\n");    #abre el archivo file1 que es una lista
  21.  
  22. my %filas;
  23. my $n = 0;
  24. my $numero1;
  25.  
  26. while ( my $linea = <LISTA> ) {
  27.     chomp($linea);
  28.     if ( $linea =~ /(\d+)/ ) {         # este patrón se puede cambiar para adaptarse a las necesidades
  29.         $numero1 = $1;
  30.         $filas{$numero1} = 0;
  31.         $n++;
  32.  
  33.         #print ("$numero1\n");
  34.     }
  35. }
  36.  
  37. open( MIA, "<$file2" )                       or die("ERROR: No puedo abrir $file: $!\n");
  38. open( ROB, ">$DIR/Ta_effectors.ncbi.fasta" ) or die("ERROR: No puedo abrir Ta_effectors.ncbi.fasta: $!\n");
  39. open( RAB, ">$DIR/Ta_ID.ncbi.txt" )          or die("ERROR: No puedo abrir Ta_effectors.ncbi.fasta: $!\n");
  40.  
  41. print RAB "\t\tGI\t\t\tID\t\t\tJGI\n";
  42.  
  43. my $i = 0;
  44. my $g = 0;
  45. my %filas2;
  46. my $effector;
  47. my $gi;
  48.  
  49. while ( my $linea = <MIA> ) {
  50.     chomp($linea);
  51.     if ( $linea =~ /\>/ ) {
  52.         $linea =~ s/\,//g;
  53.         $linea =~ s/\s/\_/g;
  54.  
  55.         #print "$linea\n";
  56.         my @a = split( /\_/, $linea );
  57.         $gi       = $a[0];
  58.         $effector = $a[4];
  59.         if ( $effector =~ /([\d]+)/ ) {
  60.             $filas2{$effector} = $1;
  61.             $i++;
  62.  
  63.             #print "$1\n";
  64.             if ( exists $filas{$1} ) { ##aquí pregunto si la comparación existe, también se puede adaptar
  65.                 $g++;
  66.                 print ROB "$linea\n";
  67.                 print RAB "$gi\t$1\tTa\n";
  68.             }
  69.         }
  70.     }
  71. }
  72.  
  73. print RAB "\n\n";
  74. print "tienes $n jgi en tu archivo\n";
  75. print "tienes $i gi en tu archivo\n";
  76. print "encontre $g jgi en tu archivo de gi\n";
  77. print RAB "tienes $n jgi en tu archivo\n";
  78. print RAB "tienes $i gi en tu archivo\n";
  79. print RAB "encontre$g jgi en tu archivo de gi\n";
  80.  
  81. close ROB;
  82. close RAB;
  83.  
Coloreado en 0.004 segundos, usando GeSHi 1.0.8.4



Lista de ID a buscar
186317
232175
232182
62551
92071
211744


Archivo FASTA a buscar:

>gi|358401874|gb|EHK51157.1| hypothetical protein TRIATDRAFT_186317 [Trichoderma atroviride IMI 206040]
MSLTLVLFLIGILGFVFNRKNIILMLISIEIMLLSITFLILVSSINLDDIIGQTYAIYIIVVAGAESAIG
LAILVAFYRLRGSIAIEYK

>gi|358401873|gb|EHK51156.1| hypothetical protein TRIATDRAFT_232175 [Trichoderma atroviride IMI 206040]
MSSVTLLFIIVSIIALLFLALNFILAPHNPYQEKYSIFECGFHSFLGQNRSQFGVKFFIFALVYLLLDLE
ILLIYPYGMSIYENGLYGLIIMLIFTFIITAGFVFELGKSALKIDSRQSYTYFYKSQKFINTFIENK

>gi|358401870|gb|EHK51155.1| hypothetical protein TRIATDRAFT_232182, partial [Trichoderma atroviride IMI 206040]
MRLLEFSDTKFSFTKDLQDKNIPQYAILSHTWGLDTEEVTYKDLIDGTGMNKAGFKKLQFCGEQAMQDGL
QYFWIDTCCIDKSNSTELNEAITSMFRWYQNATRCYVYLSDVSFPTFDSLQQFNPEVDTIFRASRWFTRG
WTLQELIAPFSVEFFTKEGKLIGNKKSLEQQIHEVTKVAIQALRGESLSEFDIEERFNWADGRQTSREED
LAYCLFGIFDVSIAALYGEGKDKAFRRLRKDI

>gi|358401869|gb|EHK51154.1| hypothetical protein TRIATDRAFT_94382 [Trichoderma atroviride IMI 206040]
MCTMRTAQAITRMVDTGISIIDSAEQIRNSEDSLSQYLDNLLKEVAQEREALSKLGNKLDSDMKKRLNPL
IEKLKQLCDKLLVSPGLKMQRQGKFKIAESAIKDKLDEEDKKSEYTVMRKEFIAFKQELMASKLDNILEI
VANKTKRTYQKLQSLDQLEANYNLIKWNDFTQRSPDEFDELSQRVATLIARLNLGFDFKGQKFDNTEQAA
YGTFDWMVRFDSSIACSTRKLEEKEEEVYKRHNEENLERRAEATHQFRSFLKDDRRVYMVLGKPGSGKST
LMKSLVESPQVKYELESWALEQKKRLIKAHFFFSVTFGSGGLQTEEAMCRDILIQGRADSKRQDQHRL

>gi|358401866|gb|EHK51153.1| hypothetical protein TRIATDRAFT_303121 [Trichoderma atroviride IMI 206040]
MDDNSTLDYVYKRRKIEDTSAQLGADEGTYKKRLRLNDDNDCEHRHKRQKPESPTACITSHASSSSSSSS
LQQLEGENVAEDGTVYILGVGSVRQTALSYHLNLGLRSQ


Archivo de salida 1 (secuencias FASTA identificadas):

>gi|358401874|gb|EHK51157.1| hypothetical protein TRIATDRAFT_186317 [Trichoderma atroviride IMI 206040]
MSLTLVLFLIGILGFVFNRKNIILMLISIEIMLLSITFLILVSSINLDDIIGQTYAIYIIVVAGAESAIG
LAILVAFYRLRGSIAIEYK

>gi|358401873|gb|EHK51156.1| hypothetical protein TRIATDRAFT_232175 [Trichoderma atroviride IMI 206040]
MSSVTLLFIIVSIIALLFLALNFILAPHNPYQEKYSIFECGFHSFLGQNRSQFGVKFFIFALVYLLLDLE
ILLIYPYGMSIYENGLYGLIIMLIFTFIITAGFVFELGKSALKIDSRQSYTYFYKSQKFINTFIENK

>gi|358401870|gb|EHK51155.1| hypothetical protein TRIATDRAFT_232182, partial [Trichoderma atroviride IMI 206040]
MRLLEFSDTKFSFTKDLQDKNIPQYAILSHTWGLDTEEVTYKDLIDGTGMNKAGFKKLQFCGEQAMQDGL
QYFWIDTCCIDKSNSTELNEAITSMFRWYQNATRCYVYLSDVSFPTFDSLQQFNPEVDTIFRASRWFTRG
WTLQELIAPFSVEFFTKEGKLIGNKKSLEQQIHEVTKVAIQALRGESLSEFDIEERFNWADGRQTSREED
LAYCLFGIFDVSIAALYGEGKDKAFRRLRKDI

Archivo de salida 2 (comparación de id que sí existen):
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
    GI             ID_NCBI               JGI
358401874          EHK51157.1           186317
358401870          EHK51155.1           232182
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Archivo de salida 3 (comparación de id que no existen):
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
tienes $numero JGI que no existen en tu fichero FASTA
92071
211744
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
alemanmd
Perlero nuevo
Perlero nuevo
 
Mensajes: 11
Registrado: 2011-11-28 20:16 @886
Arriba

Publicidad

Re: Programa para comparar dos ID y recuperar FASTA

Notapor explorer » 2012-01-31 16:10 @715

De la misma manera que guardas la cabecera en un hash, indexado por el id, puedes guardar la secuencia en otro hash, con el mismo id.
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: Programa para comparar dos ID y recuperar FASTA

Notapor alemanmd » 2012-01-31 17:05 @753

Hola explorer, gracias por el tip. ¿En que sección del código debería poner el hash? Aún no logro imprimir las secuencias con su cabecera correspondiente.

Gracias
Saludos
alemanmd
Perlero nuevo
Perlero nuevo
 
Mensajes: 11
Registrado: 2011-11-28 20:16 @886
Arriba

Re: Programa para comparar dos ID y recuperar FASTA

Notapor explorer » 2012-01-31 17:44 @781

El almacenamiento de las secuencias la puedes hacer como alternativa (else) al if() donde encuentras las líneas de cabecera (/^>/). Y para sacar las cabeceras completas, debes guardarlas, por ejemplo, en el hash.
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: Programa para comparar dos ID y recuperar FASTA

Notapor alemanmd » 2012-01-31 18:33 @814

¡Aaahhh, ya entiendo! Muchas gracias, explorer, ahora mismo trabajo en eso.
Saludos,
Hasta pronto
alemanmd
Perlero nuevo
Perlero nuevo
 
Mensajes: 11
Registrado: 2011-11-28 20:16 @886
Arriba

Re: Programa para comparar dos ID y recuperar FASTA

Notapor explorer » 2012-02-01 21:20 @930

Esta es mi versión:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. # Joaquín Ferrero. 20120201
  3. #
  4. # Extractor de secuencias FASTA entre un fichero con id y un fichero FASTA.
  5. # Las secuencias con id coincidente son escritas en un nuevo archivo, así
  6. # como los id coincidentes, en otro archivo. Las secuencias no coincidentes
  7. # son grabadas en un tercer archivo.
  8. #
  9.  
  10. use Modern::Perl;               # Somos modernos
  11. use utf8::all;                  # Turn on UTF-8. All of it.
  12.  
  13. use autodie;                    # es mejor morir que regresar con deshonor (proverbio Klingon)
  14. use File::Slurp;
  15.  
  16. #### Argumentos #############################################################
  17. @ARGV == 2  or  die "Uso: $0 <archivo con id a buscar> <archivo FASTA>\n";
  18.  
  19. my($nombre_archivo_id, $nombre_archivo_fasta) = @ARGV;
  20.  
  21. -f $nombre_archivo_id     or  die "ERROR: No puedo leer el archivo $nombre_archivo_id: $!\n";
  22. -f $nombre_archivo_fasta  or  die "ERROR: No puedo leer el archivo $nombre_archivo_fasta: $!\n";
  23.  
  24.  
  25. #### Inicialización #########################################################
  26. my $DIR = 'effectors_NCBI';
  27. -d $DIR  or  mkdir($DIR);
  28.  
  29.  
  30. #### Proceso ################################################################
  31. ### Lectura de los id
  32. my @id = read_file($nombre_archivo_id, {chomp => 1});
  33.  
  34.  
  35. ### Apertura de ficheros ####################################################
  36. open my $FASTA  , '<', $nombre_archivo_fasta;
  37. open my $NCBI_EF, '>', "$DIR/Ta_effectors.ncbi.fasta";
  38. open my $NCBI_ID, '>', "$DIR/Ta_ID.ncbi.txt";
  39. open my $NCBI_NO, '>', "$DIR/Ta_NOeffectors.ncbi.fasta";
  40.  
  41. print $NCBI_ID "GI\tID\tJGI\n";
  42.  
  43.  
  44. ### Lectura del fichero FASTA y Escritura del resultado #####################
  45. my($cabecera, $secuencia);
  46. while (my $línea = <$FASTA>) {
  47.  
  48.     if ($línea =~ /^>/  or  eof $FASTA) {
  49.         if ($cabecera   or  eof $FASTA) {
  50.             # >gi|358401874|gb|EHK51157.1| hypothetical protein TRIATDRAFT_186317 [Trichoderma atroviride IMI 206040]
  51.             my(undef, $gi, undef, $id_ncbi, $nombre) = split /\|/, $cabecera;
  52.             my($jgi) = $nombre =~ /_(\d+)/;
  53.  
  54.             if ($jgi ~~ @id) {                          # Si coincide con algún id que buscamos
  55.                 print $NCBI_ID "$gi\t$id_ncbi\t$jgi\n";
  56.                 print $NCBI_EF $cabecera;
  57.                 print $NCBI_EF $secuencia;
  58.                 @id = grep { ! /^$jgi$/ } @id;          # Lo quitamos de la lista
  59.             }
  60.         }
  61.  
  62.         $cabecera  = $línea;
  63.         $secuencia = '';
  64.     }
  65.     else {
  66.         $secuencia .= $línea;
  67.     }
  68. }
  69.  
  70. ### Escritura de las secuencias no coincidentes ############################
  71. print $NCBI_NO "Tienes ", scalar(@id), " JGI que no existen en tu fichero FASTA\n";
  72. for my $id (@id) {
  73.     print $NCBI_NO "$id\n";
  74. }
  75.  
  76.  
  77. ### Cierre de ficheros ######################################################
  78. close $FASTA;
  79. close $NCBI_EF;
  80. close $NCBI_ID;
  81. close $NCBI_NO;
  82.  
  83. __END__
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: Programa para comparar dos ID y recuperar FASTA

Notapor alemanmd » 2012-02-02 16:29 @728

Hola, explorer. Te comento que sí logré recuperar la información que necesitaba con el consejo que me diste. Después vi tu versión del código, y está genial, es más rápida y me gusta más como lo propones tu. Respecto a las librerías que usas (espero que se les llame así), ¿¿¿ son módulos de Bioperl ??? ¿¿¿ Cómo puedo saber qué módulos usar y algún consejo para empezar a utilizar esta herramienta ???

¡¡ Saludos y gracias de nuevo !!
¡¡ Cada vez disfruto más lo que hago gracias a sus consejos y ayuda !!
alemanmd
Perlero nuevo
Perlero nuevo
 
Mensajes: 11
Registrado: 2011-11-28 20:16 @886
Arriba

Re: Programa para comparar dos ID y recuperar FASTA

Notapor explorer » 2012-02-02 16:56 @747

No, no son módulos de BioPerl, son módulos normales y corrientes.

Los módulos los puedes encontrar en CPAN.org. El problema es que hay decenas de miles. Por eso, hay que usar el sistema de búsqueda. Y saber cuáles tienes que usar... pues lo más cómodo y rápido es ver código de otros programadores. Por estos foros hay multitud de códigos de los que se puede aprender mucho. Y si no, siempre puedes ir al foro Módulos y preguntar por alguno que haga algo específico.

Para instalarlos en tu sistema, depende un poco del sistema operativo en que te encuentres, pero hay un procedimiento general para hacerlo.
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: Programa para comparar dos ID y recuperar FASTA

Notapor vliholl » 2016-02-09 18:40 @819

Hola, recupero este tema ya que tengo un problema relacionado con esto.

En mi caso tengo dos archivos: uno de ellos contiene dos columnas; la primera con un código asociado a un grupo de elementos; la segunda son estos elementos separados por '|'.

El segundo archivo es uno tipo fasta donde están las secuencias de cada uno de los elementos mostrados en la columna 2 del archivo 1.

La intención es imprimir el ID de cada grupo, seguido de los elementos de ese grupo y a continuación el código y función de cada uno de ellos, esta última información procedente del archivo 2.

En primer lugar, he separado el archivo 1 en dos, una variable que contiene todas las ID y un hash que contiene los grupos asociados a cada ID.

A continuación intento comparar estos grupos con el archivo2 para imprimir las similitudes junto con la ID y el grupo completo pero no lo consigo.

Soy totalmente novato en esto y no se me ocurre dónde puede estar el error. Os adjunto el código que llevo preparado:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2.  
  3. open( IN, "archivo1.txt" );
  4. while (<IN>) {
  5.     chomp;
  6.     my ( $cluster, $proteinas ) = split( /:/, $_ );
  7.  
  8.     #print "$proteinas\n";
  9.     @listaprot = split( /\|/, $proteinas );
  10.  
  11.     #print "@listaprot\n";
  12.     $prots{$cluster} = $proteinas;
  13. }
  14.  
  15. open( IN, "archivo2.fasta" );
  16. while (<IN>) {
  17.     chomp;
  18.     $linea = <IN>;
  19.     if ( $linea =~ /^>.+/ ) {
  20.         foreach $p ( values %prots ) {
  21.  
  22.             #print "$p-$linea\n";
  23.             if ( $p =~ /$linea/ ) {
  24.                 print "$cluster\t$prots{$cluster}\t$linea\n";
  25.             }
  26.         }
  27.     }
  28. }
  29.  
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


A ver si podéis echarme una mano. Gracias.
vliholl
Perlero nuevo
Perlero nuevo
 
Mensajes: 1
Registrado: 2016-02-09 18:20 @805
Arriba

Re: Programa para comparar dos ID y recuperar FASTA

Notapor explorer » 2016-08-01 04:22 @224

Al código, en principio, no veo problemas, pero sería interesante ver un ejemplo de los archivos que procesa.

¿Puedes publicar un ejemplo de los dos, reducidos a unas pocas líneas?
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
Publicar una respuesta

Volver a Bioinformática

¿Quién está conectado?

Usuarios navegando por este Foro: No hay usuarios registrados visitando el Foro y 3 invitados

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