%choose App.1 or App.2
%list order from small accession number to large accession number
%Application 1 : 9 species bacteria
%list={'NC_000853'; 'NC_000964'; 'NC_002754'; 'NC_002951'; 'NC_003997';'NC_004193';'NC_004557';'NC_005945';'NC_007929'};
%per : Genomes after the first one are listed in the order of: same species as the first two genome, same genius (but different species),
%same family (but different genius), same order (but different family), same class (but different order), same phylum (but different class), 
%same kingdom (but different phylum), respectively.
%per={'NC_003997';'NC_005945';'NC_000964';'NC_004193';'NC_002951';'NC_007929';'NC_004557';'NC_000853';'NC_002754'};
%Application 2 : Mycoplasma group
%list={'NC_000908';'NC_000912';'NC_002771';'NC_004432';'NC_004829';'NC_005364';'NC_006360';'NC_006908';'NC_007294';'NC_007295';'NC_007332';'NC_007633';'NC_009497'};
%per : order from small genome to large genome
%per={'NC_000908','NC_006908','NC_007294','NC_000912','NC_009497','NC_006360','NC_007295','NC_007332','NC_002771','NC_004829','NC_007633','NC_005364','NC_004432'};
list={'NC_000853'; 'NC_000964'; 'NC_002754'; 'NC_002951'; 'NC_003997';'NC_004193';'NC_004557';'NC_005945';'NC_007929'};
per={'NC_003997';'NC_005945';'NC_000964';'NC_004193';'NC_002951';'NC_007929';'NC_004557';'NC_000853';'NC_002754'};
cd('C:\\BGVM')
par=load('MLE.txt');
[pn,le]=size(list); 
fre1=zeros(36,1);
fre2=zeros(36,1);
hotspots=cell(pn,pn);
for k = 1: (pn-1) 
	for j = (k+1):pn
        i=(k-1)*pn-k*(k-1)/2+(j-k); 
        r=par(11,i);
        c=par(12,i);
        data=dlmread([list{k,1},'&',list{j,1},'_sin.txt'],'\t'); 
  		id1=find(strcmp(list(k,1),per));
        id2=find(strcmp(list(j,1),per));
        if id1<id2
            index1=id1;
            index2=id2;
            a=data(:,1);
            b=data(:,2);
        else
            index1=id2;
            index2=id1;
            b=data(:,1);
            a=data(:,2);
        end	
  		for m = 1:36			
   			fre1(m)=length(find(10*m-10<a&a<10*m));
   			fre2(m)=length(find(10*m-10<b&b<10*m));
        end
        r1=strread(num2str(par(11,i)),'%s');
		c1=strread(num2str(par(12,i)),'%s');
        angle1=load('angle.txt');
        data1=[];
        for i=1:length(angle1)
            data1=[data1;(angle1(i)/180*pi)*ones(fre1(i),1)];
        end
        angle2=load('angle.txt');        
        data2=[];
        for i=1:length(angle2)
            data2=[data2;(angle2(i)/180*pi)*ones(fre2(i),1)];
        end
        if max(max(fre1),max(fre2))>20
            figure;
            subplot(1,2,1);
            polar(0:2*pi/(360-1):2*pi,40*ones([1 360]),':k');
            hold on;
            polar(0:2*pi/(360-1):2*pi,20*ones([1 360]),':k');
            rose(data1,36);
            view(90,-90);
            subplot(1,2,2);
            polar(0:2*pi/(360-1):2*pi,40*ones([1 360]),':k');
            hold on;
            polar(0:2*pi/(360-1):2*pi,20*ones([1 360]),':k');
            h=rose(data2,36);
            set(rose(data2,36),'LineWidth',0.5,'Color','red');
            view(90,-90);           
        elseif max(max(fre1),max(fre2))<=20 && max(max(fre1),max(fre2))>10
            h=figure;
            h=subplot(1,2,1);
            h=polar(0:2*pi/(360-1):2*pi,20*ones([1 360]),':k');
            hold on;
            h=polar(0:2*pi/(360-1):2*pi,10*ones([1 360]),':k');
            h=rose(data1,36);
            view(90,-90);
            h=subplot(1,2,2);
            h=polar(0:2*pi/(360-1):2*pi,20*ones([1 360]),':k');
            hold on;
            h=polar(0:2*pi/(360-1):2*pi,10*ones([1 360]),':k');
            h=rose(data2,36);
            set(rose(data2,36),'LineWidth',0.5,'Color','red');
            view(90,-90);            
        elseif max(max(fre1),max(fre2))<=10
            h=figure;
            h=subplot(1,2,1);
            h=polar(0:2*pi/(360-1):2*pi,10*ones([1 360]),':k');
            hold on;
            h=polar(0:2*pi/(360-1):2*pi,5*ones([1 360]),':k');
            h=rose(data1,36);
            view(90,-90);
            h=subplot(1,2,2);
            h=polar(0:2*pi/(360-1):2*pi,10*ones([1 360]),':k');
            hold on;
            polar(0:2*pi/(360-1):2*pi,5*ones([1 360]),':k');
            h=rose(data2,36);
            set(rose(data2,36),'LineWidth',0.5,'Color','red');
            view(90,-90);            
        end
        saveas(h,[r1{1,1},'_',c1{1,1},'.fig'])
        clc;
        close;        		
	end
end	
clear all;