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.
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
