% SKRIPT PRO VYHODNOCENÍ DETEKCE PŘEKRYVŮ NA UMĚLÝCH DATECH
%----------------------------------------------------------
% skript vychází z vytvořené funkce
% dle nastavených parametrů vyhledává nejlepší překryv mezi všemi čteními z obou konců
% nastavitelná je minimální a maximální požadovaná délka překryvů, práh detekce, numerická reprezentace a metrika vyhodnocení podobnosti
% k testování jsou použita čtení, která byla uměle vygenerováná, čtení se  překrývají po sobě
% po detekci jsou srovnáný pozice detekovaných a známých překryvů
% je vypočítána senzitivita detekce
% po uběhnutí skriptu jsou zobrazeny výsledky
%%
%NASTAVENI POCATECNICH PODMINEK
minm=10; %minimalni velikost prekryvu
maxm=50; %maximalni delka prekryvu
prah=0.7; % minimalni pozadovane skore prekyvu
%NASTAVENI NUMERICKÉ REPREZENTACE a METRIKY
NR=4;
metrika='pearson'; %% 1 -  pearson, 2 - spearman
%NAČTENÍ DAT
[hlavicka,sekvence]=fastaread('50-cteni-delka-100-chybovost 0.5 % prekryv 20.fasta');
pocet_cteni=length(sekvence);
%% PŘEVEDENÍ DO ZVOLENÉ NUMERICKÉ REPREZENTACE
[signal,r,nazev]=prevod_NR(sekvence,NR); % koeficient r udava rozmer signalu
prah=prah*r;
%% VYTVOŘENÍ DVOU PRACOVNÍCH SEKVENCÍ S PREKRYVEM, SROVNANI JEDNOTLIVYCH PREKRYVU
%inicializace
score_max1=0;
score_max2=0;
score_1=zeros(1,r);
score_2=zeros(1,r);
tabulka_score1=zeros(pocet_cteni);% tabulka, do které se bude ukládat nejlepší skóre překryvu pro každé čtení
tabulka_score2=zeros(pocet_cteni);% tabulka, do které se bude ukládat nejlepší skóre překryvu pro každé čtení
tabulka_delkaprekryvu1=zeros(pocet_cteni);% tabulka, do které se bude ukládat délka nejlépe vyhodnoceného překryvu pro každé čtení
tabulka_delkaprekryvu2=zeros(pocet_cteni);% tabulka, do které se bude ukládat délka nejlépe vyhodnoceného překryvu pro každé čtení
prekryv1=0;
prekryv2=0;
for k=1:pocet_cteni 
    s1=signal{1,k}; % 1. porovnávané čtení
    delka_cteni(k)=length(s1);
    for l=1:pocet_cteni
    s2=signal{1,l}; %2.porovnávané čtení
    if k==l || l<k %neprovádim srovnání dvou stejných sekvencí, nebo sekvencí, které se již srovnávaly
        continue
    end
    l1=length(s1);
    l2=length(s2);
      for m=minm:maxm % m označuje velikost okna (překryvu), pro které počítáme metriku podobnosti mezi čteními
       % Bere se okno o velikosti m ze začátku prvního porovnávaného čtení a m okno od konce druhého porovnávaného čtení
       % pro výpočet score 1, pro výpočet score 2 okno zabírá konec prvního čtení a začátek druhého čtení
       % velikost okna se pohybuje od minimalni po maximální nastavenou hodnotu
       % v případě trojrozměrného signálu je vypočítána metrika pro každou složku a dále se pracuje se sumou všech hodnot 
            switch metrika
                case {'spearman',2}
                    for o=1:r
        score_1(o)=1-pdist2(s1(o,1:m),s2(o,l2-m+1:l2),'spearman'); %beru okno o velikosti m ze začátku první sekvence a od konce druhé sekvence
        score_2(o)=1-pdist2(s1(o,l1-m+1:l1),s2(o,1:m),'spearman'); %beru okno o velikosti m od konce první sekvence a ze začátku druhé sekvence
                    end
                case {'pearson',1}
                    for o=1:r
        score_1pom=corrcoef(s1(o,1:m),s2(o,l2-m+1:l2)); %beru okno o velikosti m ze začátku první sekvence a od konce druhé sekvence
        score_1(o)=score_1pom(1,2);
        score_2pom=corrcoef(s1(o,l1-m+1:l1),s2(o,1:m)); %beru okno o velikosti m od konce první sekvence a ze začátku druhé sekvence
        score_2(o)=score_2pom(1,2);
                    end
            end 
        score_1=sum(score_1); % v případě vícerozměrného signálu se počítá se sumou metrik
        score_2=sum(score_2);
        if score_1>=score_max1 %srovnání skóre s dosud maximálním dosaženým
        score_max1=score_1;
        prekryv1=m; %velikost překryvu
        end
        if score_2>=score_max2
        score_max2=score_2;
        prekryv2=m;
        end                       
       end
       tabulka_score1(l,k)=score_max1; 
       tabulka_score2(l,k)=score_max2;
       tabulka_delkaprekryvu1(l,k)=prekryv1;
       tabulka_delkaprekryvu2(l,k)=prekryv2;
       score_max1=0;
       score_max2=0;
    end 
end
% Doplněni tabulek 
 tabulka_score1=tabulka_score1+tabulka_score2';
 tabulka_score2=tabulka_score1';
 tabulka_delkaprekryvu1=tabulka_delkaprekryvu1+tabulka_delkaprekryvu2';
 tabulka_delkaprekryvu2=tabulka_delkaprekryvu1';
%% VYHODNOCENI PREKRYVU 
hodnota_score1=zeros(1,pocet_cteni); %Proměnná k ukládáni nejlepšího skóre překryvu ze začátku čtení
hodnota_score2=zeros(1,pocet_cteni); %Proměnná k ukládání nejlepšího skóre překryvu konce čtení
pozice_prekryvu1=zeros(1,pocet_cteni); %Proměnná k ukládání pozice nejlepšího překryvu pro každou sekvenci
pozice_prekryvu2=zeros(1,pocet_cteni); %Proměnná k ukládání pozice nejlepšího překryvu pro každou sekvenci
pocet_prekryvu1=zeros(1,pocet_cteni); %počet nalezených prekryvů pro každou sekvenci, ukládá se ten nejdelší překryv s nejlepším skóre
pocet_prekryvu2=zeros(1,pocet_cteni); %počet nalezených překryvů pro každou sekvenci
tabulka_score1=round(tabulka_score1,3); % zaokrouhlení, pro dosažení stejného počtu desetinných míst
tabulka_score2=round(tabulka_score2,3);
delka_prekryvu1=zeros(1,pocet_cteni); 
delka_prekryvu2=zeros(1,pocet_cteni);
% Zjištění pozice maximálního překryvu, hledáme překryv s nejlepším skóre,v
% případě více shodných výsledků se detekuje ten nejdelší překryv
 for i=1:pocet_cteni
    maxval1= max(tabulka_score1(:,i));
    maxval2= max(tabulka_score2(:,i));
    hodnota_score1(i)=maxval1;
    hodnota_score2(i)=maxval2;
    idx1 = find(tabulka_score1(:,i)== maxval1);
    idx2 = find(tabulka_score2(:,i)== maxval2);
    pocet_prekryvu1(i)=length(idx1);
    pocet_prekryvu2(i)=length(idx2);
    if length(idx1)>1 && maxval1>=prah 
    [hod,poz]=max(tabulka_delkaprekryvu1(idx1,i));
    pozice_prekryvu1(i)=idx1(poz);
    delka_prekryvu1(i)=hod;
    elseif maxval1<prah
        pozice_prekryvu1(i)=0; %zadny prekryv nedetekujeme
        pocet_prekryvu1(i)=0;
    else
        pozice_prekryvu1(i)=idx1;
        delka_prekryvu1(i)=tabulka_delkaprekryvu1(idx1,i);
    end
    
    if length(idx2)>1 && maxval2>=prah
    [hod,poz]=max(tabulka_delkaprekryvu2(idx2,i));
    pozice_prekryvu2(i)=idx2(poz);
    delka_prekryvu2(i)=hod;
    elseif maxval2<prah
        pozice_prekryvu2(i)=0; %zadny prekryv nedetekujeme
        pocet_prekryvu2(i)=0;
    else
        pozice_prekryvu2(i)=idx2;
        delka_prekryvu2(i)=tabulka_delkaprekryvu2(idx2,i);
    end

 end
 %%
 %VYHODNOCENÍ VÝSLEDKŮ
 %Nastaveni srovnávacího vektoru - pozice překryvů (každé čtení se překrývá s tím nasledujícím zprava)
ind_prek=[2:pocet_cteni,0];
nula=pocet_cteni; % posledni sekvence se neprekryva s nicim
srovnej=(ind_prek==pozice_prekryvu2); %srovnani znamych prekryvu s pozicemi prekryvu nalezenych
srovnej(nula)=[]; %vynulovani na pozici, kde se nejedna o prekryv
TP=sum(srovnej); %pocet spravne nalezenych prekryvu
FN=sum(srovnej==0); %pocet spatne nalezenych prekryvu
senzitivita=TP/(TP+FN)*100; 
%%
%ZOBRAZENÍ VÝSLEDKŮ
disp ('VYHODNOCENÍ VÝSLEDKŮ')
disp(['Název použité numerické reprezentace: ' nazev])
disp(['Počet čtení ' num2str(pocet_cteni)])
disp(['Délka jednotlivých čteni= ' num2str(delka_cteni)])
disp (['Práh ' num2str(prah)])

disp ('Překryv čtení na začátku :')
disp(['Hodnota skóre ' num2str(hodnota_score1) ])
disp (['Pozice překrývajícího se čtení ' num2str(pozice_prekryvu1)])
disp(['Délka překryvu= ' num2str(delka_prekryvu1)])
disp(['Počet nalezených překryvů:' num2str(pocet_prekryvu1)])

disp ('Překryv čtení na konci :')
disp(['Hodnota skóre ' num2str(hodnota_score2) ])
disp (['Pozice překrývajícího se čtení ' num2str(pozice_prekryvu2)])
disp(['Délka překryvu= ' num2str(delka_prekryvu2)])
disp(['Pocet nalezenych prekryvu:' num2str(pocet_prekryvu2)])
disp(['Senzitivita= ' num2str(senzitivita) ' %'])
clearvars -except delka_prekryvu1 delka_prekryvu2 hodnota_score1 hodnota_score2 pozice_prekryvu1 pozice_prekryvu2 NR p1 p2 pocetprek1 pocetprek2 prah nazev sekvence signal tabulka_delkaprek1 tabulka_delkaprek2 tabulka_index1 tabulka_index2 tabulka_score1 tabulka_score2 
