Programa para comparar dos ID y recuperar FASTA
Publicado: 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.
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):
Archivo de salida 3 (comparación de id que no existen):
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.
Using perl Syntax Highlighting
- ########################################################################################################
- # Este programa fue creado para extraer las secuencias FASTA resultado de una comparación de diferentes#
- # identificadores; si el id de la lista existe en el fichero FASTA, imprime el nombre con su secuencia #
- # correspondiente, al mismo tiempo cree un archivo adicional donde imprime el Gi, id y JGI comparados #
- ########################################################################################################
- #!/usr/bin/perl
- use strict;
- use warnings;
- my $file = $ARGV[0] or die "Uso: $0 <lista de ID>\n";
- my $file2 = $ARGV[1] or die "Uso: $1 <fichero FASTA>\n";
- my $DIR = 'effectors_NCBI';
- if ( !-d $DIR ) { # si no existe ese directorio
- mkdir($DIR); # lo creamos
- }
- open( LISTA, "$file" ) || die("No puedo abrir $file\n"); #abre el archivo file1 que es una lista
- my %filas;
- my $n = 0;
- my $numero1;
- while ( my $linea = <LISTA> ) {
- chomp($linea);
- if ( $linea =~ /(\d+)/ ) { # este patrón se puede cambiar para adaptarse a las necesidades
- $numero1 = $1;
- $filas{$numero1} = 0;
- $n++;
- #print ("$numero1\n");
- }
- }
- open( MIA, "<$file2" ) or die("ERROR: No puedo abrir $file: $!\n");
- open( ROB, ">$DIR/Ta_effectors.ncbi.fasta" ) or die("ERROR: No puedo abrir Ta_effectors.ncbi.fasta: $!\n");
- open( RAB, ">$DIR/Ta_ID.ncbi.txt" ) or die("ERROR: No puedo abrir Ta_effectors.ncbi.fasta: $!\n");
- print RAB "\t\tGI\t\t\tID\t\t\tJGI\n";
- my $i = 0;
- my $g = 0;
- my %filas2;
- my $effector;
- my $gi;
- while ( my $linea = <MIA> ) {
- chomp($linea);
- if ( $linea =~ /\>/ ) {
- $linea =~ s/\,//g;
- $linea =~ s/\s/\_/g;
- #print "$linea\n";
- my @a = split( /\_/, $linea );
- $gi = $a[0];
- $effector = $a[4];
- if ( $effector =~ /([\d]+)/ ) {
- $filas2{$effector} = $1;
- $i++;
- #print "$1\n";
- if ( exists $filas{$1} ) { ##aquí pregunto si la comparación existe, también se puede adaptar
- $g++;
- print ROB "$linea\n";
- print RAB "$gi\t$1\tTa\n";
- }
- }
- }
- }
- print RAB "\n\n";
- print "tienes $n jgi en tu archivo\n";
- print "tienes $i gi en tu archivo\n";
- print "encontre $g jgi en tu archivo de gi\n";
- print RAB "tienes $n jgi en tu archivo\n";
- print RAB "tienes $i gi en tu archivo\n";
- print RAB "encontre$g jgi en tu archivo de gi\n";
- close ROB;
- close RAB;
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):
Using text Syntax Highlighting
GI ID_NCBI JGI
358401874 EHK51157.1 186317
358401870 EHK51155.1 232182
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):