I want to find DFT of a input sequence. I have run the following code but not getting desired output. Is there any mistake? please help.

if true
  % code
end
clear all;
close all;
clc
sample_dna = fastaread('sequence_AF099922.Fasta');
sample_dna = struct2cell(sample_dna);
sample_dna = cell2mat(sample_dna(2));
sample_dna = sample_dna(7021:15020);
cds = [928 1038 2527 2856 4113 4376 5464 5643 7254 7604];
for ctr = 1:2:length(cds)
actual_exons(cds(ctr):cds(ctr+1))=1;
end
actual_exons(cds(length(cds))+1:length(sample_dna)) =0;
[x] = dna_binary(sample_dna);
ln = length(x);
xk = zeros (1,ln);
% code block to find the DFT of the sequence
i = sqrt(-1);
for k = 0: ln-1
for n=0: ln-1
xk(k+1) = xk(k+1) + (x(n+1)*exp((-i)*2*pi*k*n/ln));
end
end
Y = abs(xk),^2;
P = Y/max(Y);
plot(actual_exons,'k:');
hold on
plot(P);
axis([0 8000 0 1.05])
xlabel('Nucleotide position');
ylabel('Output');
title('Detection of period three behaviour using DFT');

[Kodavati Mahendra]

Please explain more about what output you are getting, or the error you are getting :-)

[SUBHAJIT KAR]

Sir, I am getting a blank frame. The output is attached herewith.

[SUBHAJIT KAR]

The output.

[Kodavati Mahendra]

Can you please zoom in the lower frequencies ?

Answer 1

[Kodavati Mahendra]

line 26,

Y = abs(xk).^2;

you can use the "fft" function in MATLAB, Lines 18-25 can be replaced with

xk = fft(x,ln);

or even more simpler, Lines 17-25 can be replaced with

xk = fft(x);

Both are equivalent, rest seems fine. Also I do not have bioinformatics toolbox, so do not know about dna_binary function. Please add it to your code dna_binary

Regards,

Mahendra

[SUBHAJIT KAR]

dna_binary is a function which accepts DNA sequence comprising A,C,T,G and convert it to specific binary indicator sequence. Here is the coding.

if true
  % code
end
function[x] = dna_binary(sample_dna)
for i= 1:length(sample_dna)
if sample_dna(i)=='A' || sample_dna(i) =='a'
x(i) = 0.1260;
elseif sample_dna(i)=='T'||sample_dna(i) == 't'
x(i) = 0.1335;
elseif sample_dna(i) == 'C' || sample_dna(i) == 'c'
x(i) = 0.1340;
elseif sample_dna(i) = 'G' || sample_dna(i) == 'g'
x(i) = 0.0806;
end
end
end