24 de agosto de 2010

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)

19 de agosto de 2010

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)

6 de agosto de 2010

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';  
      }  
 }