• Publicidad

Promotores de Arabidopsis

Perl aplicado a la bioinformática

Promotores de Arabidopsis

Notapor dzavallo » 2013-01-15 09:19 @430

Hola a todos, estoy teniendo un problema con un script para tratar de obtener los promotores de Arabidopsis en formato gff3.

En principio los obtuve simplemente contando -500 nt río arriba del inicio de los genes, pero luego me di cuenta que muchos de ellos se solapan con genes contiguos.

De manera que solo me quiero quedar con los que no solapen con genes.

Para eso escribí un script para que me encuentre coordenadas y si me solapan las exporte a un archivo X y si no a otro Y.

Parto de dos archivos de entrada, el de los genes en formato gff3 y el de los promotores (genes -500 nt) también en gff3.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/local/bin/perl
  2. use strict;
  3. use warnings;
  4.  
  5. my $input1 = "/home/diego-ubuntu/0-archivos_deep_seq/indices/index/TAIR10_solo_prom";
  6. my $input2 = "/home/diego-ubuntu/0-archivos_deep_seq/indices/index/TAIR10_solo_gene";
  7. my $output1 = "/home/diego-ubuntu/0-archivos_deep_seq/indices/index/TAIR10_join_promotores_genes_solapados";
  8. my $output2 = "/home/diego-ubuntu/0-archivos_deep_seq/indices/index/TAIR10_join_promotores_genes_sin_solapar";
  9. open (FILE1, "<$input1.gff3") or die  "could not open $input1 file";
  10. open (OUTPUTFILE1, ">>$output1.gff") or die "could not $output1.txt output file";
  11. open (OUTPUTFILE2, ">>$output2.gff") or die "could not $output2.txt output file";
  12. while (<FILE1>) { # Chr1        TAIR10  gene    3130    3630    .       +       .       ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
  13.         my @split1 = split (/\t/, $_);
  14.         chomp $_;
  15.         my      $chr1 = $split1[0];
  16.         my      $tair1 = $split1[1];
  17.         my      $gene1 = $split1[2];
  18.         my      $coord1 = $split1[3];
  19.         my      $coord2 = $split1[4];
  20.         my      $column5 = $split1[5];
  21.         my      $strand1 = $split1[6];
  22.         my      $column7 = $split1[7];
  23.         my      $id1 = $split1[8];
  24.         open (FILE2, "<$input2.gff") or die "could not open $input2 file";
  25.         while (<FILE2>) { # Chr1        TAIR10  gene    3131    5899    .       +       .       ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
  26.                 my @split2 = split (/\t/, $_);
  27.                 chomp $_;
  28.                 my      $chr2 = $split2[0];
  29.                 my      $tair2 = $split2[1];
  30.                 my      $gene2 = $split2[2];
  31.                 my      $coord3 = $split2[3];
  32.                 my      $coord4 = $split2[4];
  33.                 my      $column55 = $split2[5];
  34.                 my      $strand2 = $split2[6];
  35.                 my      $column77 = $split2[7];
  36.                 my      $id2 = $split2[8];
  37.                 if  ( (($coord1 >= $coord3) and ($coord1 < $coord4)) or (($coord2 >= $coord3) and ($coord2 < $coord4)) )  {
  38.                         if ($chr1 eq $chr2){
  39.                                 print OUTPUTFILE1 "$chr1\t$coord1\t$coord2\t$chr2\t$coord3\t$coord4\t$id1\t$id2";
  40.                                 print "$chr1\t$coord1\t$coord2\t$chr2\t$coord3\t$coord4\t$id1\n";
  41.                         }else {
  42.                                 print OUTPUTFILE2 "$chr1\t$tair1\t\t$gene1\t$coord1\t$coord2\t$column5\t$strand1\t$column7\t$id1";
  43.                                
  44.                         }
  45.                 }
  46.         }
  47. }      
  48.  
  49. exit;
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4


Sé que mis conocimientos no son muy amplios y que debe haber una mejor manera de hacer esto, pero en general uso esta estructura para varios de mis scripts en donde tengo que hacer coincidir coordenadas.

Mis problemas son:
  1. Tuve que poner dos if() porque cuando uso solo uno con todas las condiciones no me reconoce el ($chr1 eq $chr2), no entiendo por qué ya que varias veces usé algo similar y me ha funcionado
  2. Los archivos de salida están bien, pero revisándolos encontré algunos promotores en el archivos "sin_solapar" que también aparecen en el "solapados" y de hecho sí se solapan, con lo cual algo no está funcionando. ¿Es el if() o el else{}? En teoría si un ID aparece en un archivo en el otro no debería aparecer, ¿no es así?
  3. Otro fallo que tiene es que en el archivo "sin_solapar" me aparecen mas de una vez las lineas, es evidente que tengo algo en mi bucle

Si alguien puede decirme en qué está fallando mi script se los agradecería.
dzavallo
Perlero nuevo
Perlero nuevo
 
Mensajes: 12
Registrado: 2013-01-08 11:39 @527

Publicidad

Re: Promotores de Arabidopsis

Notapor explorer » 2013-01-15 10:58 @499

Puedes ahorrar mucho código, en
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.         chomp;
  2.         my($chr1, $tair1, $gene1, $coord1, $coord2, $column5, $strand1, $column7, $id1) = split /\t/;
  3.         # borrar el resto de las líneas
  4.  
  5.  
  6.  
  7.  
  8.  
  9.  
  10.  
  11.  
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

Y en
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.                 chomp;
  2.                 my($chr2, undef, undef, $coord3, $coord4, undef, undef, undef, $id2) = split /\t/;
  3.                 # borrar el resto de las líneas
  4.  
  5.  
  6.  
  7.  
  8.  
  9.  
  10.  
  11.  
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

Como ves hemos quitado las variables $tair2, $gene2, $column55, $strand2 y $column77 porque no se usan luego. Y que el chomp() lo hacemos antes de la partición, para evitar que el último campo contenga el carácter de nueva línea.

Y en
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.                         if ($chr1 eq $chr2) {
  2.                                 my $msg = join("\n", $chr1, $coord1, $coord2, $chr2, $coord3, $coord4, $id1);
  3.                                 print OUTPUTFILE1 "$msg\t$id2\n";
  4.                                 print             "$msg\n";
  5.                         }
  6.                         else {
  7.                                 print OUTPUTFILE2 join("\t", $chr1, $tair1, '', $gene1, $coord1, $coord2, $column5, $strand1, $column7, $id1), "\n";
  8.                                
  9.                         }
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


También, es muy poco eficiente el tener que leer el segundo archivo, continuamente, por cada registro del primer archivo. Lo inteligente es leerlo solo una vez, pasando su información a un array, o mejor aún, a un hash, cuyas claves serán aquél campo que es único por cada registro del segundo archivo (¿cuál es? ¿$chr2?)

Lo del if(), viendo el código, me parece lógico: haces dos operaciones distintas según sean $chr1 igual o distinto a $chr2.

En el tema de la aparición de promotores en el archivo, cuando no deberían, todo depende de la condición puesta en el if(), y recordar que estás revisando todos los registros del segundo archivo con todos los registros del primero, por lo que si una misma condición se repite en todas esas condiciones, se repetirán los print() que sacan la información hacia afuera.

Para resolver esto, a mi me haría falta saber qué campos son únicos en cada archivo (si los hay), o qué hacer en caso de que se detecte una repetición de la condición (¿filtrar o no filtrar?).

Me falta información más detallada de este problema.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Promotores de Arabidopsis

Notapor dzavallo » 2013-01-15 11:29 @520

explorer, muchas gracias otra vez por tu respuesta.
Voy a empezar a implementar el split() como escribiste para ahorrar código.

¿A qué te refieres con campos únicos?

Los archivos tienen exactamente la misma estructura, con los mismos campos.

La primera condición if ($chr1 eq $chr2) es porque necesito que solape las coordenadas correspondientes al mismo cromosoma.

Los archivos tienen en el primer campo todos los cromosomas de Arabidopsis (Chr1, Chr2, Chr3, Chr4, etc). Para cada $chr igual en ambos archivos, tengo que hacer la segunda parte: (($coord1 > $coord3) and ($coord1 < $coord4)) or (($coord2 > $coord4) and ($coord2 > $coord3)), en donde las $coord1 y $coord2 corresponden al promotor y las $coord3 y $coord4 al gen.

Lo que intenté hacer con estas condiciones es que, si se solapan, las guarde en un archivo, y si no cumple esta condición, que guarde las coordenadas del promotor (primer archivo) en otro archivo distinto (OUTPUTFILE2).

El array lo hago para leer las líneas del segundo archivo una sola vez, guardándolas allí.

Gracias por los aportes.
dzavallo
Perlero nuevo
Perlero nuevo
 
Mensajes: 12
Registrado: 2013-01-08 11:39 @527

Re: Promotores de Arabidopsis

Notapor explorer » 2013-01-15 12:53 @578

dzavallo escribiste:¿A qué te refieres con campos únicos?
A que no sabemos si estos archivos contienen alguno de los campos, en las distintas columnas, que podamos considerar únicos (no se repiten en el resto de las líneas.

dzavallo escribiste:La primera condición if ($chr1 eq $chr2) es porque necesito que solape las coordenadas correspondientes al mismo cromosoma.
Bien, esto ya es una pista. La pregunta anterior entonces se puede replantear de la siguiente manera: ¿puede aparecer, en el archivo primero o en el segundo, más de una línea, con el mismo cromosoma? O... ¿cada cromosoma solo aparece una vez en todo el archivo primero y una sola vez en todo el archivo segundo?

dzavallo escribiste:Los archivos tienen en el primer campo todos los cromosomas de Arabidopsis (Chr1, Chr2, Chr3, Chr4, etc). Para cada $chr igual en ambos archivos, tengo que hacer la segunda parte: (($coord1 > $coord3) and ($coord1 < $coord4)) or (($coord2 > $coord4) and ($coord2 > $coord3)), en donde las $coord1 y $coord2 corresponden al promotor y las $coord3 y $coord4 al gen.
Entonces, si se repiten las condiciones de la condición if(), obtendrás varios resultados, desde luego.

dzavallo escribiste:El array lo hago para leer las líneas del segundo archivo una sola vez, guardándolas allí.
Me temo que no... estás leyendo el segundo archivo cada vez por cada línea del primero...

¿Puedes publicar dos trozos pequeños de ejemplo, de cada archivo de entrada, con los correspondientes archivos de salida, en que aparezcan todas las combinaciones posibles de las condiciones que quieres procesar? Viéndolo, será mucho más fácil buscar una solución.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Promotores de Arabidopsis

Notapor dzavallo » 2013-01-15 13:20 @597

Bien, esto ya es una pista. La pregunta anterior entonces se puede replantear de la siguiente manera: ¿puede aparecer, en el archivo primero o en el segundo, más de una línea, con el mismo cromosoma? O... ¿cada cromosoma solo aparece una vez en todo el archivo primero y una sola vez en todo el archivo segundo?


Claro, aparecen muchas veces, miles.

Aquí pego 10 líneas del archivo "solo_prom" correspondiente a los promotores:

Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 8738 9238 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 13715 14215 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 22645 23145 . + . ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
Chr1 TAIR10 gene 27999 28499 . + . ID=AT1G01046;Note=miRNA;Name=AT1G01046
Chr1 TAIR10 gene 33154 33654 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
Chr1 TAIR10 gene 37872 38372 . - . ID=AT1G01060;Note=protein_coding_gene;Name=AT1G01060
Chr1 TAIR10 gene 40945 41445 . - . ID=AT1G01070;Note=protein_coding_gene;Name=AT1G01070
Chr1 TAIR10 gene 44176 44676 . + . ID=AT1G01073;Note=protein_coding_gene;Name=AT1G01073
Chr1 TAIR10 gene 47020 47520 . - . ID=AT1G01080;Note=protein_coding_gene;Name=AT1G01080

Aquí pego 10 líneas del archivo " solo_gene" correspondiente a los genes:

Chr1 TAIR10 gene 3631 5899 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 5928 8737 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 11649 13714 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 28500 28706 . + . ID=AT1G01046;Note=miRNA;Name=AT1G01046
Chr1 TAIR10 gene 31170 33153 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
Chr1 TAIR10 gene 33379 37871 . - . ID=AT1G01060;Note=protein_coding_gene;Name=AT1G01060
Chr1 TAIR10 gene 38752 40944 . - . ID=AT1G01070;Note=protein_coding_gene;Name=AT1G01070
Chr1 TAIR10 gene 44677 44787 . + . ID=AT1G01073;Note=protein_coding_gene;Name=AT1G01073
Chr1 TAIR10 gene 45296 47019 . - . ID=AT1G01080;Note=protein_coding_gene;Name=AT1G01080
Chr1 TAIR10 gene 47485 49286 . - . ID=AT1G01090;Note=protein_coding_gene;Name=AT1G01090

y los archivos de salida son "genes_solapados":

Chr1 33154 33654 Chr1 33379 37871 ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
ID=AT1G01060;Note=protein_coding_gene;Name=AT1G01060
Chr1 47020 47520 Chr1 47485 49286 ID=AT1G01080;Note=protein_coding_gene;Name=AT1G01080
ID=AT1G01090;Note=protein_coding_gene;Name=AT1G01090
Chr1 63812 64312 Chr1 64166 67625 ID=AT1G01130;Note=protein_coding_gene;Name=AT1G01130
ID=AT1G01140;Note=protein_coding_gene;Name=AT1G01140
Chr1 72139 72639 Chr1 72339 74096 ID=AT1G01150;Note=protein_coding_gene;Name=AT1G01150
ID=AT1G01160;Note=protein_coding_gene;Name=AT1G01160
Chr1 71838 72338 Chr1 70115 72138 ID=AT1G01160;Note=protein_coding_gene;Name=AT1G01160
ID=AT1G01150;Note=protein_coding_gene;Name=AT1G01150

y "genes_sin_solapar":

Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 3130 3630 . + . ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
Chr1 TAIR10 gene 8738 9238 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 8738 9238 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 8738 9238 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 8738 9238 . - . ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
Chr1 TAIR10 gene 13715 14215 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 13715 14215 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 13715 14215 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 13715 14215 . - . ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
Chr1 TAIR10 gene 22645 23145 . + . ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
Chr1 TAIR10 gene 22645 23145 . + . ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
Chr1 TAIR10 gene 22645 23145 . + . ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
Chr1 TAIR10 gene 22645 23145 . + . ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
Chr1 TAIR10 gene 27999 28499 . + . ID=AT1G01046;Note=miRNA;Name=AT1G01046
Chr1 TAIR10 gene 27999 28499 . + . ID=AT1G01046;Note=miRNA;Name=AT1G01046
Chr1 TAIR10 gene 27999 28499 . + . ID=AT1G01046;Note=miRNA;Name=AT1G01046
Chr1 TAIR10 gene 33154 33654 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
Chr1 TAIR10 gene 33154 33654 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
Chr1 TAIR10 gene 33154 33654 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050
Chr1 TAIR10 gene 33154 33654 . - . ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050

Remarqué en color uno de los ID como ejemplo. Si observas las coordenadas del promotor AT1G01050 del archivo "solo_prom", se solapa con el gen AT1G01060 del archivo "solo_gene", porque lo que según las condiciones que quiero ajustar, debería archivarse solo en "genes_solapados", como efectivamente aparece, e imprimir los campos que quiero.

Pero no debería aparecer en "genes_sin_solapar" por la condición else{}. Sin embargo también me lo imprime allí (última línea).

Además, como se puede observar, me repite las líneas varias veces, no sé por qué.

Subo un gráfico en donde creo se ve más claramente las combinaciones que quiero hacer.

promotor.jpg

Aquí se ven las dos posibles maneras en que se solapen el promotor y un gen (si está en la hebra positiva o en la negativa). Es por ello que necesito todas esas condiciones juntas y que son excluyentes entre sí. No me debería aparecer en los dos archivos.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. if ( ($chr1 eq $chr2) and (($coord1 >= $coord3) and ($coord1 <= $coord4)) or (($coord2 >= $coord3) and ($coord2 <= $coord4))    ) {
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


Ahora que tal vez me explique mejor, y con las correcciones que me pasaste, vuelvo a mostrar mi script, salvo que en esta oportunidad, el archivo "sin_solapar" repite la primera línea indefinidamente.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/local/bin/perl
  2. use strict;
  3. use warnings;
  4.  
  5. my $input1 = "/home/diego-ubuntu/0-archivos_deep_seq/indices/index/TAIR10_solo_prom";
  6. my $input2 = "/home/diego-ubuntu/0-archivos_deep_seq/indices/index/TAIR10_solo_gene";
  7. my $output1 = "/home/diego-ubuntu/0-archivos_deep_seq/indices/index/TAIR10_join_promotores_genes_solapados1";
  8. my $output2 = "/home/diego-ubuntu/0-archivos_deep_seq/indices/index/TAIR10_join_promotores_genes_sin_solapar1";
  9. open (FILE1, "<$input1.gff3") or die  "could not open $input1 file";
  10. open (OUTPUTFILE1, ">>$output1.gff") or die "could not $output1.txt output file";
  11. open (OUTPUTFILE2, ">>$output2.gff") or die "could not $output2.txt output file";
  12. while (<FILE1>) { # Chr1        TAIR10  gene    3130    3630    .       +       .       ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
  13.         chomp;
  14.    my($chr1, $tair1, $gene1, $coord1, $coord2, $column5, $strand1, $column7, $id1) = split /\t/;       
  15.         open (FILE2, "<$input2.gff") or die "could not open $input2 file";
  16.         while (<FILE2>) { # Chr1        TAIR10  gene    3131    5899    .       +       .       ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
  17.                 chomp;
  18.         my($chr2, undef, undef, $coord3, $coord4, undef, $strand2, undef, $id2) = split /\t/;          
  19.                 if ( ($chr1 eq $chr2) and (($coord1 >= $coord3) and ($coord1 <= $coord4)) or (($coord2 >= $coord4) and ($coord2 <= $coord4))    ) {
  20.                         my $msg = join("\t", $chr1, $coord1, $coord2, $chr2, $coord3, $coord4, $id1);
  21.         print OUTPUTFILE1 "$msg\n";
  22.          print "$msg\n";
  23.                 }       else {         
  24.                                 print OUTPUTFILE2 join("\t", $chr1, $tair1, $gene1, $coord1, $coord2, $column5, $strand1, $column7, $id1), "\n";
  25.            }   
  26.         }      
  27. }      
  28.  
  29. exit;
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4


¿Qué estoy haciendo mal con mis condiciones?

Espero que haya sido más claro que mi anterior explicación, y desde ya agradezco mucho la ayuda.
dzavallo
Perlero nuevo
Perlero nuevo
 
Mensajes: 12
Registrado: 2013-01-08 11:39 @527

Re: Promotores de Arabidopsis

Notapor explorer » 2013-01-15 18:49 @826

Las condiciones están bien, salvo que has puesto dos veces $coord4 en la segunda parte. Te falta un $coord3.

Pero el principal problema no son las condiciones (la lógica de la condición es correcta). El problema está en que recorres todas las líneas de genes por todas las líneas de los promotores, imprimiendo el resultado por cada línea del archivo de genes, cuando en realidad debería imprimirse al final del bucle de los genes, sabiendo entonces si el promotor está solapado o no con alguno de los genes.

Creo que he encontrado una posible solución. Dados los dos primeros archivos de ejemplo, con este programa:
Sintáxis: (solapamiento.pl) [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;
  4. use autodie;
  5.  
  6. # Formato de los archivos
  7. #
  8. # Chr1        TAIR10  gene    3131    5899    .       +       .       ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
  9.  
  10. ## Leemos el archivo de genes
  11. my @genes;
  12. open my $GENES, '<', 'genes.gff3';
  13. while (<$GENES>) {
  14.     chomp;
  15.  
  16.     # solo nos quedaremos con los datos que realmente necesitaremos
  17.     my($chr, undef, undef, $inicio, $final, undef, undef, undef, $id) = split /\t/;
  18.  
  19.     push @genes, [ $chr, $inicio, $final, $id ];
  20. }
  21. close $GENES;
  22.  
  23. ## Recorremos los promotores, buscando solapamientos
  24. open my $SOLAPA, '>',    'solapados.gff3';
  25. open my $NOSOLA, '>', 'no_solapados.gff3';
  26.  
  27. open my $PROMO,  '<',   'promotores.gff3';
  28.  
  29. while (my $promotor = <$PROMO>) {
  30.     chomp $promotor;
  31.     my($promo_chr, undef, undef, $promo_inicio, $promo_final, undef, undef, undef, $promo_id) = split /\t/, $promotor;
  32.  
  33.     my @solapamiento;                                   # guarda cualquier solapamiento que encontremos
  34.  
  35.     for my $gen (@genes) {                              # para todos los genes
  36.         my($gen_chr, $gen_inicio, $gen_final, $gen_id) = @{$gen};
  37.  
  38.         next if $promo_chr ne $gen_chr;                 # no son cromosomas iguales, saltamos a mirar el siguiente gen
  39.  
  40.         next if $promo_inicio > $gen_final;             # no solapado (está por detrás), saltamos
  41.         last if $promo_final  < $gen_inicio;            # no solapado (está por delante), no seguimos buscando
  42.  
  43.         @solapamiento = (
  44.             $promo_chr, $promo_inicio, $promo_final,
  45.             $gen_chr,   $gen_inicio,   $gen_final,
  46.             $promo_id,  $gen_id
  47.         );
  48.        
  49.         print $SOLAPA join("\t", @solapamiento), "\n";
  50.     }
  51.  
  52.     unless (@solapamiento) {                            # si no tenemos algún @solapamiento,
  53.         print $NOSOLA "$promotor\n";                    # sacamos el promotor, tal cual
  54.     }
  55. }
  56.  
  57. close $PROMO;
  58.  
  59. close $SOLAPA;
  60. close $NOSOLA;
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4
Primero leemos el archivo de los genes a memoria, pasándolo a un array de arrays (estructura bidimensional).

Luego, hacemos el doble bucle clásico, con las siguientes características:

  • línea 38: si no coinciden los nombres de los cromosomas, saltamos a comparar al siguiente (next) gen
  • líneas 40 y 41. Estas líneas contienen la misma lógica que tu condición de solapamiento, pero de diferente manera. En lugar de preguntar si el inicio y final del promotor está dentro de las coordenadas de inicio y final del gen, pregunto por la situación contraria: ¿están las coordenadas del promotor fuera del rango de las del gen? En la línea 40 pregunto si el inicio del promotor está más allá del final del gen. Si es así, no está solapado, y paso a analizar el siguiente gen. En la línea 41 pregunto si el final del promotor está antes del inicio del gen. Si es así, tampoco está solapado, y termino el bucle de genes, porque, si supongo que los archivos están ordenados según la posición, los siguientes genes también estarán "por detrás" del final del promotor (y no es posible más solapamientos)
  • línea 43: aquí solo llegamos en caso de un solapamiento. Guardamos los datos que nos interesan en la variable (que nos servirá además, de bandera, fuero del bucle) e imprimimos el solapamiento
  • línea 52: una vez terminado de buscar por todos los genes, preguntamos si hemos encontrado un solapamiento o no. Si no hemos encontrado ninguno, imprimiremos el promotor, como gen no solapado.
Solo lo he podido probar con los archivos de ejemplo que nos has pasado, y sale esto:
Sintáxis: (solapados.gff3) [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
  1. Chr1    33154   33654   Chr1    33379   37871   ID=AT1G01050;Note=protein_coding_gene;Name=AT1G01050    ID=AT1G01060;Note=protein_coding_gene;Name=AT1G01060
  2. Chr1    47020   47520   Chr1    47485   49286   ID=AT1G01080;Note=protein_coding_gene;Name=AT1G01080    ID=AT1G01090;Note=protein_coding_gene;Name=AT1G01090
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
Sintáxis: (no_solapados.gff3) [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
  1. Chr1    TAIR10  gene    3130    3630    .       +       .       ID=AT1G01010;Note=protein_coding_gene;Name=AT1G01010
  2. Chr1    TAIR10  gene    8738    9238    .       -       .       ID=AT1G01020;Note=protein_coding_gene;Name=AT1G01020
  3. Chr1    TAIR10  gene    13715   14215   .       -       .       ID=AT1G01030;Note=protein_coding_gene;Name=AT1G01030
  4. Chr1    TAIR10  gene    22645   23145   .       +       .       ID=AT1G01040;Note=protein_coding_gene;Name=AT1G01040
  5. Chr1    TAIR10  gene    27999   28499   .       +       .       ID=AT1G01046;Note=miRNA;Name=AT1G01046
  6. Chr1    TAIR10  gene    37872   38372   .       -       .       ID=AT1G01060;Note=protein_coding_gene;Name=AT1G01060
  7. Chr1    TAIR10  gene    40945   41445   .       -       .       ID=AT1G01070;Note=protein_coding_gene;Name=AT1G01070
  8. Chr1    TAIR10  gene    44176   44676   .       +       .       ID=AT1G01073;Note=protein_coding_gene;Name=AT1G01073
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4

Como hemos dicho antes, estamos haciendo varias suposiciones, como que los archivos están perfectamente ordenados, según por cromosoma y según por la posición de inicio/final de cada gen.

Date cuenta de una cosa: este programa es prácticamente igual al tuyo, salvo por un detalle básico: solo imprimimos el gen promotor, si no está solapado, fuera del bucle de comparación con el resto de los genes, y en caso de que no haya habido ningún solapamiento con ellos.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Promotores de Arabidopsis

Notapor dzavallo » 2013-01-16 08:23 @391

Una vez más, ¡muchas gracias! Voy a probarlo ya que seguramente corra mucho más rápido que el mio.

Finalmente pude resolverlo agregando una variable (solapado) y la inicio en falso cada vez que lee el solo_prom, y cuando encuentra un solapado imprime en el archivo de salida y le pone verdadero. Al final le puse que si la variable (solapado) es igual a falso, lo imprima en no solapados. Pero me tardó mucho la corrida.


Gracias por la ayuda.
dzavallo
Perlero nuevo
Perlero nuevo
 
Mensajes: 12
Registrado: 2013-01-08 11:39 @527

Re: Promotores de Arabidopsis

Notapor explorer » 2013-01-16 11:48 @533

Dinos a qué velocidad va. Y si quieres, se puede acelerar aún más.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
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 5 invitados

cron