clear all
load BW38.mat


GeneNamesCell=cellstr(GeneNames);
%%%%%%usinięcie house keeping genes
AFFX='AFFX'
for i=1:length(GeneNamesCell)
   X=strncmp(AFFX,GeneNamesCell,4);
end

empty=find(X);
DNAdata(empty,:)=[];
GeneNamesCell(empty,:)=[];
%%%%%%normalizacja z-score
normalDNAdata=zscore(DNAdata);

TN=0; TP=0; FP=0; FN=0;
%%%%%%podział na zbiór testowy (1) i uczący(37)

for i=1:size(DNAdata,2)
    test=normalDNAdata(:,i);  %dane testowego
    probka=SampleNames(i,:); %zbiór testowy
    gr_test=probka(:,4:6); %grupa zbioru testowego
    uczacy=normalDNAdata;
    uczacy(:,i)=[];
    %selekcja
[mask,data,gene] = geneentropyfilter(uczacy,GeneNamesCell);
testowe= test(mask,:);
[mask2, data2, gene2] = generangefilter(data,gene,'ABSVALUE',log2(4));
testowe= testowe(mask2,:);

ALL='ALL';
ALM='ALM';

sample=SampleNames(:,4:6);
sample=cellstr(sample);

for j=1:38
   ALL = strcmp(sample,'ALL');
   AML = strcmp(sample,'AML');
end

DNAdataALL=normalDNAdata(:,ALL);
GeneNamesALL=SampleNames(ALL,:);

DNAdataAML=normalDNAdata(:,AML);
GeneNamesAML=SampleNames(AML,:);

       %test t-studenta
[pvaluesP,tscoresP]=mattest(DNAdataALL,DNAdataAML,...
                    'Permute', 100);

 Gene005=find(pvaluesP<0.05);
 Gene005names=GeneNamesCell(Gene005);



volcano = mavolcanoplot(DNAdataALL,DNAdataAML,pvaluesP,'Labels', GeneNamesCell);
close all
Genes_vol=0;
Genes_vol = getfield(volcano,'GeneLabels');
for i=1:length(Genes_vol),
    HKG2= strcmp(Genes_vol{i},GeneNamesCell);
    [row(i),v(i)]=find(HKG2==1);
end
row=row';
values=DNAdata(row,:);

testowe=test(row,:);


%%%%klasteryzacja hierarchiczna

[cidx, ctrs] = kmeans(values, 2,...
                      'dist','corr',...
                      'rep',5,...
                      'disp','final');

[row_a,Kl1]= find(cidx==1); %wiersze z genami należącymi do 1-szego klastra
Genes_Kl1=Genes_vol(row_a);
data_Kl1=values(row_a);

[row_b,Kl2]= find(cidx==2); %wiersze z genami należącymi do 2-go klastra
Genes_Kl2=Genes_vol(row_b);
data_Kl2=values(row_b);


%klasyfikacja

a=test(row_a,:); % dane zbioru testowego dla genów z klasy ALL
b=test(row_b,:); %AML

dif1=(a-data_Kl1);    % odjęcie wartości ekspresji genów 1 (ALL) klastra od wartości ekspresji genów zbioru testowego
dif2=(b-data_Kl2);    %odjęcie wartości ekspresji genów 2 (AML) klastra od wartości ekspresji genów zbioru testowego

if mean(dif1)<mean(dif2)
    wynik='ALL';     %przypisanie grupy ALL do zbioru testowego
    if strncmp(gr_test,wynik,3)==1    %wynik zgadza się z prawdziwą przynależnością zbioru
    TN(i)=1;
    else FN(i)=1;      %błędne dopasowane
    end
else
    wynik2='AML';   %przypisanie grupy AML do zbioru testowego
    if strncmp(gr_test,wynik2,3)==1
    TP(i)=1;
     else FP(i)=1;
    end
end



end

TN=sum(TN);
TP=sum(TP);
FP=sum(FP);
FN=sum(FN);

I = ((TP + TN)/(TP + TN+ FP + FN))*100; %jakość
SE = TP/(TP+FN);       %czułość
SP = TN/(TN+FP);  %specyficzność