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