Alineamiento con transformadas de Fourier

En 1992 el grupo de Ilya Vakser publicó en PNAS el método fundamental para poder hacer simulaciones de docking de biomoléculas de manera eficiente, disminuyendo considerablemente el tiempo de cálculo aplicando transformadas rápidas de Fourier. Diez años más tarde se publicó la primera versión del exitoso programa de alineamiento múltiple de secuencias MAFFT, del que ya hemos hablado en otras entradas, que aplica ideas similares para el problema de encajar secuencias similares, como se hace al alinear secuencias.

Figura tomada de http://www.ncbi.nlm.nih.gov/pubmed/12136088. donde se muestran dos picos del análisis de Fourier que se corresponden con dos subalineamientos locales.

En esta entrada mi intención es explicar el fundamento de esta técnica recurriendo al lenguaje Octave, la versión open source de MatLab. El código es el siguiente:

 % alineamiento (sin gaps) de secuencias proteicas usando la FFT  
 % escrito en Octave, que es compatible con Matlab   
 % requiere la libreria 'bioinfo'   
   
 clear all  
   
 %% Conversor binario de secuencias   
 %% Parametros:   
 %% 1) longitud de la secuencia S  
 %% 2) longitud del fragmento F (de S)  
 function [S , F] = aa2bits( sequence , fragment )  
   
 % cada residuo se representa por un numero del 1 al 20  
 % http://www.mathworks.com/help/toolbox/bioinfo/ref/aa2int.html  
 [ seq ]= aa2int( sequence )  
 [ frag ] = aa2int( fragment )   
   
 % codificamos residuos como cadenas binarias de ancho fijo (20 columnas)  
 S=[];  
 for i=1:length(seq)  
   % cada residuo ocp  
   S=[S ones(1,seq(i)) zeros(1,20-seq(i)) ];  
 end  
   
 F=[];  
 for i=1:length(frag)  
   F=[F ones(1,frag(i)) zeros(1,20-frag(i)) ];  
 end  
   
 % completa F hasta igualar la longitud de S  
 F=[F zeros(1,length(S)-length(F))];  
   
 endfunction % no se usa en Matlab  
   
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
   
 % define la secuencia S y el fragmento F que queremos alinear  
 sequence = 'SDEVRKNLMDMFRDRQAFSEHTWKMLLSVCRSWAAWCKLNNRKWFPAEPEDVRDYLLY';  
 fragment = 'KMLNSVCRSWWWWC';  
   
 % convertimos ambas secuencias en dos segnales binarias,  
 % de manera que cada residuo se representa por un natural del 1 al 20  
 [S , F] = aa2bits(sequence,fragment);  
   
 % calculamos el alineamiento de F con S optimizando la funcion objetivo Fobj,  
 % que en este ejemplo es simplemente el producto (como un AND binario)  
 % obtenido al desplazar F de izq a der sobre S:   
 % producto(n) = sum { F(i) * S(i-n) }  
 %        i   
 % aprovechamos que Fobj se puede expresar por medio de transformadas  
 % de Fourier para reducir las operaciones necesarias (de N^2 a 4NlogN)  
 % https://sites.google.com/site/cartografiaygeodesia/prac7.pdf  
 % http://www.pnas.org/content/89/6/2195.abstract  
   
 FTsec = fft(S);  
 FTfrag = fft(F);  
   
 FTproducto = conj(FTfrag) .* FTsec; % producto termino a termino  
   
 % Deshacemos transformacion y localizamos valor maximo,   
 % la posicion que optimiza el alineamiento   
 producto = ifft(FTproducto);   
 [valor maxpos] = max(producto);  
   
 % imprime alineamiento sobre secuencicas binarias   
 tit = sprintf("Optimal alignment, optimal position = %d",maxpos);  
 plot(S*0.75)  
 title(tit)  
 xlabel('sequence position')  
 hold on  
 plot([zeros(1,maxpos) F*0.9] , 'r');  
 axis([0 length(S) 0 1]);  
 legend('sequence','fragment');  
 print -dpng figure_align.png;  
   
 % alineamiento sobre secuencias originales, cambiando maxpos de escala  
 align = printf("Alignment:\nS %s\nF %s%s\n",sequence,blanks((maxpos-1)/20),fragment);  

Como resultado obtendremos un alineamiento como éste, que :

S SDEVRKNLMDMFRDRQAFSEHTWKMLLSVCRSWAAWCKLNNRKWFPAEPEDVRDYLLY
F                        KMLNSVCRSWWWWC
 

Que se corresponde con éste producto de Fourier:

En Perl podríamos usar el módulo Math::FFT para este fin, o la GNU Scientific Library que ya revisamos en otra entrada. Un saludo,
Bruno