% ======================================================================================================================================= %
%                                                  The Relative R-Squared Method (RRSM)                                                   % 
% ======================================================================================================================================= %
% (1) Choose a "p0" value and a "s" value as thresholds for p-value and the relative R^2 value,                                           %
%     respectively, for runing RRSM.                                                                                                      %
%     eg. p0=0.5 , s=0.99                                                                                                                 %
% (2) There are two options for the data format, normalized data and original data.                                                       %
%     Set ToF=1 for nomalize data and ToF=0 for original data                                                                             %
% (3) You need to upload the data first                                                                                                   %
% >    C=as.matrix(read.table('D:/DATA/C.xls'))                                                                                           %  
% >    X=as.matrix(read.table('D:/DATA/X.xls'))                                                                                           % 
% >    Z=as.matrix(read.table('D:/DATA/Z.xls'))                                                                                           %
% >    glbl_HGNC=as.matrix(read.table('D:/DATA/glbl_HGNC.txt'))                                                                           %
% >    mlbl=as.matrix(read.table('D:/DATA/mlbl.txt'))                                                                                     %
%                	                                                                                                                      %               
% (4) Save the codes and data in the same folder                                                                                          % 
%                                                                                                                                         %
%     Run the codes:                                                                                                                      %
%     eg. [finalSelectName,numtotal] = RRSM(0.5,0.99,0,C,X,Z,glbl_HGNC,mlbl)                                                              %
% ========================================================================================================================================%

function [finalSelectName,numtotal] = RRSM(p0,s,ToF,C,X,Z,glbl_HGNC,mlbl)  

     mRNAnum=size(C,1);   % Number of mRNAs 
     miRNAnum=size(C,2);  % Number of microRNAs 
     tissuenum=size(X,2); % Number of tissue/cell types 

     nonzero_mRNA=find(sum(C,2)); % we can know which mRNA being targeted 
     nonzero_mRNAnum=size(nonzero_mRNA,1);

     rr=[];      % save relative R^2 value
     rsquare1=[];% save the original R^2 value
     rsquare2=[];% After regressing the mRNA on the selected miRNAs by 1st criterion p-value, save the R^2value. 
     
     if ToF==1
         L1=ones(mRNAnum,1);
         Xmean=mean(X);                      % column mean of mRNA*tissues expression  (1*tissuenum)
         Xstd=std(X);                        % column standard deviation of mRNA*tissues expression (1*tissuenum)
         Xnormalize=(X-L1*Xmean)./(L1*Xstd); % standardize mRNA*tissues expression  (mRNAnum*tissuenum)

         Zmean=mean(Z);                      % column mean of miRNA*tissues expression (1*tissuenum)
         Zstd=std(Z);                        % column standard deviation of miRNA*tissues expression (1*tissuenum) 
         L2=ones(114,1);
         Znormalize=(Z-L2*Zmean)./(L2*Zstd); % standardize miRNA*tissues expression (miRNAnum*tissuenum)
         X=Xnormalize;
         Z=Znormalize;
     else
         X=X;
         Z=Z;
     end

     Selectrue=zeros(mRNAnum,miRNAnum);

     for i=1:nonzero_mRNAnum % just for loop nonzero_mRNAnum(890), not all mRNA number(16063)
   
       INDEPEND=[];      % save temporary miRNA across tissues corresponding their mRNA   
       miRvalue=[];      % note corresponding numbers of miRNAs   
       count1=0;count2=0;

       for k=1:miRNAnum        
           if C(nonzero_mRNA(i),k)==1           
               INDEPEND=[INDEPEND;Z(k,:)];          
               miRvalue=[miRvalue;k];            
               count1=count1+1;        
           else
               INDEPEND=INDEPEND;        
           end

       end

       YY=(X(nonzero_mRNA(i),:))';  
       V=(INDEPEND)';    
       S=V'*V;   
       b=inv(S)*V'*YY; % using least square error method to estimate the regression parameters  
       YYbar=mean(YY); 
        
       %SST=(YY-YYbar*I)'*(YY-YYbar*I);    
       SST=YY'*YY;              % total sum of squares (No intercept here)   
       SSE1=(YY-V*b)'*(YY-V*b); % residual sum of squares
      
       rsquare1(i)=1-SSE1/SST; % solve the R^2 value
      
       inn=2*(1-normcdf(abs(b./(diag(inv(S))).^(1/2)))); % calculate each b's p-value   
       rlabel=inn < p0;                                  % label which b's p-value < threshold
        
       INDEPEND2=[]; % After applying the measurement of p-value, save temporary miRNA across tissues corresponding their mRNA   
       miRvalue2=[]; % After applying the measurement of p-value, note corresponding numbers of miRNAs
        
       if rlabel==0        
           rsquare2(i)=0;  
       else
           miRvalue=miRvalue.*rlabel;     
        
           for j=1:count1           
               if miRvalue(j)~=0           
                   count2=count2+1;          
                   miRvalue2(count2)=miRvalue(j);         
                   INDEPEND2=[INDEPEND2;Z(miRvalue2(count2),:)];               
               end
           end
           
           V2=(INDEPEND2)';        
           S2=V2'*V2;        
           b2=inv(S2)*V2'*YY;      
           SSE2=(YY-V2*b2)'*(YY-V2*b2);
             
           rsquare2(i)=1-SSE2/SST;   % After applying the measurement of p-value, solve the R^2 value     
       end
       rr(i)=rsquare2(i)/rsquare1(i);% relative R^2 value
       
       if rr(i) < s       
           Selectrue=Selectrue;    
       else
           for h=1:count1           
               if rlabel(h)==1               
                   Selectrue(nonzero_mRNA(i), miRvalue(h))=1;            
               else
                   Selectrue=Selectrue;           
               end
           end
       end
end

rr=rr';
rsquare1= rsquare1';
rsquare2= rsquare2';

numtotal=sum(Selectrue(:)); %Sum all elements of Selectrue to show how many high-confidence targets being selected?(mRNAnum*miRNAnum)
[a b]=find(Selectrue);
finalSelect=[b,a];          % Figure out the code names of miRNA vs. mRNA 
finalSelectName=[mlbl(b) glbl_HGNC(a)]; % return the miRNAs and their high-confidence targets

