• Publicidad

Secuencias aleatorias de ADN

Perl aplicado a la bioinformática

Secuencias aleatorias de ADN

Notapor j_f_r_85 » 2010-12-06 06:37 @317

¡¡¡Buen día compañeros!!!

Estoy iniciando en la bioinformática, y pues Perl es el primer peldaño.

Tengo un pequeño problema con un script que comencé por curiosidad. Deseo generar secuencias aleatorias de ADN, con una longitud especificada con el usuario, así como el numero de secuencias. Igualmente el usuario puede definir el % de GC, es decir, qué porcentaje de C y G hay en la cadena.

He creado mi script, pero no me da la salida que deseo; yo deseo algo como:

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
CATCTAAAAG
TGCGTAGATC
GTGTACCAAC
AGGAGGCTTG
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Pero me da algo así:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
CATCTAAAAG
CATCTAAAAGTGCGTAGATC
CATCTAAAAGTGCGTAGATCGTGTACCAAC
CATCTAAAAGTGCGTAGATCGTGTACCAACAGGAGGCTTG
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Sí, me crea las secuencias aleatorias, con el porcentaje de GC deseado, pero no en el patrón que deseo. ¿Me pueden ayudar a encontrar qué hice incorrectamente? ¡Gracias de antemano!

El script:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. #random_nt.pl
  3. #Script generates "n" random DNA sequences given a %GC and length of sequence
  4.  
  5. use strict; use warnings;
  6.  
  7. die "usage: random_nt.pl <gc %> <length> <# sequences>" unless @ARGV == 3;
  8.  
  9. my ($gc, $length, $nseq) = @ARGV;
  10.  
  11. my $seq = '';
  12.  
  13. foreach (my $i = 0; $i < $nseq;   $i++){
  14.         for (my $j = 0; $j < $length; $j++){
  15.                 my $newnt = rndmbp ($gc);
  16.                 $seq .= $newnt;
  17.         }
  18.         print "$seq\n";
  19. }
  20.  
  21. sub rndmbp {
  22. my ($gc)      = @_;
  23. my $GC        = ($gc/2);
  24. my $AT        = (0.5-$GC);
  25. my ($base)    = '';
  26. my $rndm      = rand;
  27. if ($rndm     > ((2*$AT)+$GC))
  28.         {$base   .= "G";}
  29. elsif ($rndm  > (2*$AT))
  30.         {$base   .= "C";}
  31. elsif ($rndm  > $AT)
  32.         {$base   .= "A";}
  33. elsif ($rndm  > 0)
  34.         {$base   .= "T";}
  35. return $base;
  36. }
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4
j_f_r_85
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2010-12-06 06:24 @309

Publicidad

Re: Secuencias aleatorias de ADN

Notapor wanako » 2010-12-06 11:31 @522

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. ## my $seq = ''; ahorramos 8 bytes o_O
  2. for ( 1 .. $nseq ) {
  3.     my $seq;
  4.     for ( 1 .. $length ) {
  5.         my $newnt = rndmbp( $gc );
  6.         $seq .= $newnt;
  7.     }
  8.     print "$seq\n";
  9. }
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4
wanako
Perlero nuevo
Perlero nuevo
 
Mensajes: 27
Registrado: 2010-09-23 11:27 @519

Re: Secuencias aleatorias de ADN

Notapor j_f_r_85 » 2010-12-06 12:14 @551

¡Muchas gracias, wanako! Sabía que la respuesta seguro sería algo simple, pero la verdad no la encontraba, ¡y además ganamos 8 bytes! En serio, ¡muchas gracias!

No sé si se pueda optimizar el script, pero lo importante es que ahora funciona como quiero.

Ahora intentaré hacer algo con esa salida con un nuevo script. Gracias.
j_f_r_85
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2010-12-06 06:24 @309

Re: Secuencias aleatorias de ADN

Notapor explorer » 2010-12-06 13:02 @585

Perdona, j_f_r_85, pero creo que el programa no funciona muy bien.

Si ejecuto ./code_23733_2.pl .50 10 10, es decir, diez secuencias de diez bases, con un 50% de aparición de G y C, los resultados que me salen son:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
CACAATGCTA 4/10
ATGTCTGATT 3/10
CCTACATCGC 6/10
GTTTGTCGGG 6/10
ATGAATTGGA 3/10
AGAAGACTTA 3/10
TCTAGCGACG 6/10
GAACTAAAAG 3/10
CAGGAAACTT 4/10
CCTTTCGCCT 6/10
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
que, como ves, ninguno llega al 50% de aparición.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Secuencias aleatorias de ADN

Notapor explorer » 2010-12-06 15:53 @703

Esta es una solución que he preparado, que sí genera presencias de Cy G acorde a lo que nos pide el usuario.
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. #
  3. # Genera secuencias de ADN de una determinada longitud,
  4. # en la que se indica el porcentaje de aparición de C y G
  5. #
  6. # JF 20101206.2153
  7.  
  8. use common::sense;              # siempre hay que tener sentido común
  9.  
  10. ## Argumentos al programa
  11. @ARGV == 3  or  die "Uso: $0 GC% <longitud de las secuencias> <número de secuencias>\n";
  12.  
  13. my ($GC_porcentaje, $longitud, $número_secuencias) = @ARGV;
  14.  
  15. $longitud          > 0  or  die "ERROR: la longitud debe ser mayor que 0\n";
  16. $número_secuencias > 0  or  die "ERROR: el número de secuencias debe ser mayor que 0\n";
  17. $GC_porcentaje     > 0  or  die "ERROR: el GC% debe estar entre 0 y 1\n";
  18. $GC_porcentaje     < 1  or  die "ERROR: el GC% debe estar entre 0 y 1\n";
  19.  
  20.  
  21. ## Calculamos el número de caracteres de cada parte
  22. my $GC_longitud_caracteres = int ($longitud * $GC_porcentaje);
  23. my $AT_longitud_caracteres = $longitud - $GC_longitud_caracteres;
  24.  
  25.  
  26. ## Bucle por todas las secuencias
  27. for (1 .. $número_secuencias) {
  28.  
  29.     ## Generamos la secuencia
  30.     my @ATCG
  31.         = ( (map { ('C','G')[rand 2] } 1 .. $GC_longitud_caracteres)
  32.           , (map { ('A','T')[rand 2] } 1 .. $AT_longitud_caracteres)
  33.           )
  34.         ;  
  35.  
  36.     my $ATCG;
  37.     while (@ATCG) {
  38.         $ATCG .= splice @ATCG, rand @ATCG, 1;   # extraemos un elemento al azar
  39.     }
  40.  
  41.     ## comprobación
  42.     my $cuenta_CG = $ATCG =~ tr/CG/CG/;         # contamos las C y G en $ATCG
  43.  
  44.     say "[$ATCG] $cuenta_CG/$longitud";
  45. }
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4
que sale, para, por ejemplo ./code_23733_2.pl .32 100 10:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
[ATTTGTAGCTTTCTAAAACTATATACAAGAGAGGTATTCATTTCCCATTCAAGCATAGGTGTAGCACTGAAACTTAACCGAGTTACCAAATTAATATTTG] 32/100
[CACGTAATCTGTAGTTATAGATCGTTCGTTGAGTTTACAGTGCAAAGGTTTACTAAGAATTCTAACGGGTTTCTCATTTAACTATCATTTTATAAGTTGT] 32/100
[TTAAAGATAAGAAGCTTAAGAATTTAAACTGGAGTCTTCGTAACATGTACGAAAAAACGCCGAGCTTATGTTTTTGTCATTTATGATCAAGGCTATGATT] 32/100
[GAAGTCTTCGGGCATGCCCGTATCTGAACGTGTTATAGAATTTATTAACTATAGACATCTATATTTATAACTAAAACTCGTCTATAATCTAGTCATTCAA] 32/100
[AACGATTTGCTGATAAAAGTGACCGGTTTCACATATTGCATACATTTTATTTTTAGACTATTCCCGTTGAATGACTTCGTGGATTGTTACTTTACTAATA] 32/100
[TAACATAGTAATCAAGGTCATGTTCAGGTAAGAATCGGTTTTCAAATTATTGAGGTTGTTATGGTCCTGTCCCTATACATCTTTTTATGAAATATCAGTT] 32/100
[AAAAGTTGATATTCAAAAAATACGAGCGGGGATACAATTTGTCAAATGAAAATAATCTACCATTGCTTTATCGTTAGACATACTTTAAGCGCTCCAGTTT] 32/100
[AAAACATAACTAAATCAATTCCATCTGTAGTAGGTATGTATTCACAGACACGACACAGCGATTTGAACGTCTATTTGTTATAAATCATTAAGGGTCAAAT] 32/100
[CAAGTAATCATGTACCTATAAGAAATAAAATTATTGACTTGTATTCAGTTTTACTATTCCCTAGCAGTCAAAAAATGGTGCGGCAAGTAATAACGTGAGA] 32/100
[TTTTTGCAAGTTATTAAGAATTAAGTAAGTGGTATAAAAGTGAACACGTTTTTATAACTTTAGAGTGAGCCTCGACTGCGTGATTAAATAGAAGGGAACA] 32/100
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4

La generación de la secuencia se hace con un array, lo cual no es lo recomendable si queremos generar secuencias de varios cientos de millones de bases. Pero bueno, para cosas pequeñas, puede valer :)

De todas maneras, creo que esto lo hace de forma directa el módulo Bio::Matrix::PSM::SiteMatrix.

Y me sonaba que había otro módulo de matemáticas que podía generar secuencias según la frecuencia que le indiquemos. Lo he intentando encontrar, pero nada...
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Secuencias aleatorias de ADN

Notapor j_f_r_85 » 2010-12-06 16:01 @709

Mi estimado explorer:

Inmediatamente se nota tu sapiencia al respecto del tema tu script luce muy profesional, y como tu dices ¡con sentido común!

¡Muchas gracias!

Tengo que seguir practicando.
j_f_r_85
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2010-12-06 06:24 @309

Re: Secuencias aleatorias de ADN

Notapor explorer » 2010-12-06 16:28 @728

common::sense es un módulo muy curioso.

Ya sé que los puristas del lenguaje me dirán que no hay que usarlo, y menos en unos foros, y menos en el foro básico, pero el caso es que tiene justo lo que quiero para todos mis programas: utf8, las nuevas características de Perl v5.10 y siguientes, strict, warnings, etc, etc.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
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 0 invitados

cron