clear all;
close all;
clc;
%*************************************************************************
%*  Peter Mičík
%*  Program Shannon entropy - výpočet entropie sekvencie
%*  vstup - DNA sekvencia (FASTA, EMBL, GENBANK)
%*  výstup - graf (figure)
%*************************************************************************
reply1 = input('Zvožte typ formátu DNA sekvencie (FASTA,EMBL,GENBANK):', 's');
switch reply1
    case 'FASTA'
        sequence = input('\nZadajte názov sekvencie s priponou fasta: ','s');
        [header seq_in] = fastaread (sequence);   % Načítanie DNA sekvencie fasta formatu
        N = length(seq_in); % Dĺka sekvencie
        alpha = 'ACGT'; % abeceda sekvencie
    case 'EMBL'
        sequence = input('\nZadajte názov sekvencie s priponou: ','s');
        seq_in=emblread (sequence, 'SequenceOnly',true); 
        N = length(seq_in); % načítanie z EMBL
        alpha = 'acgt';
    case 'GENBANK'
        sequence = input('\nZadajte názov sekvencie s priponou: ','s');
        seq=genbankread(sequence); % načítanie z genbank
        seq_in=seq.Sequence;
        N = length(seq_in);
        alpha = 'acgt';
end

reply2 = input('\nChcete zada úsek sekvencie?(A/N): ', 's');
if reply2=='A'
    str1=sprintf('\nZvožte interval sekvencie od 1 do %d. Od: ',N);
    a=input(str1);
    str2=sprintf('\nDo: ');
    b=input(str2);
    seq_in=seq_in(a:b);
    N=length(seq_in);
end

bases=basecount(seq_in); % Počet jednotlivých nukleotidov
[dimer, percent]=dimercount(seq_in); % Počet dvojic a ich pravdepodobnos

seqint = nt2int(seq_in(1:N));
N=length(seqint);
pBi = [bases.A/N bases.C/N bases.G/N bases.T/N]; % Určenie prvkov Shannonovej entropie
pBiSj = percent;
    
    for i = 1:1:N-1
        Bi = seqint(i);
        Sj = seqint(i+1);
        entropy(i+1) = entropy(i)+(pBiSj(Bi,Sj)*log(pBiSj(Bi,Sj)/pBi(Bi))); % Výpočet entropie
    end
entropy=-entropy;    
entropia=entropy(end);

reply2 = input('\nChcete vykresli Shannonovu entropiu?(A/N): ', 's');
if reply2=='A'
    figure(1)       % Vykreslenie Shannonovej entropie
    plot(1:N,entropy);
    title('Vývoj entropie v závislosti na dĺke sekvencie');
    xlabel('Dĺka sekvencie (počet báz)');
    ylabel('Shannon entropy')
    str=sprintf('\nShannonova entropia je: %d ',entropia);
    disp(str);
else
    str=sprintf('\nShannonova entropia je: %d ',entropia);
    disp(str);
end