clc; %Clear command window
clear all;
t0 = clock; %Current date and time as date vector
mkdir dataDC; %create fold storing data files

%计算数据集里的蛋白质数目(即蛋白质网络中的节点数)
allorthology='Nucleus_protein_allorthology.txt';
DATAfile =allorthology;
fid = fopen(DATAfile,'r'); %open the file 
Df = textscan(fid,'%s%s%d%s'); %read file
fclose(fid); %close file, success 0, fail -1 
protein=cat(2,Df{1},Df{2},Df{4});
essentialproteinnumbers = sum(Df{3});  %关键蛋白质数目

%打开亚细胞定位的数据文件
subcell='subcellular.txt';
DATAfile_s = subcell;
fid_s=fopen(DATAfile_s,'r');
Df_s = textscan(fid_s,'%f');
fclose(fid_s);
sub=Df_s{1}; 



%打开Nucleus_PPI.txt文件，计算蛋白质总数
PIN='Nucleus_PPI.txt';
DATAfile_r=PIN;
fid_r = fopen(DATAfile_r,'r'); %open file
Df_r = textscan(fid_r,'%s%s'); %read file
fclose(fid_r); %close file success 0,fail -1
Pro = union(Df_r{1},Df_r{2}); %用union取并集，Pro存放网络中的所有蛋白质
number = length(Pro); %number是蛋白质总数


%创建邻接矩阵
[bool,edge(:,1)] = ismember(Df_r{1},Pro); %edge以蛋白质编号的形式存放网络中的所有边
[bool,edge(:,2)] = ismember(Df_r{2},Pro);
AdjMatrix = zeros(number); %AdjMatrix为邻接矩阵
for i = 1:length(Df_r{1}) %length(Df_r{1})是网络的边数
    AdjMatrix(edge(i,1),edge(i,2)) = 1;
    AdjMatrix(edge(i,2),edge(i,1)) = 1;
end
SparseAdjMatrix = sparse(AdjMatrix); %sparse求稀疏矩阵，用于保存，节省空间
save dataDC/AdjMatrix.mat SparseAdjMatrix;

%计算节点度dc
degree = sum(AdjMatrix,2); %对AdjMatrix的每行求和，degree存放每个节点的度
%maxdegree = max(degree); % 求degree的最大值

%计算节点度dc的值
%dc = degree/maxdegree;
%dc= dc+sub;
%dc=dc/max(dc);




%计算ic
%B = diag(degree)+ones(number)-AdjMatrix; %diag(degree)是以degree为对角线元素，其他元素为零的矩阵
%M = pinv(B); %pinv求广义逆矩阵
%ic = zeros(number,1); %ic存放每个节点的ic值
%for i = 1:number
%    for j = 1:number
%        if j ~= i
%            ic(i) = (M(i,i)+M(j,j)-M(i,j)) + ic(i); %计算ic值
%        end
%    end
%end
%ic = number./ic;
%计算bc
%[bc,E] = betweenness_centrality(SparseAdjMatrix); %matlabBGL中的betweenness_centrality函数；bc返回点介数，E返回边介数
%bc=bc/max(bc);

%计算 SOECC
ECC= zeros(number,number);  %ECC为边聚集系数矩阵

for i = 1:number
    for j = 1:number
    neighbornumber=0;
     if ( AdjMatrix(i,j) ~=0)
      for k=1:number
       if(k~=j)&&(k~=i)&&(AdjMatrix(i,k)~=0)&&(AdjMatrix(j,k)~=0)
       neighbornumber=neighbornumber+1;
       end
      end
      if((degree(i)>1) && (degree(j)>1))
      ECC(i,j)=neighbornumber/min(sum(AdjMatrix(i,:))-1 ,sum(AdjMatrix(j,:))-1);
      end
     end 
    end
   end
SOECC = sum(ECC,2); %对ECC的每行求和  
ecc=SOECC/max(SOECC);







%计算ec
%[V,D] = eig(AdjMatrix); %对邻接矩阵求其特征值和特征向量；对角矩阵D的对角线上元素是特征值，并且是由小到大排列的，矩阵V的每一列是对应于特征值的特征向量
%if sum(V(:,number)<0) < sum(V(:,number)>0)
%   ec = V(:,number); %最大的特征值对应的特征向量，即为ec；由于特征值是由小到大排列的，所以最大的特征值在最后一列，即第num列
%else
%   ec = -V(:,number); %如果特征向量的num个分量中绝大部分都是负值(负值个数大于正值个数)，则对特征向量取负值(即乘以-1)仍然是对应于相应特征值的特征向量，若不取负变正的话，结果很糟糕
%end

%计算sc
%sc = zeros(number,1); %sc存放每个节点的sc值
%for i = 1:number
%    for j = 1:number
%        sc(i) = exp(D(j,j))*(V(i,j))^2 + sc(i); %计算sc值
%    end
%end
%sc=sc/max(sc);

%综合DC和SC的信息
value=zeros(number:11);
for i=0:10
    value(:,i+1)=i/10*ecc+(1-i/10)*sub;
end
 

ratio=99; %99 reference organisms
alpha=0.7;




tNodes= zeros(number,1); 
for i=1:number
    tNodes(i)=(1-alpha)*str2num(protein{i,3})/sum(ratio);
end

H_0 = zeros(number,11);
for j=1:11
  for i = 1:number
        H_0(i,j)=alpha*value(i,j);
  end
end

%计算score 
score = zeros(number,11);
for j=1:11
  for i=1:number
    score(i,j) = tNodes(i)+H_0(i,j);
  end
end

%比较结果
ESpro = importdata('Essential_yeast.txt');
filter =[51,255,510,764,1019,1274];
%filter =[100,200,300,400,500,600]; % 不同的输入数据, 数值不同
ESnum = zeros(6,11);
value1 = zeros(number,11); 
inx = zeros(number,11);
params_0 = score;

for i=1:11
  [value1(:,i),inx(:,i)] = sort(score(:,i),'descend'); %按分值按降序排列
end
%分别取降序排列的前51，255，510，764，1019，1274个蛋白质统计其中的关键蛋白质数目并存放在ESnum的第一列中
for j=1:11
  for i = 1:6
      f_num = filter(i);
      i_num = inx(1:f_num,j);
      p_num = protein(inx(1:filter(i),j),2);
      sum1 = 0;
      for x = 1:length(p_num)
          [bool,is_num]=ismember(p_num(x),ESpro);
          sum1=sum1+bool;
      end
      ESnum(i,j)=sum1;  
      
    %ESnum(i,j) = sum(ismember(protein(inx(1:filter(i),j),2),ESpro)); 
  end
end

essentialproteinnumbersum=0;
for j=1:11
  essentialproteinnumbersum=sum(ESnum(:,j));
  if essentialproteinnumbersum>2186
       fprintf('good result %d\n', j);
  end
end



%计算程序运行时间
elapsed_time = etime(clock,t0);
save dataDC/elapsed_time.mat elapsed_time;

%输出到文件		
save dataDC/proteinscore.mat score;
save dataDC/essentialproteinsnumber.mat ESnum;
%OUT=strcat('locs/',locs,'_son.mat'); 
%save (OUT,'results')
%save Data/OUT results '-ascii' '-tab';