Foro - Perl en Español

La comunidad de programadores en Perl de habla hispana.

Obviar


  • Publicidad

Modelos ocultos de Márkov

Perl aplicado a la bioinformática
Publicar una respuesta

Re: Modelos ocultos de Márkov

Notapor zipetardix » 2012-02-10 19:42 @862

Hola buenas,

Pues yo tengo un problema que en principio tiene que ver con el tema aquí planteado, aunque no es lo que aquí se
ha tratado de resolver. Me explico.

Mi modelo de cadena de Márkov de orden 1, trata de, en una secuencia de ADN, estudiar la frecuencia con la que
cada nucleótido tiende a preceder al siguiente con el objetivo de comprobar diferencia de éstas frecuencias entre
la zona llamada "isla CpG" y la que no.

Pues bien, ya he conseguido escribir todo el script que me da dichas frecuencias. Ahora en el siguiente
paso tengo que implementar un algoritmo de búsqueda en otro script en el que, leyendo las probabilidades
calculadas en el anterior (asignadas tanto para "isla CpG", como para no), calcule la puntuación para islas CpG de
una longitud dada a lo largo de la secuencia. Así pues, el script debe aceptar el input de la
longitud de dicha "isla", que lo escribe el usuario.

Adjunto el script que he hecho para leer las secuencias y que me devuelva las frecuencias:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;
  4.  
  5. my $file1 = $ARGV[0];
  6.  
  7. my %chain1;                     # background seq
  8. my %chain2;                     # CpG island
  9. my %first_nocpg_first;          # 1st nucleotide of the 1st background seq
  10. my %first_nocpg_second;         # 1st nucleotide of the 2nd background seq
  11. my %first_yescpg;               # 1st nucleotide of the CpG island
  12.  
  13. open(my $FILE1, '<', $file1) or die "cannot open $file1: $!\n";
  14.  
  15. while (<$FILE1>) {
  16.     chomp;
  17.    
  18.    
  19.  
  20.     my ($first_col, $second_col, $seq) = split("\t", $_);
  21.     my $cpg_start = $first_col; #cpg_start = primer nt que ES cpg
  22.     my $cpg_end = $second_col-1;        #cpg_end = último nt que ES cpg                    
  23.     my @nt_seq = split(//, $seq);
  24.     my $len = scalar @nt_seq;
  25.                 my @firstnt_yescpg = $nt_seq[$cpg_start];
  26.                 my @firstnt_nocpg_first = $nt_seq[0];
  27.                 my @firstnt_nocpg_second = $nt_seq[$second_col];
  28.                
  29.         my @not_cpg = @nt_seq[0..499, $second_col..$len];
  30.         my @yes_cpg = @nt_seq[$cpg_start..$cpg_end];
  31.  
  32.         my $i=1;
  33.     for my $nt_nocpg (@not_cpg) {
  34.         $chain1{$not_cpg[$i]}{$not_cpg[$i-1]}++;
  35.         $i++;                  
  36.                  }
  37.         my $f=1;          
  38.     for my $nt_yescpg (@yes_cpg) {
  39.         $chain2{$yes_cpg[$f]}{$yes_cpg[$f-1]}++;
  40.         $f++;
  41.                 }
  42.  
  43.         my $j=0;
  44. for my $firstnt_nocpg (@firstnt_nocpg_first) {
  45.         $first_nocpg_first{$firstnt_nocpg_first[$j]} += @firstnt_nocpg_first;
  46.         $j++;
  47.         }
  48.  
  49.         my $k=0;
  50. for my $firstnt_nocpg (@firstnt_nocpg_second) {
  51.         $first_nocpg_second{$firstnt_nocpg_second[$k]} += @firstnt_nocpg_second;
  52.         $k++;
  53.         }
  54.  
  55.         my $l=0;
  56. for my $firstntyescpg (@firstnt_yescpg) {
  57.         $first_yescpg{$firstnt_yescpg[$l]} += @firstnt_yescpg;
  58.         $l++;
  59.         }
  60.  
  61.  
  62. }
  63.  
  64. my $nAA = $chain1{A}{A} / ($chain1{A}{A} + $chain1{C}{A} + $chain1{G}{A} + $chain1{T}{A});
  65. my $nCA = $chain1{C}{A} / ($chain1{A}{A} + $chain1{C}{A} + $chain1{G}{A} + $chain1{T}{A});
  66. my $nGA = $chain1{G}{A} / ($chain1{A}{A} + $chain1{C}{A} + $chain1{G}{A} + $chain1{T}{A});
  67. my $nTA = $chain1{T}{A} / ($chain1{A}{A} + $chain1{C}{A} + $chain1{G}{A} + $chain1{T}{A});
  68.  
  69. my $nAC = $chain1{A}{C} / ($chain1{A}{C} + $chain1{C}{C} + $chain1{G}{C} + $chain1{T}{C});
  70. my $nCC = $chain1{C}{C} / ($chain1{A}{C} + $chain1{C}{C} + $chain1{G}{C} + $chain1{T}{C});
  71. my $nGC = $chain1{G}{C} / ($chain1{A}{C} + $chain1{C}{C} + $chain1{G}{C} + $chain1{T}{C});
  72. my $nTC = $chain1{T}{C} / ($chain1{A}{C} + $chain1{C}{C} + $chain1{G}{C} + $chain1{T}{C});
  73.  
  74. my $nAG = $chain1{A}{G} / ($chain1{A}{G} + $chain1{C}{G} + $chain1{G}{G} + $chain1{T}{G});
  75. my $nCG = $chain1{C}{G} / ($chain1{A}{G} + $chain1{C}{G} + $chain1{G}{G} + $chain1{T}{G});
  76. my $nGG = $chain1{G}{G} / ($chain1{A}{G} + $chain1{C}{G} + $chain1{G}{G} + $chain1{T}{G});
  77. my $nTG = $chain1{G}{G} / ($chain1{A}{G} + $chain1{C}{G} + $chain1{G}{G} + $chain1{T}{G});
  78.  
  79. my $nAT = $chain1{A}{T} / ($chain1{A}{T} + $chain1{C}{T} + $chain1{G}{T} + $chain1{T}{T});
  80. my $nCT = $chain1{C}{T} / ($chain1{A}{T} + $chain1{C}{T} + $chain1{G}{T} + $chain1{T}{T});
  81. my $nGT = $chain1{G}{T} / ($chain1{A}{T} + $chain1{C}{T} + $chain1{G}{T} + $chain1{T}{T});
  82. my $nTT = $chain1{G}{T} / ($chain1{A}{T} + $chain1{C}{T} + $chain1{G}{T} + $chain1{T}{T});
  83.  
  84. my $preA = $first_nocpg_first{A} + $first_nocpg_second{A};
  85. my $preC = $first_nocpg_first{C} + $first_nocpg_second{C};
  86. my $preG = $first_nocpg_first{G} + $first_nocpg_second{G};
  87. my $preT = $first_nocpg_first{T} + $first_nocpg_second{T};
  88.  
  89. my $nA = $preA / ($preA + $preC + $preG + $preT);
  90. my $nC = $preC / ($preA + $preC + $preG + $preT);
  91. my $nG = $preG / ($preA + $preC + $preG + $preT);
  92. my $nT = $preT / ($preA + $preC + $preG + $preT);
  93. print "\nFrequencies outside CpG island:\n";
  94. print "A: $nA", "\n";
  95. print "C: $nC", "\n";
  96. print "G: $nG", "\n";
  97. print "T: $nT", "\n\n";
  98.  
  99. print "A|A: $nAA", "\t", "A|G: $nAG","\n";
  100. print "C|A: $nCA", "\t", "C|G: $nCG","\n";
  101. print "G|A: $nGA", "\t", "G|G: $nGG","\n";
  102. print "T|A: $nTA", "\t", "T|G: $nTG","\n";
  103. print "A|C: $nAC", "\t", "A|T: $nAT","\n";
  104. print "C|C: $nCC", "\t", "C|T: $nCT","\n";
  105. print "G|C: $nGC", "\t", "G|T: $nGT","\n";
  106. print "T|C: $nTC", "\t", "T|T: $nTT","\n\n";
  107.  
  108.  
  109. my $yAA = $chain2{A}{A} / ($chain2{A}{A} + $chain2{C}{A} + $chain2{G}{A} + $chain2{T}{A});
  110. my $yCA = $chain2{C}{A} / ($chain2{A}{A} + $chain2{C}{A} + $chain2{G}{A} + $chain2{T}{A});
  111. my $yGA = $chain2{G}{A} / ($chain2{A}{A} + $chain2{C}{A} + $chain2{G}{A} + $chain2{T}{A});
  112. my $yTA = $chain2{T}{A} / ($chain2{A}{A} + $chain2{C}{A} + $chain2{G}{A} + $chain2{T}{A});
  113.  
  114. my $yAC = $chain2{A}{C} / ($chain2{A}{C} + $chain2{C}{C} + $chain2{G}{C} + $chain2{T}{C});
  115. my $yCC = $chain2{C}{C} / ($chain2{A}{C} + $chain2{C}{C} + $chain2{G}{C} + $chain2{T}{C});
  116. my $yGC = $chain2{G}{C} / ($chain2{A}{C} + $chain2{C}{C} + $chain2{G}{C} + $chain2{T}{C});
  117. my $yTC = $chain2{T}{C} / ($chain2{A}{C} + $chain2{C}{C} + $chain2{G}{C} + $chain2{T}{C});
  118.  
  119. my $yAG = $chain2{A}{G} / ($chain2{A}{G} + $chain2{C}{G} + $chain2{G}{G} + $chain2{T}{G});
  120. my $yCG = $chain2{C}{G} / ($chain2{A}{G} + $chain2{C}{G} + $chain2{G}{G} + $chain2{T}{G});
  121. my $yGG = $chain2{G}{G} / ($chain2{A}{G} + $chain2{C}{G} + $chain2{G}{G} + $chain2{T}{G});
  122. my $yTG = $chain2{G}{G} / ($chain2{A}{G} + $chain2{C}{G} + $chain2{G}{G} + $chain2{T}{G});
  123.  
  124. my $yAT = $chain2{A}{T} / ($chain2{A}{T} + $chain2{C}{T} + $chain2{G}{T} + $chain2{T}{T});
  125. my $yCT = $chain2{C}{T} / ($chain2{A}{T} + $chain2{C}{T} + $chain2{G}{T} + $chain2{T}{T});
  126. my $yGT = $chain2{G}{T} / ($chain2{A}{T} + $chain2{C}{T} + $chain2{G}{T} + $chain2{T}{T});
  127. my $yTT = $chain2{G}{T} / ($chain2{A}{T} + $chain2{C}{T} + $chain2{G}{T} + $chain2{T}{T});
  128.  
  129. my $ypreA = $first_yescpg{A};
  130. my $ypreC = $first_yescpg{C};
  131. my $ypreG = $first_yescpg{G};
  132. my $ypreT = $first_yescpg{T};
  133.  
  134. my $yA = $ypreA / ($ypreA + $ypreC + $ypreG + $ypreT);
  135. my $yC = $ypreC / ($ypreA + $ypreC + $ypreG + $ypreT);
  136. my $yG = $ypreG / ($ypreA + $ypreC + $ypreG + $ypreT);
  137. my $yT = $ypreT / ($ypreA + $ypreC + $ypreG + $ypreT);
  138.  
  139. print "\nFrequencies inside CpG island:\n";
  140. print "A: $yA", "\n";
  141. print "C: $yC", "\n";
  142. print "G: $yG", "\n";
  143. print "T: $yT", "\n\n";
  144.  
  145. print "A|A: $yAA", "\t", "A|G: $yAG","\n";
  146. print "C|A: $yCA", "\t", "C|G: $yCG","\n";
  147. print "G|A: $yGA", "\t", "G|G: $yGG","\n";
  148. print "T|A: $yTA", "\t", "T|G: $yTG","\n";
  149. print "A|C: $yAC", "\t", "A|T: $yAT","\n";
  150. print "C|C: $yCC", "\t", "C|T: $yCT","\n";
  151. print "G|C: $yGC", "\t", "G|T: $yGT","\n";
  152. print "T|C: $yTC", "\t", "T|T: $yTT","\n\n";
  153.  
Coloreado en 0.008 segundos, usando GeSHi 1.0.8.4


Como podéis comprobar, separo los primeros nucleótidos de cada inicio de bloque al no estar precedidos por ninguno.

Y adjunto un par de ejemplos de secuencias a analizar ($último nt que no es CpG, $último nt que es CpG, $secuencia)

Quizá pida demasiado. En todo caso toda ayuda que me podáis proporcionar, ya sean pautas que he de seguir, pistas,
consejos, etc, ya me vendrán muy bien, así que lo que podáis.

Muchas gracias.
Adjuntos
code_29665_1.txt
Secuencias a analizar
(5.14 KiB) 100 veces
zipetardix
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2012-02-04 16:27 @727
Arriba

Publicidad

Re: Modelos ocultos de Márkov

Notapor explorer » 2012-02-11 07:52 @369

zipetardix escribiste:Ahora en el siguiente paso tengo que implementar un algoritmo de búsqueda en otro script en el que, leyendo las probabilidades calculadas en el anterior (asignadas tanto para "isla CpG", como para no), calcule la puntuación para islas CpG de una longitud dada a lo largo de la secuencia. Así pues, el script debe aceptar el input de la longitud de dicha "isla", que lo escribe el usuario.
Yo no dispongo de suficiente información como para saber cómo ha de ser ese algoritmo, pero viendo el código mostrado, si te das cuenta, hay muchas líneas que son prácticamente iguales. Eso quiere decir que existe otra forma de hacer este primer programa, más corta.

Además, salen errores:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
Use of uninitialized value $not_cpg[1000] in hash element at code_29665.pl line 34, <$FILE1> line 1.
Use of uninitialized value within @not_cpg in hash element at code_29665.pl line 34, <$FILE1> line 1.
Use of uninitialized value $not_cpg[1000] in hash element at code_29665.pl line 34, <$FILE1> line 1.
Use of uninitialized value within @yes_cpg in hash element at code_29665.pl line 39, <$FILE1> line 1.
Use of uninitialized value $not_cpg[1000] in hash element at code_29665.pl line 34, <$FILE1> line 2.
Use of uninitialized value within @not_cpg in hash element at code_29665.pl line 34, <$FILE1> line 2.
Use of uninitialized value $not_cpg[1000] in hash element at code_29665.pl line 34, <$FILE1> line 2.
Use of uninitialized value within @yes_cpg in hash element at code_29665.pl line 39, <$FILE1> line 2.
Use of uninitialized value $not_cpg[1000] in hash element at code_29665.pl line 34, <$FILE1> line 3.
Use of uninitialized value within @not_cpg in hash element at code_29665.pl line 34, <$FILE1> line 3.
Use of uninitialized value $not_cpg[1000] in hash element at code_29665.pl line 34, <$FILE1> line 3.
Use of uninitialized value within @yes_cpg in hash element at code_29665.pl line 39, <$FILE1> line 3.
Use of uninitialized value $not_cpg[1000] in hash element at code_29665.pl line 34, <$FILE1> line 4.
Use of uninitialized value within @not_cpg in hash element at code_29665.pl line 34, <$FILE1> line 4.
Use of uninitialized value $not_cpg[1000] in hash element at code_29665.pl line 34, <$FILE1> line 4.
Use of uninitialized value within @yes_cpg in hash element at code_29665.pl line 39, <$FILE1> line 4.
Use of uninitialized value in addition (+) at code_29665.pl line 84, <$FILE1> line 4.
Use of uninitialized value in addition (+) at code_29665.pl line 85, <$FILE1> line 4.
Use of uninitialized value in addition (+) at code_29665.pl line 87, <$FILE1> line 4.
...
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4

El problema está en la línea 29:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.     my @not_cpg = @nt_seq[0..499, $second_col..$len];
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

El valor de $len es absoluto (la longitud de @nt_seq), pero ahí lo estás usando como índice para extraer los valores de dentro de @nt_seq, así que, en realidad, estás metiendo un valor demás (e indefinido) en @not_cpg, saltando el error en la línea 34.

Luego vienen los tres bucles for de las líneas 44, 50 y 56. ¿Es correcto que solo den una vuelta cada uno de ellos? Fíjate, por ejemplo en el primero, que @firstnt_nocpg_first siempre contiene un solo valor: $nt_seq[0]. Y lo mismo para los otros dos... ¿eso es correcto? ¿Hacer un bucle de una vuelta para sumar un solo valor?

Sería interesante saber una descripción del algoritmo que quieres implementar aquí, para el cálculo de frecuencias. Sé que cuentas las veces que aparece un nucleótido después de otro, pero a partir de ahí no veo muy claro lo que quieres hacer...
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: Modelos ocultos de Márkov

Notapor zipetardix » 2012-02-11 19:04 @836

Hola explorer,

Se trata de un programilla que trata de analizar las secuencias dadas del siguiente modo:

- En la secuencia hay tres subtipos: una "isla CpG", que empieza y termina en los nucleótidos especificados por la primera y segunda columna respectivamente y a ambos flancos secuencia "background" o "no CpG" de 500nt cada una.

- Hay que analizar las frecuencias en las que un nt da paso a otro en concreto tanto en un subtipo como en el otro.

- A parte, debido a la estructura del hash, hay que contar a parte los nt que comienzan dichos subtipos, ya que estos no están precedidos por ninguno. En este caso:

- - @firstnt_nocpg_first = $nt_seq[0]; me construye un array de todos los primeros nt de la primera secuencia flanqueante.
- - @firstnt_nocpg_second = $nt_seq[$second_col]; lo mismo pero con los primeros nt de la segunda secuencia flanqueante. Lo normal es hacer ambos juntos puesto que luego las sumo pero inexplicablemente me duplicaba los resultados así que lo he dejado así.
- - @firstnt_yescpg = $nt_seq[$cpg_start]; primer nucleótido de cada isla CpG.

- Para el siguiente paso, siguiendo el ejemplo expuesto aquí, había pensado en hacer un output de una matriz para después poder trabajar con ella para calcular scores en otras secuencias. Eso lo he conseguido hacer pero ¿podría hacerse sin necesidad de escribir un archivo? Y ¿cómo podría trabajar con dicha matriz donde tengo primer nucleótido en filas y segundo nucleótido en columnas de manera que al leer la secuencia problema, el programa sepa recorrer la matriz hasta el resultado adecuado?


¡Muchas gracias por la ayuda hasta ahora y por la que pueda venir!
zipetardix
Perlero nuevo
Perlero nuevo
 
Mensajes: 5
Registrado: 2012-02-04 16:27 @727
Arriba

Re: Modelos ocultos de Márkov

Notapor explorer » 2012-02-11 21:09 @923

Entonces, las líneas 43 a 59 se pueden reducir a:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.     $first_nocpg_first { $nt_seq[          0] }++;
  2.     $first_nocpg_second{ $nt_seq[$second_col] }++;
  3.     $first_yescpg      { $nt_seq[ $cpg_start] }++;
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4
Es decir, van sumando las apariciones de los nt en esas posiciones.
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: Modelos ocultos de Márkov

Notapor explorer » 2012-02-14 18:27 @810

Esta es mi versión, pero no ha pasado los controles de calidad (me refiero a que no sabemos si realiza bien los cálculos) ya que no sabemos a dónde pertenecen los nucleótidos indicados por las dos primeras columnas de datos.
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. #
  3. # Estadísticas de frecuencias de aparición entre dos nucleótidos consecutivos, en dos zonas diferentes.
  4. #
  5. use Modern::Perl;                       # Somos modernos
  6. use utf8::all;                          # Activar todo el soporte UTF-8
  7. use List::Util 'sum';                   # Operaciones de lista de elementos
  8.  
  9.  
  10. ### Constantes
  11. my @NT = qw( A C G T );
  12.  
  13.  
  14. ### Variables
  15. my %nt_seq_no_CpG;                      # background seq
  16. my %nt_seq_si_CpG;                      # CpG island seq
  17. my %firstnt_no_CpG_1;                   # 1st nucleotide of the 1st background seq
  18. my %firstnt_no_CpG_2;                   # 1st nucleotide of the 2nd background seq
  19. my %firstnt_si_CpG;                     # 1st nucleotide of the CpG island
  20.  
  21.  
  22. ### Proceso
  23. # Recibimos el archivo a procesar por la entrada estándar, o
  24. # por el nombre del archivo a procesar, pasado como argumento
  25. while (<>) {
  26.     chomp;
  27.  
  28.     if (3 != (my($start, $stop, $seq) = split)) {       # Cada línea debe tener 3 campos
  29.         say "ERROR en línea $.";
  30.         next;
  31.     }
  32.     else {
  33.         my @seq = split //, $seq;
  34.         unshift @seq, ' ';                              # Insertamos nt en posición 0,
  35.                                                         # para que todos los cálculos estén basados en 1
  36.  
  37.         my @no_CpG = @seq;                                      # La sub secuencia no CpG es igual a toda la secuencia,
  38.         my @si_CpG = splice @no_CpG, $start, $stop - $start;    # menos la parte de la isla
  39.  
  40.         ## Frecuencia de nt consecutivos, fuera de la isla
  41.         for (my $i = 0; $i < $#no_CpG; $i++) {
  42.             $nt_seq_no_CpG{ $no_CpG[$i] }{ $no_CpG[$i+1] }++;
  43.         }
  44.  
  45.         ## Frecuencia de nt consecutivos, en la isla
  46.         for (my $i = 0; $i < $#si_CpG; $i++) {
  47.             $nt_seq_si_CpG{ $si_CpG[$i] }{ $si_CpG[$i+1] }++;
  48.         }
  49.  
  50.         $firstnt_no_CpG_1{ $seq[1     ] }++;
  51.         $firstnt_si_CpG  { $seq[$start] }++;
  52.         $firstnt_no_CpG_2{ $seq[$stop ] }++;
  53.     }
  54. }
  55.  
  56. # Cálculo de las estadísticas fuera de la isla
  57. my %no_freq;
  58. for my $nt1 (@NT) {
  59.     my $suma_nt1 = sum map { $nt_seq_no_CpG{$_}{$nt1} } @NT;
  60.  
  61.     for my $nt2 (@NT) {
  62.         $no_freq{$nt2}{$nt1} = $nt_seq_no_CpG{$nt2}{$nt1} / $suma_nt1;
  63.     }
  64. }
  65.  
  66. my %pre_no_freq;
  67. for my $nt (@NT) {
  68.     $pre_no_freq{$nt}  = $firstnt_no_CpG_1{$nt} // 0;
  69.     $pre_no_freq{$nt} += $firstnt_no_CpG_2{$nt} // 0;
  70. }
  71.  
  72. my $suma_pre_no_freq = sum map { $pre_no_freq{$_} } @NT;
  73.  
  74. #use Data::Dumper::Simple;
  75. #say Dumper %firstnt_no_CpG_1;
  76. #say Dumper %firstnt_no_CpG_2;
  77. #say Dumper %pre_no_freq;
  78. #say $suma_pre_no_freq;
  79.  
  80. for my $nt (@NT) {
  81.     $pre_no_freq{$nt} /= $suma_pre_no_freq;
  82. }
  83.  
  84. say 'Frecuencias fuera de la isla CpG:';
  85. for my $nt (@NT) {
  86.     say "$nt: $pre_no_freq{$nt}";
  87. }
  88.  
  89. for my $nt1 (@NT) {
  90.     for my $nt2 (@NT) {
  91.         printf "%s|%s: %3.3f\n", $nt2, $nt1, $no_freq{$nt2}{$nt1};
  92.     }
  93. }
  94.  
  95. # Cálculo de las estadísticas en la isla
  96. my %si_freq;
  97. for my $nt1 (@NT) {
  98.     my $suma_nt1 = sum map { $nt_seq_si_CpG{$_}{$nt1} } @NT;
  99.  
  100.     for my $nt2 (@NT) {
  101.         $si_freq{$nt2}{$nt1} = $nt_seq_si_CpG{$nt2}{$nt1} / $suma_nt1;
  102.     }
  103. }
  104.  
  105. my %pre_si_freq;
  106. for my $nt (@NT) {
  107.     $pre_si_freq{$nt} = $firstnt_si_CpG{$nt} // 0;
  108. }
  109.  
  110. my $suma_pre_si_freq = sum map { $pre_si_freq{$_} } @NT;
  111.  
  112. for my $nt (@NT) {
  113.     $pre_si_freq{$nt} /= $suma_pre_si_freq;
  114. }
  115.  
  116. say 'Frecuencias en la isla CpG:';
  117. for my $nt (@NT) {
  118.     say "$nt: $pre_si_freq{$nt}";
  119. }
  120.  
  121. for my $nt1 (@NT) {
  122.     for my $nt2 (@NT) {
  123.         printf "%s|%s: %3.3F\n", $nt2, $nt1, $si_freq{$nt2}{$nt1};
  124.     }
  125. }
  126.  
  127. __END__
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4

Se podría reducir aún más, ya que los cálculos estadísticos son casi los mismos tanto fuera como dentro de la isla.
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
Anterior
Publicar una respuesta

Volver a Bioinformática

¿Quién está conectado?

Usuarios navegando por este Foro: No hay usuarios registrados visitando el Foro y 0 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