• Publicidad

Modelos ocultos de Márkov

Perl aplicado a la bioinformática

Modelos ocultos de Márkov

Notapor oshaoran » 2012-01-24 12:19 @555

¡¡Hola!! Soy estudiante de ing. biológica, y actualmente curso la asignatura de bioinformática, debía crear una rutina para que dada la secuencia de un péptido el programa me identifique si dicho péptido es o no un péptido antimicrobiano.

La idea general es: Los péptidos antimicrobianos de mi trabajo, son péptidos con características especificas que los diferencian de los demás péptidos. Tienen varios aminoácidos que se conservan en posiciones muy especificas (motivos). Entonces, usando el modelo estadístico de cadenas ocultas de Márkov, yo le asigno una probabilidad de ocurrencia a cada aminoácido en cada posición de la secuencia.

Dada una secuencia, el programa debe tomarla reconociendo las posiciones, e ir leyendo posición por posición, leer el aminoácido de la secuencia, y compararlo con las probabilidades de que realmente ese aminoácido esté en esa posición (probabilidades tomadas de una matriz (aminoácido X posición), tomar el dato de probabilidad de cada posición y sumarlo, y al final entregar la sumatoria de dichas probabilidades.

Lo anterior debe repetirse para un número X de secuencias que el programa debe leer de un archivo en formato fasta. Luego, usando la suma de cada secuencia, arrojar finalmente como salida, las secuencias cuya probabilidad se encuentre en un intervalo suministrado y así puedo concluir que dichas secuencias son o no un péptido antimicrobiano.

Ya es complicado decirlo en palabras, ¿¿no?? ¡je,je,je!, bueno mis conocimientos son bastante básicos en el tema. Inicialmente la rutina que quería armar, era haciendo un ciclo for() que me leyera cada posición y de acuerdo al aminoácido que encontrara, asignara el valor que requiero tomándolo de una matriz diseñada en Excel u otro formato, pero no tengo idea de cómo implementar la matriz con los valores a asignar :cry: . Luego encontré una rutina ya hecha que ejecuta el proceso de las cadenas ocultas de Márkov, pero aplicado al procesamiento de textos.
Espero, puedan ayudarme a que esa rutina, en vez de identificarme palabras en un texto, me identifique los aminoácidos que quiero.

Llevo más de dos semanas intentando de todo, tratando de hacer la rutina yo mismo, y luego tratando de editar la que encontré, pero nada me funciona :(

Les agradezco cualquier ayuda o explicación que me puedan brindar Uds. que son bastante conocedores sobre programación.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl -w
  2.  
  3.  
  4. use strict;
  5.  
  6. my @words;    # words on a line
  7. my %wordlist; # key: prefix, value: anon hash (k: suffix,
  8.               #                                v: frequency)
  9.  
  10. my $pref_len = shift @ARGV || 2;
  11. my $maxwords = shift @ARGV || 100;
  12.  
  13. my $entries  = 0;
  14.  
  15. # build word list
  16. #
  17. # 'Blessed is the man that walketh not in the counsel'
  18. # %wordlist = ( 'blessed is' => { 'the' => 1, },
  19. #               'is the'     => { 'man' => 1, },
  20. #               'the man'    => { 'that'=> 1, },
  21. #             );
  22. #
  23. while (<>) {
  24.   my $suf;
  25.  
  26.   push @words, split;
  27.  
  28.   while ( @words > $pref_len )  {
  29.     # build prefix of $pref_len words
  30.     # join(' ', @array) is faster than qq(@array) or "@array"
  31.     #
  32.     my $pref = join(' ', @words[0..($pref_len-1)]);
  33.  
  34.     # add suffix to list
  35.     #
  36.     $suf = $words[$pref_len];
  37.  
  38.     $wordlist{$pref}{$suf}++;
  39.  
  40.     shift @words; # next word on this line
  41.  
  42.     $entries++;
  43.   }
  44. }
  45.  
  46. # change frequency count to a percentage
  47. # (with help from pcb, recipe 2.10)
  48. #
  49. foreach my $href ( values %wordlist ) {
  50.   foreach ( values %$href ) {
  51.     $_ /= $entries;
  52.   }
  53. }
  54.  
  55. # starting point
  56. #
  57. my $pref = (keys %wordlist)[rand keys %wordlist];
  58.  
  59. print "$pref";
  60.  
  61. # dump out listings
  62. #
  63. for (0..($maxwords-1)) {
  64.   last unless (exists $wordlist{$pref});
  65.  
  66.   my $suf = weighted_suffix();
  67.  
  68.   print ' '. $suf;
  69.  
  70.   print "\n" if ( $_ % 10 == 0);
  71.  
  72.   # skip past first word in prefix
  73.   #
  74.   $pref =~ s/^[^ ]+ (.+)$/$1 $suf/;
  75. }
  76.  
  77. exit;
  78.  
  79. # from pcb (recipe 2.10)
  80. #
  81. sub weighted_suffix {
  82.   my ($suf,$weight,$rand);
  83.  
  84.   while (1) {
  85.     $rand = rand;
  86.  
  87.     while ( ($suf,$weight) = each %{ $wordlist{$pref} } ) {
  88.       return $suf if ($rand -= $weight) < 0;
  89.     }
  90.   }
  91. }
Coloreado en 0.007 segundos, usando GeSHi 1.0.8.4
Avatar de Usuario
oshaoran
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2011-09-14 09:12 @425

Publicidad

Re: Modelos ocultos de Márkov

Notapor pvaldes » 2012-01-24 18:12 @800

No pillo la idea... Vas a tener que explicarte mejor porque nos falta mucha información, el código que pones cortapegado no me sugiere nada, tampoco

Lo primero de todo elimina de tu código todo lo innecesario o no relacionado con tu problema concreto, usa nombres de variable apropiados a su contenido y depura los comentarios eliminando todo lo innecesario.

El programa que pones corre con dos argumentos, ¿qué esperas pasarle dentro de estos argumentos?

$_ /= $entries;

Probablemente esto no hace lo que esperas que haga.

Vamos a suponer que tenemos la siguiente secuencia de aminoácidos:

Valina-Leucina-Treonina-Lisina-Lisina-Leucina

¿cuál es el resultado que esperas?

> tomándolo de una matriz diseñada en Excel

Tu programa no conecta con Excel. De todos modos puedes usar una variable hash para eso...
pvaldes
Perlero nuevo
Perlero nuevo
 
Mensajes: 129
Registrado: 2011-01-22 12:56 @580

Re: Modelos ocultos de Márkov

Notapor oshaoran » 2012-01-24 20:06 @879

Ok, quizá subir ese código no fue la mejor idea; solo trataba de mostrar cómo es que puede ejecutarse lo de los modelos de Márkov.

Bueno, trataré de explicarme mejor... Tengo la siguiente secuencias de letras, por ejemplo:

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
1.ASDFGHJK
2.SDFGAHKJ
3.HGFDSKAS
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Yo tengo en Excel una matriz que me dice de acuerdo a un análisis estadístico ya hecho, cuál es la probabilidad de que cada letra, realmente se encuentre en dicha posición; es algo así, por ejemplo (Letra vs. posición)
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
        1       2       3       4       5       6       7       8      
A       0.1     0.4     0.9     0.1     0.1     0.1     0.4     0.7    
S       0.5     0.5     0.5     0.5     0.5     0.5     0.5     0.5    
D       0.7     0.7     0.7     0.7     0.2     0.7     0.9     0.7    
F       0.7     0.7     0.5     0.7     0.7     0.7     0.4     0.7    
G       0.6     0.6     0.6     0.6     0.6     0.6     0.6     0.6    
H       0.4     0.4     0.4     0.4     0.4     0.4     0.4     0.4    
J       0.3     0.3     0.3     0.2     0.3     0.4     0.3     0.1    
K       0.5     0.4     0.8     0.7     0.4     0.5     0.9     0.7
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


El programa debería leerme cada secuencia, tomar cada probabilidad y sumarla
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
ASDFGHJK: 0.1+ 0.5+0.7+0.7+0.6+0.4+0.3+0.3=3.6
SDFGAHKJ: 0.5+0.7+0.5+0.6+0.1+0.4+0.9+0.1=3.8
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


El fundamento es que con un orden distinto el puntaje final va ser distinto también... entonces yo tengo un intervalo, por ejemplo (3.0-3.6), entonces el programa debería decirme: la secuencia 1 está dentro del intervalo, y con eso yo concluyo que dicha secuencia es un péptido antimicrobiano.

Mi problema en esencia es que no sé cómo programar esa parte de la rutina para que el programa me lea la secuencia y pueda realizar las sumas tomando la matriz. Mis conocimientos son demasiado básicos. El único programa que he logrado hacer fue uno para hallar ORF, y solo usé if() y for(); no sabría cómo usar ese hash que tu me propones.

Lo poco que he podido programar fue esto:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. ###########################
  3. #nombre de secuencia de entrada (aminoacidos)
  4. $secuencia= 'secprueba.fasta'; 
  5. open (SECUENCIA, $secuencia);
  6. @secuencia= <SECUENCIA>;                                                                                               
  7. close SECUENCIA;
  8. print @secuencia;
  9.  
  10. ####################################
  11. #$seccy=join ("",@cyclotide);                                                                          
  12. #print @cyclotide;
  13. #$seccy=~s/\n//g;
  14. #print length $seccy, "\n";
  15. #################################
  16. $secuencia[o]="";
  17. $proteina=join("",@secuencia);
  18. $proteina=~ s/\n//g;
  19. $test= substr($proteina, 0,21);
  20.  
  21. for ($j=1; $j< scalar @secuencia; $j= $j+2) {
  22. print "la linea $j es:\n";
  23. $seq = $secuencia [$j];
  24.  
  25. print $seq, "\n";
  26.  
  27. @test=split('',$seq);
  28.  
  29.         for ($i=0; $i< scalar @test; ++$i) {
  30.         if (substr($test,$i,1) eq 'G') {$pos=$pos+1.8;}
  31.         elsif (substr($test,$i,1) eq 'C') {$pos=$pos+2.5;}
  32.         elsif (substr($test,$i,1) eq 'V') {$pos=$pos+3.5;}
  33.         elsif (substr($test,$i,1) eq 'S') {$pos=$pos+3.5;}
  34.         #
  35.         #
  36.         #
  37.         #
  38.         #un elsif para cada aminoacido
  39.         #
  40.        
  41.         else {print "error posicion $i\n";}
  42.         }
  43. print $pos, "\n";
  44.  
  45.         }      
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4


Pero lo que hace es sumarme el valor que le doy cuando encuentra la letra, sin importar la posición. Lo que yo deseo es, como he dicho, que si encuentra por ejemplo, S en la primera posición, sume un valor, pero si en vez de S en dicha posición encuentra V, entonces sume un valor distinto, así para cada posición.
Avatar de Usuario
oshaoran
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2011-09-14 09:12 @425

Re: Modelos ocultos de Márkov

Notapor explorer » 2012-01-24 23:40 @028

Bienvenido a los foros de Perl en español, oshaoran.

Aquí tienes un ejemplo de lo que se puede hacer.

Supongamos que tenemos un archivo FASTA así:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
>98448|fgenesh3_pg.4971__1
JHJJDJKA
>23766|gw1.611.1.1
SJJKGFSD
>28683|gw1.3814.2.1
KKJAJGDH
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
y un archivo con las estadísticas, así:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
        1       2       3       4       5       6       7       8      
A       0.1     0.4     0.9     0.1     0.1     0.1     0.4     0.7    
S       0.5     0.5     0.5     0.5     0.5     0.5     0.5     0.5    
D       0.7     0.7     0.7     0.7     0.2     0.7     0.9     0.7    
F       0.7     0.7     0.5     0.7     0.7     0.7     0.4     0.7    
G       0.6     0.6     0.6     0.6     0.6     0.6     0.6     0.6    
H       0.4     0.4     0.4     0.4     0.4     0.4     0.4     0.4    
J       0.3     0.3     0.3     0.2     0.3     0.4     0.3     0.1    
K       0.5     0.4     0.8     0.7     0.4     0.5     0.9     0.7
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
Entonces, con este programa:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use autodie;                    # Mejor morir que regresar con deshonor --Proverbio Klingon
  3. use Data::Dumper;               # Para ver las estructuras complejas
  4.  
  5. ### Constantes
  6. my @intervalo = (3.0, 3.6);
  7.  
  8. ### Variables
  9. my %estadisticas;               # Estructura bidimensional.
  10.                                 # Cada fila está indexada por una clave, que será una letra.
  11.                                 # El valor correspondiente es una referencia a un array, que
  12.                                 # almacenará los distintos valores de probabilidad según la
  13.                                 # posición dentro de la secuencia. Ejemplo
  14.                                 #
  15.                                 # $estadisticas{'A'}->[3] : probabilidad de 'A' cuando está
  16.                                 #                           en cuarta posición
  17.  
  18.  
  19.  
  20. ### Lectura de las estadísticas
  21. open my $excel, '<', 'excel.txt';
  22. <$excel>;                               # leemos y despreciamos la primera línea, la de los números
  23. while (<$excel>) {
  24.     chomp;                              # quitamos el fin de línea
  25.     my @columnas = split;               # separamos por los espacios en blanco
  26.  
  27.                                         # Construimos la matriz
  28.     my $letra = shift @columnas;        # La primera columna es la de las letras
  29.     $estadisticas{$letra} = \@columnas; # el resto, las probabilidades
  30. }
  31. close $excel;
  32.  
  33. #print Dumper(\%estadisticas);          # Vemos el aspecto de nuestra criatura
  34.  
  35. ### Lectura del fichero fasta
  36. my $nombre_secuencia;
  37. my $secuencia;
  38.  
  39. open my $fasta, '<', 'fasta.txt';
  40. while (my $linea = <$fasta>) {
  41.     chomp $linea;
  42.    
  43.     if ($linea =~ /^>(.+)/) {           # Si es una línea de cabecera
  44.  
  45.         if ($secuencia) {               # Si teníamos una secuencia anterior
  46.             procesa($secuencia, $nombre_secuencia);     # la procesamos
  47.         }
  48.  
  49.         $nombre_secuencia = $1;         # nos quedamos con el nombre de la nueva secuencia
  50.         $secuencia = '';                # y reiniciamos la secuencia a leer
  51.     }
  52.     else {                              # sino, es parte de una secuencia
  53.         $secuencia .= $linea;           # la agregamos como secuencia que estamos leyendo
  54.     }
  55. }
  56.  
  57. if ($secuencia) {                       # en caso de llegar al final del fichero fasta
  58.     procesa($secuencia, $nombre_secuencia);     # procesamos la última secuencia
  59. }
  60.  
  61. close $fasta;
  62.  
  63. sub procesa {
  64.     my ($seq, $nombre) = @_;            # leemos los argumentos
  65.  
  66.     my $total = 0;                      # total acumulado
  67.     my $pos;                            # posición de la letra en la secuencia
  68.  
  69.     for ($pos = 0; $pos < length $seq; $pos++) {        # para toda la secuencia
  70.  
  71.         my $letra = substr $seq, $pos, 1;               # sacamos una letra cada vez
  72.  
  73.         if (not $estadisticas{$letra}) {                # comprobación
  74.             die "ERROR: No tenemos estadísticas para letra $letra\n";
  75.         }
  76.  
  77.         if (not exists $estadisticas{$letra}->[$pos]) { # más comprobaciones
  78.             die "ERROR: No tenemos estadísticas para letra $letra en posición ", $pos+1, "\n";
  79.         }
  80.  
  81.         $total += $estadisticas{$letra}->[$pos];        # acumulamos
  82.     }
  83.  
  84.     if ($total >= $intervalo[0]  and  $total <= $intervalo[1]) {
  85.         print "$nombre_secuencia está dentro del intervalo: $total\n";
  86.     }
  87.     else {
  88.         print "$nombre_secuencia no está dentro del intervalo: $total\n";
  89.     }
  90. }
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4
tenemos la salida
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
98448|fgenesh3_pg.4971__1 está dentro del intervalo: 3.4
23766|gw1.611.1.1 no está dentro del intervalo: 4.3
28683|gw1.3814.2.1 está dentro del intervalo: 3.5
Coloreado en 0.000 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: Modelos ocultos de Márkov

Notapor pvaldes » 2012-01-25 05:14 @260

Me daba la impresión, de que había un problema con eso, sí.

La probabilidad de tener ASDFGHJK no es 0.1+ 0.5+0.7+0.7+0.6+0.4+0.3+0.3 = 3.6

Probabilidad de tener A + S = probabilidad de tener A * probabilidad de tener S ... etc, luego

0.1*0.5*0.7*0.7*0.6*0.4*0.3*0.3 = 0.0005292

Si no, podríamos alcanzar probabilidades superiores al 100% fácilmente, (en éste ejemplo el 360%) lo cual no tiene sentido estadístico.
pvaldes
Perlero nuevo
Perlero nuevo
 
Mensajes: 129
Registrado: 2011-01-22 12:56 @580

Re: Modelos ocultos de Márkov

Notapor explorer » 2012-01-25 06:49 @326

Pero es que, justamente, el modelo oculto de Márkov, es la suma de las probabilidades de ocurrencia de un determinado evento cuando ocurre otro evento relacionado en un determinado momento o espacio. Las probabilidades ya las ha calculado antes, con la Excel.

Lo interesante para mi sería hacer el cálculo de probabilidades en el mismo programa.

Si oshaoran nos da más detalles del cálculo, lo podríamos intentar resolver.

Y me parece que el ejemplo que puso de la cadena ASDFGHJK, no hace bien la cuenta... el último 0.3 debería ser 0.7, creo.
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: Modelos ocultos de Márkov

Notapor oshaoran » 2012-01-26 22:09 @964

¡¡Hey!!! :D Les agradezco muchísimo, realmente la rutina me fue muy útil, pude completar a la perfección mi trabajo, nuevamente muchas gracias, realmente valoro que aun existan personas, dispuestas a compartir su conocimiento.

Bueno, continuando con ello, lo de hallar las probabilidades, quizá sea algo más complicado. Para hallarlas, tuvimos que hacer un alineamiento de secuencias de unos 30 péptidos antimicrobianos, para descubrir qué posiciones se conservan, y poder sacar las probabilidades. Fue algo como esto:

Imagen

Primero sacamos un promedio de las bases en cada posición, para el alineamiento mostrado sería P(S)=2/16, P(G)= 2/16... En nuestro modelo obviamos los gaps (simbolizados por un guion) pero pues también se podrían incluir: P(gap)= 12/16. Eso sería para posición 1... y así para todas las posiciones sacando la probabilidad de cada aminoácido.

Con solo el promedio no podríamos montar un algoritmo correcto, porque no estaríamos teniendo en cuenta las cadenas de Márkov, que para este caso, como conocemos unas posiciones que se conservan (mírese las C, que son constantes para dichos péptidos), podremos dilucidar por Márkov cuál sería el aminoácido que debería ir antes (esto fundamentado con cuestiones de polaridad de aminoácidos, conformación de proteínas, etc, etc).

Para hallar dichas probabilidades de Márkov, usamos una formula matemática, muy aproximada, demasiado creo yo ¡je,je!... que era sacar el promedio, luego a todos los valores sumarles un valor arbitrario (para evitar que existan ceros) pues luego al valor promedio le sacábamos el logaritmo. En esencia, llenamos nuestra matriz de probabilidades con valores de este calculo. En resumen fue algo así... :lol: Nuestro trabajo era solo mostrar un prototipo, una aproximación, de ahí que no desarrollamos más lo del modelo Márkov (que, honestamente, para mi es algo complejo).

Siguiendo tu idea pues, el programa debería entonces leer las posiciones de todas las secuencias del alineamiento, ejecutar la fórmula matemática y llenar la matriz que se usaría para los cálculos posteriores...
Avatar de Usuario
oshaoran
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2011-09-14 09:12 @425

Re: Modelos ocultos de Márkov

Notapor explorer » 2012-01-27 22:17 @970

Una pregunta... ¿los espacios (gap) se pusieron a mayores? es decir, no estaban en las secuencias originales, ¿verdad?
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: Modelos ocultos de Márkov

Notapor oshaoran » 2012-01-28 02:04 @128

Nop, los gaps los genera el programa de alineamiento; el archivo con las secuencias en FASTA no los tiene.
Avatar de Usuario
oshaoran
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2011-09-14 09:12 @425

Re: Modelos ocultos de Márkov

Notapor explorer » 2012-01-28 19:43 @863

Todo el problema que has comentado en este hilo, parece extraído del libro "Genomic Perl: From Bioinformatics Basics to Working Code".

Quiero decir que en ese libro se comentan y publican los códigos para resolver este problema. Desde el alineamiento de múltiples secuencias usando variantes del algoritmo de Needleman-Wunsch, hasta la construcción de las matrices PAM (Porcentaje de Mutación Aceptada).

Y todo en código Perl :)
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

Siguiente

Volver a Bioinformática

¿Quién está conectado?

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