Código en C dentro de un programa Perl

Hola, siguiendo el hilo de los comentarios a la entrada anterior de Álvaro, hoy pongo un ejemplo de como incluir una o más subrutinas escritas en C dentro de un texto en Perl, recurriendo al módulo Inline:

 #!/usr/bin/perl -w 
 # programa 7 
 use Inline C; 
 use strict; 
 my $BASE = 0.987654321; 
 my $EXP = 1000; 
 print "# $BASE^$EXP = \n"; 
 print "# ".potencia_perl($BASE,$EXP)." (perl)\n"; 
 print "# ".potencia_C($BASE,$EXP)." (C)\n"; 
 sub potencia_perl 
 { 
    my ($base,$exp) = @_; 
    my $result = $base; 
    for(my $p=1; $p<$exp; $p++) { $result *= $base; } 
    return $result; 
 } 
 __END__  
 ### alternativa en C 
 __C__ 
 #include <stdio.h> 
 float potencia_C( double base, double exp )  
 { 
    double result = base; 
    int p; 
    for(p=1; p<exp; p++) 
    {  
       result *= base; 
    }    
    return result; 
 } 

La salida obtenida pone en evidencia que diferentes lenguajes suelen tener diferentes precisiones a la hora de representar números y operar con ellos:

$ perl prog7.pl 
# 0.987654321^1000 = 
#  4.02687472213071e-06 (perl)
#  4.02687464884366e-06 (C)

Generar todas las posibles combinaciones posibles de n nucleótidos o aminoácidos

Erróneamente decimos que queremos generar todas las posibles 'combinaciones' de n letras, nucleótidos o aminoácidos cuando el término correcto sería 'variaciones con repetición' (Permutations-Variations-Combinations). Por ejemplo, las 16 posibles variaciones con repetición de 2 nucleótidos serían: AA, AC, AG, AT, CA, CC, CG, CT, GA, GC, GG, GT, TA, TC, TG, TT. El número de variaciones con repetición de n nucleótidos es 4^n y de n aminoácidos 20^n.

El siguiente código, es una adaptación de un código en PHP para generar todas las variaciones con repetición de n elementos, en nuestro caso, 2 aminoácidos y 4 nucleótidos.

 #!/usr/bin/perl -w  
 # Generate all the variations with repetition of n letters  
 my @aminoacids = ('A','R','N','D','C','Q','E','G','H','L','I','K','M','F','P','S','T','W','Y','V');  
 my @nucleotides = ('A','C','G','T');  
 # Generate all the variations of 2 aminoacids  
 my $aas = variations(\@aminoacids,2);  
 # Generate all the variations of 4 nucleotides  
 my $nts = variations(\@nucleotides,4);  
 # Print the results  
 print "\nAll the variations with repetition of 2 aminoacids: ";  
 foreach my $vari (@{$aas}){  
      foreach my $elem (@{$vari}){  
           print "$elem";  
      }  
      print " ";  
 }  
 print "\n";  
 print "\nAll the variations with repetition of 4 nucleotides: ";  
 foreach my $vari (@{$nts}){  
      foreach my $elem (@{$vari}){  
           print "$elem";  
      }  
      print " ";  
 }  
 print "\n\n";  
 sub variations {  
      my ($letters,$num) = @_;  
      my $last = [map $letters->[0] , 0 .. $num-1];  
      my $result;  
      while (join('',@{$last}) ne $letters->[$#{$letters}] x $num) {  
           push(@{$result},[@{$last}]);  
           $last = char_add($letters,$last,$num-1);  
           print '';  
      }  
      push(@{$result}, $last);  
      return $result;  
 }  
 sub char_add{  
      my ($digits,$string,$char) = @_;  
      if ($string->[$char] ne $digits->[$#{$digits}]){  
           my ($match) = grep { $digits->[$_] eq $string->[$char]} 0 .. $#{$digits};  
           $string->[$char] = $digits->[$match+1];  
           return $string;  
      } else {  
           $string = changeall($string,$digits->[0],$char);  
           return char_add($digits,$string,$char-1);  
      }  
 }  
 sub changeall {  
      my ($string,$char,$start,$end) = @_;  
      if (!defined($start)){$start=0;}  
      if (!defined($end)){$end=0;}  
      if ($end == 0) {$end = $#{$string};}  
      for(my $i=$start; $i<=$end; $i++){  
           $string->[$i] = $char;  
      }  
      return $string;  
 }  

Para terminar una breve lección de estadística, fuente: Aulafacil.com

a) Combinaciones:
Determina el número de subgrupos de 1, 2, 3, etc. elementos que se pueden formar con los "n" elementos de una nuestra. Cada subgrupo se diferencia del resto en los elementos que lo componen, sin que influya el orden.

Por ejemplo, calcular las posibles combinaciones de 2 elementos que se pueden formar con los números 1, 2 y 3.
Se pueden establecer 3 parejas diferentes: (1,2), (1,3) y (2,3). En el cálculo de combinaciones las parejas (1,2) y (2,1) se consideran idénticas, por lo que sólo se cuentan una vez.

b) Variaciones:
Calcula el número de subgrupos de 1, 2, 3, etc.elementos que se pueden establecer con los "n" elementos de una muestra. Cada subgrupo se diferencia del resto en los elementos que lo componen o en el orden de dichos elementos (es lo que le diferencia de las combinaciones).

Por ejemplo, calcular las posibles variaciones de 2 elementos que se pueden establecer con los número 1, 2 y 3.
Ahora tendríamos 6 posibles parejas: (1,2), (1,3), (2,1), (2,3), (3,1) y (3,3). En este caso los subgrupos (1,2) y (2,1) se consideran distintos.

c) Permutaciones:
Cálcula las posibles agrupaciones que se pueden establecer con todos los elementos de un grupo, por lo tanto, lo que diferencia a cada subgrupo del resto es el orden de los elementos.

Por ejemplo, calcular las posibles formas en que se pueden ordenar los número 1, 2 y 3.
Hay 6 posibles agrupaciones: (1, 2, 3), (1, 3, 2), (2, 1, 3), (2, 3, 1), (3, 1, 2) y (3, 2, 1)

Script para conocer el formato de un archivo

En el campo de la bioinformática se trabaja con numerosos formatos de archivos de datos biológicos (FASTA, GenBank, PDB...). Habitualmente tenemos que conocer el formato de los datos de un archivo que usamos como entrada en un script. Sin embargo, esto no es tarea fácil puesto que diferentes formatos pueden contener el mismo tipo de datos de forma poco estandarizada. En estos casos tenemos que usar nuestros conocimientos previos sobre el formato para establecer el tipo de archivo. El siguiente código de perl permite introducir información sobre los formatos más habituales en forma de expresiones regulares para posteriormente leer un fichero y asignar una puntuación a cada uno de los formatos definidos. De acuerdo a la puntuación final conoceremos el formato de fichero más probable.

 # Obtains the format of a file  
 sub obtain_file_format {  
      my ($file) = @_;  
      open(FILE,$file)|| die "# $0 : cannot read $file\n";  
      my %score;  
      while(<FILE>){  
           if (/^#/){  
                next;  
           }  
           if (/^>/){  
                $score{'fasta'}++;  
           }  
           if (/^(LOCUS|DEFINITION|ACCESSION|FEATURES)/){  
                $score{'genbank'}++;  
           }  
           if (/^(HEADER|TITLE|COMPND|SEQRES|ATOM)/){  
                $score{'pdb'}++;  
           }  
           if (/^(DE|VV|XX)/){  
                $score{'transfac'}++;  
           }  
      }  
      close(FILE);  
      my @sorted_scores;  
      if (@sorted_scores = sort { $score{$b} <=> $score{$a} } (keys %score)) {  
           return $sorted_scores[0];  
      } else {  
           return 'unknown';  
      }  
 }  

Bioinformática en la nube a partir de 10000 usuarios

La comunidad bioinformática tiene ya una cierta tradición en poner a disposición de los usuarios el software desarrollado a través de la red. Esto sirve para dar a conocer las últimas ideas y algoritmos que se van desarrollando en los laboratorios, de manera que se ponen a prueba en escenarios reales en manos de usuarios que no son sus creadores. Por supuesto, también sirve para ir ampliando el repertorio de herramientas que cualquiera con conexión a Internet puede usar, normalmente de forma gratuita.

Para el usuario final éstas son a menudo simplemente aplicaciones accesibles desde un navegador web, pero para sus creadores son programas más o menos complejos que viven en servidores web (como máquinas virtuales) que consumen luz y hay que administrar, que son a veces objeto de ataques informáticos, y que dependiendo del uso que tengan exigirán ampliaciones de hardware.

Qué coste real tienen estos servidores? El coste absoluto de mantener un servidor es difícil de promediar, porque dependerá de costes (luz, hardware, la carga de usuarios) y beneficios (evaluación peer to peer de los algoritmos, corrección de errores, retroalimentación de los usuarios, incluso citas bibliográficas) variables.
Otra cuestión  es si es mejor alojar estos servidores en nuestras propias máquinas o aprovecharnos de las arquitecturas de computación en la nube, que son escalables. Para ayudarnos a tomar esta decisión un reciente artículo en Bioinformatics del grupo de Cedric Notredame, donde prueban la versión paralelizada de su famoso alineador múltiple Tcoffee (basado en el principio de consistencia), sugiere que con los costes actuales sólo compensa alojar este tipo de servicios web en arquitecturas en nube a partir de 10000 usuarios al mes.

Nos vemos a la vuelta de las vacaciones.

Alineamiento perfil-perfil teniendo en cuenta la estructura secundaria

Con el alineamiento múltiple de secuencias de proteínas relacionadas se puede construir un perfil. A parte de su uso en técnicas de fold recognition o para predicción de estructura secundaria, este se puede utilizar para compararlo con las secuencias de otras proteínas e incluso con otro grupo de proteínas relacionadas entre sí, en este caso mediante una comparación perfil-perfil.

En el siguiente ejemplo, escrito en python y del que hay un ejemplo en perl aquí, vamos a realizar la comparación entre 2 perfiles, y apoyaremos el alineamiento con la información disponible de la estructura secundaria predicha para ambos. Utilizaremos un algoritmo de programación dinámica, por lo que si no se está familiarizado con la técnica se puede acudir a uno de los muchos libros sobre algoritmos o al excelente artículo de Sean R Eddy What is dynamic programming?, el cual se basa precisamente en ejemplos de alineamiento de secuencias biológicas.

Básicamente, necesitaremos un archivo PSSM y otro PSIPRED, perfil y estructura secundaria, respectivamente, para cada uno de los 2 grupos de proteínas. Puedes encontrar archivos de ejemplo en este curso, donde ya enlazamos para el ejemplo en perl.

Remarcar, entre los aspectos más importantes del código, que la puntuación entre perfiles se calcula como la media entre los valores de las 2 matrices y que se asigna una puntuación extra al caso de que 2 posiciones alineadas tengan el mismo estado de estructura secundaria.

 #!/usr/bin/python   
   
 from __future__ import with_statement  
 import re  
   
 __title__ = "Secondary structure and pairwise DP profile-profile alignment"  
 __author__ = "Carlos P Cantalapiedra"  
 __institution__ = "CSIC AulaDei Laboratory of Computational Biology"  
   
 ### MAIN STARTS UNDERNEATH THE FUNCTIONS ###  
   
 # Function that parses PSSM files  
 #    The format of the array is:  
 #    pssm[  
 #    {'aa':aa1, 'weights', [w1, w2, ...]},  
 #    {'aa':aa2, 'weights', [w1, w2, ...]},  
 #    .  
 #    .  
 #    .  
 #    ]  
 def read_PSSM(pssm_file):  
   pssm = [] # Returning dictionary  
   aa = ''  
   weights = []  
     
   # Finds the aminoacid (\w) and his weights ((\s*-?\d*.\d*)+)  
   reg_exp = re.compile(r"^\s*\d+\s+(\w)((\s*-?\d*.\d*)+)")  
     
   # NOTE: with the with_statement there's no need to close the file  
   with open(pssm_file, 'r') as file:  
     for line in file:  
       result = reg_exp.search(line)  
         
       if result:  
         try:  
           aa = result.group(1)  
           weights = str(result.group(2)).split()  
   
         except IndexError:  
           raise IOError(\  
           "There is a problem with the format of the PSSM file "+str(pssm_file))  
           
         pssm.append({'aa' : aa, 'weights' : weights})  
   return pssm  
   
 # Function that parses "Pred:" line of PSIPRED files,  
 #    to obtain the predicted secondary structures  
 def read_PSIPRED(psipred_file):  
   psipred = []  
     
   reg_exp = re.compile(r"^\s*Pred:\s*(\w+)")  
     
   # NOTE: with the with_statement there's no need to close the file  
   with open(psipred_file, 'r') as file:  
     for line in file:  
       result = reg_exp.search(line)  
         
       if result:  
         try:  
           [psipred.append(a) for a in result.group(1)]  
             
         except IndexError:  
           raise IOError(\  
           "There is a problem with the format of the PSIPRED file "+str(psipred_file))  
         
   return psipred  
   
 # Function that returns the bonus score associated to secondary structure  
 def get_secondary_score(st1, st2):  
   global secondary_score  
   score = 0  
     
   if st1 == st2 and st2 != 'C':  
     score = secondary_score  
   
   return score  
   
 # Function that creates de Dynamic Programming matrix  
 #     The format of the matrix is going to be:  
 #     matrix[  
 #       row(0) [('r', gap), ('r', gap), ..., ('r', gap)] # len(row) = len(prot2) + 1  
 #       row(1) [('c', gap), (path, score), ..., (path,score)]  
 #      row(2) [('c', gap), (path, score), ..., (path,score)]  
 #       .  
 #       .  
 #       .  
 #       row(len(prot1)) [('c', gap), (path, score), ..., (path, score)]  
 #       # len(col) = len(prot1) + 1  
 #    ]  
 #   
 def create_matrix_DP(pssm1, pssm2, psipred1, psipred2):  
   global aa_order  
   global gapcost  
     
   matrix_DP = [] # Returning matrix  
   scores = []  
   score = 0  
   score1 = 0 # score1 and score2 serve several purposes  
   score2 = 0  
   aa1_count = 0 # 2 indexes  
   aa2_count = 0  
   path = '' # 'd': diagonal, 'c': column, 'r': row  
     
   try:  
       
     # Fisrt row is set up with gap values  
     init_row = [('r', i * gapcost) for i in range(0, len(psipred2) + 1)]  
     matrix_DP.append(init_row)  
       
     for aa1 in pssm1:  
       aa1_count+=1  
         
       index1 = aa_order.index(aa1['aa']) # aa1 position in PSSM columns  
         
       scores = []  
       aa2_count = 0  
         
       for aa2 in pssm2:  
         aa2_count+=1  
           
         # First column is set up with gap values  
         if aa2_count == 1:  
           scores.append(('c', aa1_count * gapcost))  
           
         index2 = aa_order.index(aa2['aa'])  
           
         score1 = int(aa2['weights'][index1])  
         score2 = int(aa1['weights'][index2])  
           
         score = (score1 + score2) / 2.0  
         score += get_secondary_score(psipred1[aa1_count-1], psipred2[aa2_count-1])  
             
         # Now does apply the function for matrix positions:  
         #    Sij = max (score(-1,-1)+score, score(-1,0)+gap, score(0,-1)+gap)  
         # [1] is the score in the tuple (path, score)  
         score += matrix_DP[aa1_count - 1][aa2_count - 1][1]   
         score_1 = matrix_DP[aa1_count - 1][aa2_count][1] + gapcost  
         score_2 = scores[aa2_count - 1][1] + gapcost  
             
         if score_1 >= score and score_1 >= score_2:  
           path = 'c' # Column  
           score = score_1  
               
         elif score_2 >= score:  
           path = 'r' # Row  
           score = score_2  
               
         else:  
           path = 'd' # Diagonal  
             
         # Adds the score to the list of current row  
         scores.append((path, score))  
             
       # Adds a row of scores  
       matrix_DP.append(scores)  
         
   except ValueError, value:  
     raise Exception("Error. Maybe there is a rare aminoacid in the PSSM.\n"+value)  
   except IndexError, value:  
     raise Exception("Error. "+str(value))  
   
   return matrix_DP  
   
 # Recursive function that recovers the alignment  
 def visit_alignment(i, j, output):  
   global matrix_DP # This could be changed to be a parameter by reference  
   global prot1_pssm, prot2_pssm, prot1_psipred, prot2_psipred  
     
   if i == 0 and j == 0: return # END OF RECURSIVE SEARCH  
     
   join = ' '  
   if prot1_pssm[i - 1]['aa'] == prot2_pssm[j - 1]['aa']: join = '|'  
     
   path = matrix_DP[i][j][0]  
   score = matrix_DP[i][j][1]  
     
   if path == 'd':  
     visit_alignment(i - 1, j - 1, output)  
     output.append((prot1_psipred[i - 1], prot1_pssm[i - 1]['aa'],\  
      join, prot2_pssm[j - 1]['aa'], prot2_psipred[j - 1]))  
       
   elif path == 'c':  
     visit_alignment(i - 1, j, output)  
     output.append((prot1_psipred[i - 1], prot1_pssm[i - 1]['aa'], join, '-', '-'))  
       
   elif path == 'r':  
     visit_alignment(i, j - 1, output)  
     output.append(('-', '-', join, prot2_pssm[j - 1]['aa'], prot2_psipred[j - 1]))  
       
   else:  
     raise Exception("Invalid path "+path)  
     
   return  
   
   
 ############################ MAIN ##################################  
   
 # Global variables: program parameters  
 aa_order = ('ARNDCQEGHILKMFPSTWYV'); # This must copy aa order in PSSM columns  
 gapcost = -8  
 secondary_score = 2  
   
 # File locations: please, change this variables to suit your needs  
 prot1_pssm_file = '../pssm/1ngk_A.pssm'  
 prot2_pssm_file = '../pssm/1s69_A.pssm'  
 prot1_psipred_file = '../psipred/1ngk_A.psipred'  
 prot2_psipred_file = '../psipred/1s69_A.psipred'  
   
 # Printing title and author  
 print  
 print "*** "+__title__+" ***"  
 print "Author: "+__author__+" at "+__institution__  
   
 # Printing parameters  
 print  
 print "Aminoacids and order (PSSM files): "+aa_order  
 print "Gap penalty = "+str(gapcost)  
 print "Equal secondary structure score = "+str(secondary_score)  
   
 ### 1) Opening and parsing PSSM and PSIPRED files  
   
 prot1_pssm = []  
 prot2_pssm = []  
 prot1_psipred = []  
 prot2_psipred = []  
   
 try:  
   print  
   prot1_pssm = read_PSSM(prot1_pssm_file)  
   prot2_pssm = read_PSSM(prot2_pssm_file)  
   print "> PSSM files successfully opened."  
     
   prot1_psipred = read_PSIPRED(prot1_psipred_file)  
   prot2_psipred = read_PSIPRED(prot2_psipred_file)  
   print "> PSIPRED files successfully opened."  
     
 except IOError, value:  
   print "There has been an error reading input files:"  
   print value  
   
 ### 2) Creating the matrix of DP  
   
 matrix_DP = [[]]  
 len_prot1 = len(prot1_psipred)  
 len_prot2 = len(prot2_psipred)  
   
 try:  
   matrix_DP = create_matrix_DP(prot1_pssm, prot2_pssm, prot1_psipred, prot2_psipred)  
     
 except Exception, value:  
   print value  
   
 ### 3) Printing the maximum score  
   
 print  
 print "Maximum score "+str(matrix_DP[len_prot1][len_prot2][1])  
   
 ### 4) Recovering and printing the alignment  
   
 output = [] # [(ss1, aa1, align, aa2, ss2)]  
 # Lines to print  
 ss1 = ''  
 aa1 = ''  
 align = ''  
 aa2 = ''  
 ss2 = ''  
   
 try:  
   # Recursively recovering the alignment  
   visit_alignment(len_prot1, len_prot2, output)  
     
   # Printing the alignment  
   print  
   for nts in output:  
     ss1 += str(nts[0])  
     aa1 += str(nts[1])  
     align += str(nts[2])  
     aa2 += str(nts[3])  
     ss2 += str(nts[4])  
   print ss1  
   print aa1  
   print align  
   print aa2  
   print ss2  
     
 except Exception, value:  
   print "Error while tracing the alignment."  
   print value  
 else: # END OF PROGRAM  
   print  
   print "*** Alignment finished successfully ***"  

Esta sería la salida del programa, donde es importante fijarse en el alineamiento obtenido (que se muestra aquí en 2 fragmentos por motivo de espacio) y cómo, a pesar de que pocos aminoácidos coinciden, los estados de estructura secundaria se muestran alineados, lo que podría interpretarse como que la estructura está más conservada que la secuencia.

 *** Secondary structure and pairwise DP profile-profile alignment ***  
 Author: Carlos P Cantalapiedra at CSIC AulaDei Laboratory of Computational Biology  
   
 Aminoacids and order (PSSM files): ARNDCQEGHILKMFPSTWYV  
 Gap penalty = -8  
 Equal secondary structure score = 2  
   
 > PSSM files successfully opened.  
 > PSIPRED files successfully opened.  
   
 Maximum score 444.0  
   
 CCHHHHHCCHHHHHHHHHHHHHHHHCCHHHHHHCCCCCHHHHHHHHHHHHHHHHCCCCCCCCCCCCCCH  
 KSFYDAVGGAKTFDAIVSRFYAQVAEDEVLRRVYPEDDLAGAEERLRMFLEQYWGGPRTYSEQRGHPRL  
   |  ||  | | || | |     | |    ||  ||  |       
 STLYEKLGGTTAVDLAVDKFYERVLQDDRIKHFFADVDMAKQRAHQKAFLTYAFGGTDKYDGRYMREAH  
 CCHHHHHCCHHHHHHHHHHHHHHHHCCHHHHHHCCCCCHHHHHHHHHHHHHHHHCCCCCCCCCCHHHHH  
   
 HHCCCCCCCCHHHHHHHHHHHHHHHHHHCCCCCCHHHHHHHHHHHHHHHH--HHHCCCC  
 RMRHAPFRISLIERDAWLRCMHTAVASIDSETLDDEHRRELLDYLEMAAH--SLVNSPF  
        ||            |    ||   |    
 KELVENHGLNGEHFDAVAEDLLATLKEMG---VPEDLIAEVAAVAGAPAHKRDVLNQ--  
 HCCCCCCCCCHHHHHHHHHHHHHHHHHCC---CCHHHHHHHHHHHHHHHHHHHHHCC--  
   
 *** Alignment finished successfully ***