• Publicidad

Identificar agrupación de probables sitios de unión

Perl aplicado a la bioinformática

Identificar agrupación de probables sitios de unión

Notapor K-lixto » 2011-05-16 14:54 @662

Hola tengo un archivo con el siguiente estilo de filas por columnas:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
  1. 2L      bound_moiety=MED(+)     intron  4336355 4336361 3.8     +       .       ID=intron_FBgn0020762:1_FBgn0020762:3;Name=Atet-in;Parent=FBtr0077427;parent_type=mRNA
  2. 2L      bound_moiety=MED(+)     intron  4334989 4334995 3.0     +       .       ID=intron_FBgn0020762:1_FBgn0020762:3;Name=Atet-in;Parent=FBtr0077427;parent_type=mRNA
  3. 2L      bound_moiety=MED(+)     intron  4334462 4334468 2.9     +       .       ID=intron_FBgn0020762:1_FBgn0020762:3;Name=Atet-in;Parent=FBtr0077427;parent_type=mRNA
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Quisiera poder saber si existen grupos de coordenadas a distancia menores de 50 bp teniendo en cuenta las coordenadas de las columnas 4 (inicio) y 5 (fin), y si es positivo poder guardar un archivo que me indique qué motivo (columna 2) está en grupo en la región (columna 3).

Le agradezco desde ya su ayuda.

Saludos
K-lixto
Perlero nuevo
Perlero nuevo
 
Mensajes: 9
Registrado: 2011-05-13 14:46 @657

Publicidad

Re: Identificar agrupación de probables sitios de unión

Notapor explorer » 2011-05-16 15:09 @672

Bienvenido a los foros de Perl en español, K-lixto.

No me queda claro el cálculo que hay que hacer...

¿Puedes poner un ejemplo de cálculo de las coordenadas?

Aquí hay una posible solución:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use Modern::Perl;                       # somos modernos
  3. use autodie;                            # es mejor morir que regresar con deshonor (proverbio Klingon)
  4.  
  5. open my $FH, q[<], 'kk.txt';
  6.  
  7. while (<$FH>) {
  8.     my @campos = split;
  9.    
  10.     if ($campos[3] - $campos[4] > 50) {
  11.         say "$campos[1] $campos[2]";
  12.     }
  13. }
  14.  
  15. close $FH;
  16. __END__
Coloreado en 0.003 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

Re: Identificar agrupación de probables sitios de unión

Notapor K-lixto » 2011-05-19 13:38 @610

¡¡Gracias por la ayuda!! La idea básica es comparar dentro del mismo archivo (que contiene coordenadas de motivos cortos dentro de una secuencia larga cromosómica, ej: 1..1000 000) cuáles de estos motivos, según sus coordenadas, se pueden agrupar como un grupo si entre ellos hay menos de 50 de distancia entre sus coordenadas...
Yo he avanzado con el siguiente código, no soy experto pero me ha resultado:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/env perl
  2.  
  3. use warnings;
  4. use strict;
  5.  
  6. my @arg;
  7. my $n;
  8.  
  9. my $dist;
  10.  
  11. open (OUTPUT, ">>prueba_identificar_cluster_motifs_altenativa.gff" );
  12. print OUTPUT "Chr\tFT\tFeature\tstart\tend\t#Cluster\tCluster\tDistancia\tDescripcion\n";  
  13. open (ARQ, $ARGV[0]) || die "can not open the file"; #HMM_MMZB_all-chromosome-genes_expressed-AS_sorted.gff
  14. while ( my $lastline =<ARQ>){
  15. my ($chr1, $BS1, $FT1, $s1, $e1, $scr1, $a1,  $dir1, $gene1)= split /\t/, $lastline;
  16.  
  17. $n =0; 
  18.  
  19.        open (LIST, $ARGV[1]) || die "can not open the file";#HMM_MMZB_all-chromosome-genes_expressed-AS_sorted.gff
  20.         while ( my $thisline =<LIST>){
  21.              
  22.                chomp $thisline;
  23.                 my ($chr, $BS, $FT, $s,  $e, $scr, $a, $dir, $gene) = split /\t/, $thisline;
  24.                 $dist = $s - $e1;
  25.               unless ( $thisline eq $lastline ){
  26.                      if ( $chr1 eq $chr && $BS1 eq $BS && $dist >0 && $dist <= 50 && $s != $s1 && $e != $e1) {
  27.                      $n++;
  28.                      open (OUTPUT, ">>prueba_identificar_cluster_motifs_altenativa.gff" );      
  29.                      print OUTPUT "$chr1\t$BS1\t$FT\t$s\t$e\t$n\t$s1-$e1\t$dist\t$gene\n";                                                            
  30.                      }else{next;}
  31.               }
  32.        }      
  33.  }
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4
K-lixto
Perlero nuevo
Perlero nuevo
 
Mensajes: 9
Registrado: 2011-05-13 14:46 @657

Re: Identificar agrupación de probables sitios de unión

Notapor explorer » 2011-05-20 17:36 @775

El programa no es óptimo: estás releyendo un fichero por cada línea del primer fichero.

Lo ideal sería mantener en memoria toda la información, leyendo los ficheros solo una vez.
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

Re: Identificar agrupación de probables sitios de unión

Notapor K-lixto » 2011-05-23 14:58 @665

Ciertamente, no es óptimo porque demora muchísimo. ¿Cómo sería la opción de guardar el archivo en memoria y leerlo solo una vez?

Saludos.
K-lixto
Perlero nuevo
Perlero nuevo
 
Mensajes: 9
Registrado: 2011-05-13 14:46 @657

Re: Identificar agrupación de probables sitios de unión

Notapor explorer » 2011-05-23 16:05 @712

Esta es una posible solución, pero como no tengo los ficheros de entrada, no sé si estará bien la salida.
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;
  4.  
  5. ### Cabecera
  6. open my $OUTPUT, '>>', 'prueba_identificar_cluster_motifs_altenativa.gff';
  7. print   $OUTPUT "Chr\tFT\tFeature\tstart\tend\t#Cluster\tCluster\tDistancia\tDescripcion\n";
  8.  
  9. ### Leer el primer archivo
  10. open my $ARQ, '<', $ARGV[0] or die "can not open the file $ARGV[0]";
  11. my @arq = <$ARQ>;
  12. close $ARQ;
  13.  
  14. ### Leer el segundo archivo
  15. open my $LIST, '<', $ARGV[1] or die "can not open the file $ARGV[1]";
  16. my @list = <$LIST>;
  17. close $LIST;
  18.  
  19. ### Comparación
  20. for my $lastline (@arq) {                       # para todas las líneas del primer archivo
  21.  
  22.     my $n = 0;
  23.  
  24.     for my $thisline (@list) {                  # para todas las líneas del segundo archivo
  25.  
  26.         if ( $thisline ne $lastline ) {         # solo miramos si las líneas son distintas
  27.  
  28.             my ($chr1, $BS1, $FT1,  $s1, $e1, undef, undef, undef, $gene1) = split /\t/, $thisline;
  29.             my ($chr2, $BS2, undef, $s2, $e2, undef, undef, undef, undef ) = split /\t/, $lastline;
  30.  
  31.             my $dist  = $s1 - $e2;
  32.  
  33.             if (
  34.                  $chr1  eq  $chr2
  35.             and  $BS1   eq  $BS2
  36.             and  $dist  >   0    
  37.             and  $dist  <=  50  
  38.             and  $s1    !=  $s2  
  39.             and  $e1    !=  $e2  
  40.             ) {
  41.                 $n++;
  42.                 print $OUTPUT join("\t", $chr1, $BS1, $FT1, $s1, $e1, $n, $s2-$e2, $dist, $gene1);
  43.             }
  44.         }
  45.     }
  46. }
  47.  
  48. close $OUTPUT;
  49.  
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4


Notar que $gene1 lleva el carácter de retorno de carro, porque no se lo hemos quitado antes (no hay chomp).
También notar que solo creamos las variables que vamos a usar en la comparación en la impresión del resultado. El resto de campos los llevamos a 'undef', es decir, quedan descartados automáticamente.
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

Re: Identificar agrupación de probables sitios de unión

Notapor K-lixto » 2011-05-24 16:24 @725

Muy agradecido de la mejora del script, pero tengo un tema sin poder resolver: tengo un archivo del cual solo quiero algunas líneas; por ejemplo del siguiente archivo las líneas de la 1 a la 4 es igual a las líneas de la 5 a la 8, respectivamente, aunque la última columna (desde "ID=...") sea diferente. Yo solo quisiera guardar las 4 primeras ya que la repetición de las líneas está dada por los números que vienen a continuación de la columna que dice "intron":

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
  1. 1-2L      bound_moiety=ZEN        intron  8422674 8422679 1       8422665-8422670 4       ID=intron_CG13388:1_CG13388:2;Name=Akap200-in;Parent=FBtr0079664,FBtr0079667;parent_type=mRNA
  2. 2-2L      bound_moiety=ZEN        intron  8422694 8422699 2       8422665-8422670 24      ID=intron_CG13388:1_CG13388:2;Name=Akap200-in;Parent=FBtr0079664,FBtr0079667;parent_type=mRNA
  3. 3-2L      bound_moiety=ZEN        intron  8422703 8422708 3       8422665-8422670 33      ID=intron_CG13388:1_CG13388:2;Name=Akap200-in;Parent=FBtr0079664,FBtr0079667;parent_type=mRNA
  4. 4-2L      bound_moiety=ZEN        intron  8422717 8422722 4       8422665-8422670 47      ID=intron_CG13388:1_CG13388:2;Name=Akap200-in;Parent=FBtr0079664,FBtr0079667;parent_type=mRNA
  5. 5-2L      bound_moiety=ZEN        intron  8422674 8422679 5       8422665-8422670 4       ID=intron_CG13388:7_CG13388:2;Name=Akap200-in;Parent=FBtr0079665;parent_type=mRNA
  6. 6-2L      bound_moiety=ZEN        intron  8422694 8422699 6       8422665-8422670 24      ID=intron_CG13388:7_CG13388:2;Name=Akap200-in;Parent=FBtr0079665;parent_type=mRNA
  7. 7-2L      bound_moiety=ZEN        intron  8422703 8422708 7       8422665-8422670 33      ID=intron_CG13388:7_CG13388:2;Name=Akap200-in;Parent=FBtr0079665;parent_type=mRNA
  8. 8-2L      bound_moiety=ZEN        intron  8422717 8422722 8       8422665-8422670 47      ID=intron_CG13388:7_CG13388:2;Name=Akap200-in;Parent=FBtr0079665;parent_type=mRNA
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Gracias por la ayuda
K-lixto
Perlero nuevo
Perlero nuevo
 
Mensajes: 9
Registrado: 2011-05-13 14:46 @657

Re: Identificar agrupación de probables sitios de unión

Notapor explorer » 2011-05-24 16:40 @736

Bueno, es muy sencillo, y además está indicado en la documentación de Perl: usando un hash se puede saber si una línea está repetida o no.
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use Modern::Perl;                       # somos modernos
  3. #use utf8;                               # este programa está escrito en utf8
  4. use autodie;                            # es mejor morir que regresar con deshonor (proverbio Klingon)
  5.  
  6. my %campo3_vistos;
  7.  
  8. open my $FH, '<', 'code_26274.txt';
  9.  
  10. while (<$FH>) {
  11.     # imprime la línea si el campo que está en la cuarta columna no lo hemos visto antes
  12.     print if not $campo3_vistos{ (split)[3] }++;
  13. }
  14.  
  15. close   $FH;
  16.  
  17. __END__
  18. 1-2L      bound_moiety=ZEN        intron  8422674 8422679 1       8422665-8422670 4       ID=intron_CG13388:1_CG13388:2;Name=Akap200-in;Parent=FBtr0079664,FBtr0079667;parent_type=mRNA
  19. 2-2L      bound_moiety=ZEN        intron  8422694 8422699 2       8422665-8422670 24      ID=intron_CG13388:1_CG13388:2;Name=Akap200-in;Parent=FBtr0079664,FBtr0079667;parent_type=mRNA
  20. 3-2L      bound_moiety=ZEN        intron  8422703 8422708 3       8422665-8422670 33      ID=intron_CG13388:1_CG13388:2;Name=Akap200-in;Parent=FBtr0079664,FBtr0079667;parent_type=mRNA
  21. 4-2L      bound_moiety=ZEN        intron  8422717 8422722 4       8422665-8422670 47      ID=intron_CG13388:1_CG13388:2;Name=Akap200-in;Parent=FBtr0079664,FBtr0079667;parent_type=mRNA
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


Volver a Bioinformática

¿Quién está conectado?

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

cron