clear all; close all; clc

%% načtení fasta souboru
% [header seq] = fastaread('X64775.fasta'); % oryza sativa s mikrosatelity
% [header seq] = fastaread('AL133493_human_chr21.fasta'); % human chr21 s minisatelity
% [header seq] = fastaread('randomseq_minisatelit.fasta'); % náhodná sekvence s minisatelity
% [header seq] = fastaread('neisseria gonorrhoeae clone pNG273.fasta'); % rozptýlené repetice bakterie n. gonorrhoeae 152, 25bp
% [header seq] = fastaread('Rabbit short interspersed C repeat.fasta'); % rozptýlené repetice typu SINE (13bp)
% [header seq] = fastaread('Caenorhabditis elegans Cosmid F56F11.fasta'); % gen F56F11 C. elegans s pěti exony
% [header seq] = fastaread('Mycoplasma pneumoniae.fasta'); % část genomu 1-100000bp
% [header seq] = fastaread('C. carsonella ruddii.fasta'); % celý genom 160000bp

seq = randseq(60000); % vygenerování náhodné sekvence o délce 60000 nukleotidů
   
% seq(15:15:end) = 'A'; % každý 15. nukleotid je A
% seq(13:13:end) = 'T'; % každý 13. nukleotid je T
% seq(9:9:end) = 'G';   % každý 9. nukleotid je G
% seq(11:11:end) = 'C'; % každý 11. nukleotid je C

seq(10000:3:10150) = 'C'; % trinukleotidová repetice CGG opakující se v délce 150 bp - mikrosatelit
seq(10001:3:10151) = 'G';
seq(10002:3:10152) = 'G';

% seq(10000:5:10150) = 'A'; % pětinukleotidová repetice AAATT opakující se v délce 150 bp - mikrosatelit
% seq(10001:5:10151) = 'A';
% seq(10002:5:10152) = 'A';
% seq(10003:5:10153) = 'T';
% seq(10004:5:10154) = 'T';

%% numerická reprezentace
L = length(seq); % délka sekvence

num_repr = zeros(4,L); 
num_repr(1,:)= seq == 'A'; % převod symbolického zápisu bází na numerickou formu
num_repr(2,:)= seq == 'T';
num_repr(3,:)= seq == 'C';
num_repr(4,:)= seq == 'G';

%% vytvoření spektrogramu
w = 1000; % window size - délka okna
o = 500;  % overlap - posun oken
pom = 1;  % pomocná pro vkládání segmentů do sloupců vedle sebe ve spektrogramu

% barevné mapovací vektory:
tr = 1; tg = 0; tb = 0; % Red   [1 0 0]
cr = 0; cg = 1; cb = 0; % Green [0 1 0]
ar = 0; ag = 0; ab = 1; % Blue  [0 0 1]
gr = 1/3; gg = 1/3; gb = 1/3; % Grey  [1/3 1/3 1/3]; Black [0 0 0]; Yellow [1 1 0]

for u = 1:o:L-w % procházení numerickou reprezentací 'num_repr' po sloupcích
    for v = 1:4 % procházení numerickou reprezentací 'num_repr' po řádcích = každý pro jiný nukleotid
        x =(num_repr(v,u:u+w)-mean(num_repr(v,u:u+w))); % procházíme vždy jeden řádek num. repr. ze všech čtyř pro všechny hodnoty sloupců 'u'
%         x = (hann(length(x1)))'.*x1; % vynásobení Hannovým oknem
        spektrum = DFT(x); % volání funkce 'DFT.m'
        spektrum_celkem(v,:)= spektrum.^2; % čtyři spektra pro A,T,C,G uložíme do matice o čtyřech řádcích
    end
    
    % mapování do RGB prostoru:
    R = ar.*(spektrum_celkem(1,:))+tr.*(spektrum_celkem(2,:))+cr.*(spektrum_celkem(3,:))+gr.*(spektrum_celkem(4,:));
    G = ag.*(spektrum_celkem(1,:))+tg.*(spektrum_celkem(2,:))+cg.*(spektrum_celkem(3,:))+gg.*(spektrum_celkem(4,:));
    B = ab.*(spektrum_celkem(1,:))+tb.*(spektrum_celkem(2,:))+cb.*(spektrum_celkem(3,:))+gb.*(spektrum_celkem(4,:));
    R=R'; % transponování vektorů
    G=G';
    B=B';
    
    % normalizace hodnot pixelů:
    MaxR = max(max(R)); MinR = min(min(R)); % normalizace hodnot RGB na rozsah 0-1
    nR = R - MinR;
    nnR = nR/(MaxR-MinR);
    
    MaxG = max(max(G)); MinG = min(min(G));
    nG = G - MinG;
    nnG = nG/(MaxG-MinG);
    
    MaxB = max(max(B)); MinB = min(min(B));
    nB = B - MinB;
    nnB = nB/(MaxB-MinB);
    
    % vkládání segmentů (sloupců) do spektrogramu:
    RGB(:,pom,1) = nnR;
    RGB(:,pom,2) = nnG;
    RGB(:,pom,3) = nnB;
                          
    pom = pom+1; % přechod na další segment
end

%% vykreslení spektrogramu
image(RGB)
periody = round((w./(1:w/2))*100)/100; % převedení svislé osy z frekvence na periodu
krok_osy_y = 30;
y_tick = [linspace(1,w/2,w/2)];
set(gca,'ytick',y_tick(1:krok_osy_y:end),'yTicklabel',periody(1:krok_osy_y:end))

[rows colms layers] = size(RGB); % upravení vodorovné osy zobrazující pozice bází
krok_osy_x = 10;
osa_x = fix(linspace(1,L,krok_osy_x));
set(gca,'xtick',[linspace(1,colms,krok_osy_x)],'xTicklabel',osa_x)

title('randomseq repetice CGG') % název sekvence
xlabel('pozice v sekvenci')
ylabel('perioda')
