% Identification of protein coding regions using the Modified
% Gabor-wavelet transform.
% by J. P. Mena-Chalco, H. Carrer, Y. Zana and R. M. Cesar-Jr.
%
% Copyright (C) 2005, 2006: J. P. Mena-Chalco
% jmena@vision.ime.usp.br.
% http://www.vision.ime.usp.br/~jmena/MGWT/
%
function varargout = MGWT(varargin)
% MGWT M-file for MGWT.fig
% MGWT, by itself, creates a new MGWT or raises the existing
% singleton*.
%
% H = MGWT returns the handle to a new MGWT or the handle to
% the existing singleton*.
%
% MGWT('Property','Value',...) creates a new MGWT using the
% given property value pairs. Unrecognized properties are passed via
% varargin to MGWT_OpeningFcn. This calling syntax produces a
% warning when there is an existing singleton*.
%
% MGWT('CALLBACK') and MGWT('CALLBACK',hObject,...) call the
% local function named CALLBACK in MGWT.M with the given input
% arguments.
%
% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES
% Edit the above text to modify the response to help MGWT
% Last Modified by GUIDE v2.5 15-Jun-2006 20:23:35
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @MGWT_OpeningFcn, ...
'gui_OutputFcn', @MGWT_OutputFcn, ...
'gui_LayoutFcn', [], ...
'gui_Callback', []);
if nargin & isstr(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
% --- Executes just before MGWT is made visible.
function MGWT_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin unrecognized PropertyName/PropertyValue pairs from the
% command line (see VARARGIN)
% Choose default command line output for MGWT
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes MGWT wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = MGWT_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
varargout{1} = handles.output;
% --- Executes on button press in loadSequence.
function loadSequence_Callback(hObject, eventdata, handles)
[fastaFile,dirName]=uigetfile('*.fasta','Select fasta file');
if fastaFile==0
warndlg('Input fasta file must be selected.',' Warning ')
end
set(handles.sequenceName,'String',char(strcat(dirName,fastaFile)));
% --- Executes during object creation, after setting all properties.
function thresholdSlider_CreateFcn(hObject, eventdata, handles)
usewhitebg = 1;
if usewhitebg
set(hObject,'BackgroundColor',[.9 .9 .9]);
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on slider movement.
function thresholdSlider_Callback(hObject, eventdata, handles)
t = get(handles.thresholdSlider,'value');
set(handles.thresholdValue,'String',num2str(t));
% --- Executes on button press in execute.
function execute_Callback(hObject, eventdata, handles)
% a1 = 0.20; % firstScale
% a2 = 0.70; % lastScale
% ns = 10; % numberOfScales
% tv = 0.85; % thresholdValue
a1 = str2double( get(handles.firstScale,'String') );
a2 = str2double( get(handles.lastScale,'String') );
ns = str2double( get(handles.numberOfScales,'String') );
tv = str2double( get(handles.thresholdValue,'String') );
tamanho = 1200;
escalas = exp(linspace(log(a1),log(a2),ns));
fastaFile = char( get(handles.sequenceName,'String') );
if isempty(fastaFile)
warndlg('DNA sequence must be selected.',' Warning ')
else
fastaArray = textread(fastaFile,'%s');
i=1; while (strmatch('>',fastaArray(i))) i=i+1; end; % ">"
sequence = buildDNAStrand(fastaArray(i:end));
sequenceLen = length(sequence);
[uA,uC,uG,uT] = DNA2sequences(sequence);
twA = zeros(ns, sequenceLen);
twC = zeros(ns, sequenceLen);
twG = zeros(ns, sequenceLen);
twT = zeros(ns, sequenceLen);
for index = 1:ns
a = escalas(index);
mgw = mgabor(tamanho,a);
twA(index,:) = abs( convn(uA,mgw,'same') ).^2;
twC(index,:) = abs( convn(uC,mgw,'same') ).^2;
twG(index,:) = abs( convn(uG,mgw,'same') ).^2;
twT(index,:) = abs( convn(uT,mgw,'same') ).^2;
end;
tw = twA + twC + twG + twT;
CDS = sum(tw);
axes(handles.spectrogramA);
imagesc(twA);
set(handles.spectrogramA,'XTickLabel','');
set(handles.spectrogramA,'YTickLabel','');
axes(handles.spectrogramC);
imagesc(twC);
set(handles.spectrogramC,'XTickLabel','');
set(handles.spectrogramC,'YTickLabel','');
axes(handles.spectrogramG);
imagesc(twG);
set(handles.spectrogramG,'XTickLabel','');
set(handles.spectrogramG,'YTickLabel','');
axes(handles.spectrogramT);
imagesc(twT);
set(handles.spectrogramT,'XTickLabel','');
set(handles.spectrogramT,'YTickLabel','');
axes(handles.spectrogram);
imagesc(tw);
set(handles.spectrogram,'YTick',[1 ns]);
set(handles.spectrogram,'YTickLabel',str2mat(num2str([a1 a2]')));
axes(handles.projectionA);
plot(sum(twA));
set(handles.projectionA,'XTickLabel','');
set(handles.projectionA,'YAxisLocation','right');
axes(handles.projectionC);
plot(sum(twC));
set(handles.projectionC,'XTickLabel','');
set(handles.projectionC,'YAxisLocation','right');
axes(handles.projectionG);
plot(sum(twG));
set(handles.projectionG,'XTickLabel','');
set(handles.projectionG,'YAxisLocation','right');
axes(handles.projectionT);
plot(sum(twT));
set(handles.projectionT,'XTickLabel','');
set(handles.projectionT,'YAxisLocation','right');
axes(handles.projection);
plot(CDS);
set(handles.projection,'YAxisLocation','right');
% final
CDSlist = seq2CDS(DNAthreshold(CDS, tv));
axes(handles.final);
hold on;
for i = 1:2:length(CDSlist)
fill ([CDSlist(i) CDSlist(i) CDSlist(i+1) CDSlist(i+1)],[0,1,1,0],[0.784 0.784 0.784],'LineStyle','none');
end;
plot(CDS/max(CDS));
set(handles.final,'YAxisLocation','right');
box on;
% list of Protein coding regions (CDS)
CDSstring = '';
for i = 1:2:length(CDSlist)
CDSstring= strvcat( CDSstring, strcat(num2str(floor((i+1)/2)),' : ', num2str(CDSlist(i)),' - ', num2str(CDSlist(i+1)),' (', num2str(CDSlist(i+1)-CDSlist(i)+1),'bp)' ) );
end;
set(handles.codingRegions,'String',CDSstring);
end
% --- Executes on button press in close.
function close_Callback(hObject, eventdata, handles)
close all;
% --- Executes during object creation, after setting all properties.
function firstScale_CreateFcn(hObject, eventdata, handles)
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
function firstScale_Callback(hObject, eventdata, handles)
% --- Executes during object creation, after setting all properties.
function lastScale_CreateFcn(hObject, eventdata, handles)
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
function lastScale_Callback(hObject, eventdata, handles)
% --- Executes during object creation, after setting all properties.
function numberOfScales_CreateFcn(hObject, eventdata, handles)
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
function numberOfScales_Callback(hObject, eventdata, handles)
%% ----------------------------------------------------------------------- %%
%% ----------------------------------------------------------------------- %%
%% ----------------------------------------------------------------------- %%
%% Modification of Gabor wavelet fuction
%% ----------------------------------------------------------------------- %%
function out = mgabor(limit,a)
w0 = 1*limit/3; % frequency
w = linspace(0,limit,limit);
W = sqrt(2*pi)*( exp( -(a^2)*((w-w0).^2)/2) ); % gaussian function centered in N/3
out= ifftshift( ifft(W) );
%% ----------------------------------------------------------------------- %%
%% DNA2sequences: Make the four binary sequences %%
%% ----------------------------------------------------------------------- %%
function [uA, uC, uG, uT] = DNA2sequences(sequence)
N = length(sequence);
% Create initialized indicators for each nucleotide
uA(N) = 0;
uC(N) = 0;
uG(N) = 0;
uT(N) = 0;
% Assign the value 1 to the corresponding indexes
uA(find(upper(sequence)=='A')) = 1;
uC(find(upper(sequence)=='C')) = 1;
uG(find(upper(sequence)=='G')) = 1;
uT(find(upper(sequence)=='T')) = 1;
%% ----------------------------------------------------------------------- %%
%% DNAthreshold %%
%% ----------------------------------------------------------------------- %%
function matrix = DNAthreshold(matrix, percent)
if (percent==0)
return;
end
[m,n] = size(matrix);
for i = 1:m
seqL = sort(matrix(i,:));
bias = seqL(floor(percent*n)+1);
for j=1:n
if (matrix(i,j)<=bias)
matrix(i,j)=0;
end;
end;
end;
%% ----------------------------------------------------------------------- %%
%% seq2CDS %%
%% ----------------------------------------------------------------------- %%
function CDS = seq2CDS(sequence)
CDS = [];
i=1;
while i<=length(sequence)
if (sequence(i)~= 0 )
CDS = [CDS , i];
while (i<=length(sequence) & sequence(i)~=0)
i= i+1;
end;
CDS = [CDS , i-1];
end;
i = i+1;
end;
%% ----------------------------------------------------------------------- %%
%% buildDNAStrand %%
%% ----------------------------------------------------------------------- %%
function sequence = buildDNAStrand(DNAMatrix)
DNAMatrix = char(DNAMatrix);
[nRows, nBases] = size(DNAMatrix);
for j = 0:nRows*nBases-1
if (DNAMatrix(floor(j/nBases)+1, mod(j,nBases)+1) ~= ' ')
sequence(j+1) = DNAMatrix(floor(j/nBases)+1, mod(j,nBases)+1);
else
break;
end;
end;
%end;
