Transformada de Burrows-Wheeler (para comprimir secuencias y buscar patrones)

Hola,
como posiblemente haya vacaciones en el blog estas próximas semanas, os dejo una entrada algo más compleja, donde se muestran los fundamentos de la transformación de Burrows-Wheeler aplicada a la búsqueda de patrones en secuencias biológicas, como continuación natural de los vectores de sufijos y los árboles binarios, que vimos en una entrada anterior.

(fuente: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.12.1158)

Vaya por delante que las implementaciones reales de este algoritmo (como Bowtie o BWA) son bastante más complejas (porque permiten búsquedas inexactas) y eficientes (porque comprimen el corpus para que se pueda cargar en la memoria RAM) que el código de esta entrada, pero su filosofía es similar:

1) se precalcula la transformada de un corpus de secuencias (y se comprime)
2) se buscan patrones exactos/inexactos sobre el corpus precalculado (comprimido)

Para aprender más, además de la lectura del artículo original de Burrows y Wheeler, un material recomendable para iniciarse en estos temas puede ser el libro de Gabriel Valiente  'Combinatorial pattern matching algorithms in computational biology using Perl and R'.

El siguiente código en Perl muestra un prototipo con estas ideas, empleando para la búsqueda el árbol binario implícito en el vector de sufijos:

 use strict;  
   
 my $VERBOSE = 1;  
 my @ALPHABET = ( '$', 'a' .. 'z' ); # orden lexicografico, $ marca el final   
   
 my ($text,$pattern) = ('mississippi','ssi');           
            
 printf("# original text: '%s'\n\n",$text);  
   
 my ($bwt,$start,$F) = BWT($text,$VERBOSE);   
   
 printf("# BWT('%s') = ('%s',%d) | F: %s\n",$text,$bwt,$start,$F);   
   
 # ahora se pueden comprimir $bwt y $F, por ejemplo, con run-length coder:  
 # 'mmmmmsssiiii' (12 bytes) quedaría como 'm5s3i4' (6 bytes)  
   
 my ($orig_text,$ref_array_W) = revBWT($bwt,$start,$VERBOSE);  
   
 printf("# BWTinv('%s',%d) = '%s'\n\n",$bwt,$start,$orig_text);   
   
 print "# searching for '$pattern' matches in '$text' [base 0]:\n";  
 match_pattern_BWT($pattern,$F,$start,$ref_array_W);  
   
   
   
 # Calcula un vector de sufijos de manera que las diferentes subcadenas   
 # de una secuencia queden preordenadas lexicograficamente  
 sub make_suffix_array  
 {  
   my ($seq,$verbose) = @_;       
   
   my $l = length($seq);  
   my @suffix = (0 .. $l-1); # en base 0  
   
   # ordena los length($seq) posibles sufijos lexicográficamente (cmp),  
   # este paso consume RAM y es uno de los pasos que se puede optimizar  
   @suffix = sort {substr($seq,$a) cmp substr($seq,$b)} (@suffix);  
   if($verbose)  
   {  
     print "# suffix array for $seq :\n";  
     foreach my $suf (@suffix)  
     {   
       printf("%3d %s%s\n",  
         $suf,substr($seq,$suf),substr($seq,0,$l-($l-$suf)));   
     }  
     print "\n";  
   }  
   return @suffix;  
 }  
   
 # Devuelve: 1) string BWT, 2) índice $start, que indica la fila del vector  
 # de sufijos que contiene la cadena original, y 3) string $F, la columna   
 # F1..Fn concatenada.  
 # BWT realmente es la columna de caracteres que preceden a la primera (F)  
 # de un vector de sufijos, o la última (L) si en vez de sufijos la cadena  
 # se permuta circularmente  
 # F1 .............$ L1  
 # F2 ............$. L2  
 # Fn ......$....... Ln  
 sub BWT  
 {  
   my ($seq,$verbose) = @_;  
       
     $seq = $seq . '$'; # marca el final   
       
   my @suffix_rotation = make_suffix_array($seq,$verbose);  
   my ($bwt,$start,$F) = ('',0);  
   foreach my $c (0 .. length($seq)-1)  
   {  
     $bwt .= substr($seq,$suffix_rotation[$c]-1,1);  
         $F .= substr($seq,$suffix_rotation[$c],1);  
     if($suffix_rotation[$c] == 0)  
     {  
       # fila del vector de sufijos que contiene cadena original    
       $start = $c;   
     }  
   }  
   return ($bwt,$start,$F);  
 }  
   
 # Devuelve: 1) la cadena original correspondiente a una cadena  
 # previamente transformada con BWT y 2) una referencia a un array @W  
 # que contiene, para cada posicion del vector L (cadena transformada), la que  
 # le sigue en la secuencia original; ojo, el orden de @W se conserva en @F  
 sub revBWT  
 {  
   my ($bwt,$start,$verbose) = @_;  
   my (@L,@C,@K,@M,@W,@T,$original_string);  
   my ($c,$l,$r,$total,$last_char);  
     my $last_alphabet = scalar(@ALPHABET)-1;  
   
   # convierte cadena $bwt en array @L  
   @L = split(//,$bwt);  
     $last_char = scalar(@L)-1;  
   
   # calcula vectores auxiliares C y K, contienen las observaciones   
   # de cada letra del alfabeto en $bwt (@L)  
   foreach $l (0 .. $last_alphabet){ $K[$l] = 0 }  
       
   foreach $c (0 .. $last_char)  
   {   
     foreach $l (0 .. $last_alphabet)  
     {  
       if($L[$c] eq $ALPHABET[$l])  
       {  
         $C[$c] = $K[$l];    
         $K[$l]++;  
       }   
     }   
   }  
       
     foreach $l (0 .. $last_alphabet)  
   {  
     if($verbose && $bwt =~ /$ALPHABET[$l]/)  
     {  
       print "# $ALPHABET[$l] K=$K[$l]\n";  
     }               
   }  
   
   # calcula vector auxiliar M, que contiene la primera aparición  
   # de cada letra del alfabeto en vector implícito F, que es   
   # la primera columna del suffix array subyacente  
   $total = 0;  
   foreach $l (0 .. $last_alphabet)  
   {  
     $M[$l] = $total;  
     $total += $K[$l];    
     if($verbose && $bwt =~ /$ALPHABET[$l]/)      
     {  
       print "# $ALPHABET[$l] M=$M[$l]\n" if($verbose);  
     }  
   }  
   
   # calcula vector decodificador @W (forward), destruye @M  
   foreach $c (0 .. $last_char)  
   {  
     foreach $l (0 .. $last_alphabet)  
     {  
       if($L[$c] eq $ALPHABET[$l])  
       {  
         $W[$M[$l]] = $c;     
         $M[$l]++;   
       }    
     }    
   }  
   
   # decodifica cadena de texto original (@T)  
   $original_string = '';  
   $l = $start;  
   foreach $c (0 .. $last_char)  
   {  
     $original_string .= $L[$l] if($L[$l] ne '$');  
     $l = $W[$l]  
   }  
   
   return ($original_string,\@W);  
 }  
   
 # Funcion equivalente al operador cmp para comparar una cadena $pattern   
 # con la BWT de un texto, dada una fila $row del suffix_array.  
 # Devuelve < 0 si $pattern < $bwt (orden lexicográfico), > 0 en caso contrario,  
 # y 0 si la fila $r comienza con $pattern   
 sub BWTstrcmp  
 {  
   my ($pattern,$F,$row,$ref_W) = @_;  
     
   my @F = split(//,$F);  
   my $m = length($pattern);   
   my ($c,$match) = (0,'');    
       
   while($m >= 0 and $F[$row] eq substr($pattern,$c,1))  
   {  
         $match .= $F[$row];  
     $row = $ref_W->[$row];  
     $m--;  
     $c++;  
   }    
   
   if($m == 0){return 0 } # match  
   else{ return substr($pattern,$c,1) cmp $F[$row] }  
 }  
   
 # Procedimiento que busca ocurrencias de un patron recorriendo el vector de sufijos  
 # implicito en la BWT, por medio del vector @F. Imprime las posiciones encontradas  
 sub match_pattern_BWT  
 {  
   my ($pattern,$F,$start,$ref_W) = @_;  
     
     my ($med,$submed,$low,$high,$comp);  
     my $last = length($F);  
   $low = 0;  
   $high = $last;  
   
   while ($low <= $high)   
   {  
     $med = int (($low+$high)/2);     
     $comp = BWTstrcmp($pattern,$F,$med,$ref_W);  
   
     if($comp < 0){ $high = $med-1 } # retrocedemos   
     elsif($comp > 0){ $low = $med+1 } # avanzamos  
     else   
     {  
       $submed = $med; # filas anteriores  
             while($submed > 0 &&   
         BWTstrcmp($pattern,$F,$submed,$ref_W) == 0)  
       {   
         printf("# match at position %d\n",  
           get_real_position($submed,$start,$ref_W));  
         $submed--;   
       }  
         
             $submed = $med + 1; # filas posteriores  
       while ($submed < $last &&   
         BWTstrcmp($pattern,$F,$submed,$ref_W) == 0)  
       {  
         printf("# match at position %d\n",  
           get_real_position($submed,$start,$ref_W));  
         $submed++;   
       }  
       last;  
     }  
   }  
   
 }  
   
 # Devuelve la posicion en el texto original de una fila de @F  
 sub get_real_position  
 {  
   my ($Fpos,$start,$ref_W) = @_;   
       
   my $real_pos = 0;  
   while($start != $Fpos)  
   {  
     $start = $ref_W->[$start];  
     $real_pos++;       
   }  
   return $real_pos;  
 }  
   

Extraer de un texto líneas con el contenido deseado

A continuación os dejo un sencillo código para extraer las líneas deseadas de un texto dándole como parámetros el texto y un array con las palabras o patrones que deben estar presentes en las líneas a extraer.

Modificaciones de este código nos pueden ser muy útiles en diversos problemas bioinformáticos, cuando tenemos archivos con muchas líneas de datos y sólo queremos utilizar o visualizar unas pocas.

En próximas entradas utilizaremos este código para extraer por ejemplo los átomos deseados de un fichero de coordenadas atómicas en formato PDB.

 my $source_text = "Lorem ipsum dolor sit amet, consectetur adipisicing elit, sed do eiusmod tempor incididunt ut labore et dolore magna aliqua.\n  
 Ut enim ad minim veniam, quis nostrud exercitation ullamco laboris nisi ut aliquip ex ea commodo consequat.\n  
 Duis aute irure dolor in reprehenderit in voluptate velit esse cillum dolore eu fugiat nulla pariatur.\n  
 Excepteur sint occaecat cupidatat non proident, sunt in culpa qui officia deserunt mollit anim id est laborum.";  
   
 my @patterns;  
 push(@patterns, 'adipisicing');  
 push(@patterns, '/dolor/');  
   
 my $extracted_text = join('',extract_lines_from_text($source_text,\@patterns));  
   
 print "$extracted_text\n";  
   
 # Extract lines from text with the desired patterns  
 sub extract_lines_from_text {  
   
      my ($text, $patterns) = @_;  
   
      my @data;  
      my @lines = split("\n",$text);  
   
      foreach my $line (@lines){  
           foreach my $pattern (@{$patterns}){  
                if ($pattern =~ /^\/(.+)\/$/){  
                     if ($line =~ /$1/){  
                          push(@data,$line);  
                          last;  
                     }  
                } else {  
                     if ($line =~ /\Q$pattern\E/){  
                          push(@data,$line);  
                          last;  
                     }  
                }  
           }  
      }  
   
      return @data;  
   
 }