
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% EE 675 - Fall 2015
%%%% Simple Code to show How to Open Genomic Data Files
%%%%
%%%% Your project should try to Find Differences between EBOLA and MERS
%%%% Possible tool: Look at their Spectra in the Frequency Domain
%%%%
%%%% Data For All Genomes is Given in Standard FASTA format
%%%%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

range= 1:1000;

%% Note: for EBOLA virus, many samples are given - Just choose One to Use
x= fastaread('AB050936v1.fa');
x.Sequence          % show the sequence in its original form
ebola= double(x.Sequence);  % convert to genomic signal
Ebola= x;

y= fastaread('mers.fa') 
y.Sequence          % show the sequence in its original form
mers = double(y.Sequence);  % convert to genomic signal
Mers= y;

%%% ASCII codes for ACGT are 65, 67, 71, 84
%%%%% Convert all numbers into 0,1,2,3 (2-bits) 
m= find(ebola == 65);
ebola(m)= 0;
m= find(ebola == 67);
ebola(m)= 1;
m= find(ebola == 71);
ebola(m)= 2;
m= find(ebola == 84);
ebola(m)= 3;

m= find(mers == 65);
mers(m)= 0;
m= find(mers == 67);
mers(m)= 1;
m= find(mers == 71);
mers(m)= 2;
m= find(mers == 84);
mers(m)= 3;


figure(1)
subplot(2,1,1), plot(ebola(range)), title('EBOLA')
subplot(2,1,2),plot(mers(range)), title('MERS')


N= 1024*32;
febola= abs(fftshift(fft(ebola- mean(ebola),N)));
fmers= abs(fftshift(fft(mers- mean(mers),N)));

[mx1, lx1] = max(febola)
[mx2, lx2] = max(fmers)

figure(2)
subplot(2,1,1), plot(febola), title('EBOLA Spectrum')
subplot(2,1,2),plot(fmers), title('MERS Spectrum')


size1= size(ebola);
Nebola= size1(2);
size1= size(mers);
Nmers= size1(2);
 
%%%% spectrogram
M= 1024;
M2= M/8;
spectrogram_ebola=  zeros(ceil((Nebola-M+1)/M2),M);
spectrogram_mers=   zeros(ceil((Nmers -M+1)/M2),M);
ind1= 1;
for k1= M:M2:Nebola
    sig= ebola(k1-M+1: k1);
    sig= sig- mean(sig);
    spectrogram_ebola(ind1,:)= abs(fftshift(fft(sig))).^2;
    ind1= ind1+1;
end
ind1= 1;
for k1= M:M2:Nmers
    sig= mers(k1-M+1: k1);
    sig= sig- mean(sig);
    spectrogram_mers(ind1,:)= abs(fftshift(fft(sig))).^2;
    ind1= ind1+1;
end




figure(3)
subplot(2,1,1), imagesc((spectrogram_ebola')), title('EBOLA Spectrogram')
subplot(2,1,2), imagesc((spectrogram_mers')), title('MERS Spectrogram')
colormap(gray)    

figure(4)
subplot(2,1,1), imagesc(log(spectrogram_ebola'+1)), title('EBOLA Spectrogram')
subplot(2,1,2), imagesc(log(spectrogram_mers'+1)), title('MERS Spectrogram')
colormap(gray)    
    
figure(5)
subplot(2,2,1), ntdensity(Ebola)

figure(6)
subplot(2,2,2), ntdensity(Mers)

bases_ebola = basecount(Ebola)
bases_mers = basecount(Mers)

figure(7)
subplot(1,2,1)
basecount(Ebola,'chart','pie');
title('Distribution of Bases for Ebola Genome');
subplot(1,2,2)
basecount(Mers,'chart','pie');
title('Distribution of Bases for Mers Genome');

figure(8)
dimers = dimercount(Ebola,'chart','bar')
title('Ebola Genome Dimer Histogram');
figure(9)
dimers = dimercount(Mers,'chart','bar')
title('Mers Genome Dimer Histogram');

figure(10)
codoncount(Ebola, 'FIGURE', true, 'GENETICCODE',1)
figure(11)
codoncount(Mers, 'FIGURE', true, 'GENETICCODE',1)
%codoncount_dif= codoncount_ebola- codoncount_mers

