E=pcreff(ng,cp)

function E=pcreff(ng,cp)
%PCREFF Set the Efficiency of a RT-PCR to use in the relative
%quantification of transcripts.
%
%Reverse transcription(RT) followed by PCR is a powerful tool for the
%detection and quantification of mRNA. It is the most sensitive method for
%the detection and quantification of gene expression levels, in particular
%for low abundance mRNA.
%The relative quantification is based on the expression ratio of a target
%gene versus a reference gene. Some mathematical models have already been
%developed to calculate the relative expression ratios, with or without
%efficiency correction. Normally the PCR efficiency is set at 2 (the max 
%possible value) for the reference and target gene, but a difference in 
%PCR efficiency of 0.03 between the target and reference gene, the falsely
%calculated difference in expression ratio is 46% in case of Et<Er and 209%
%in the case of Et>Er. The difference will increase dramatically by higher
%efficiency differences: i.e. DE=0.05 (27% and 338%) and DE=0.1 (7.2% and
%1083%)
%This function computes the efficiency of PCR reaction and is based on MYREGR 
%function. If it is not present on the computer, pcreff will try to download
%it from FEX
%
% Syntax: 	pcreff(ng,cp)
%      
%     Inputs:
%           NG - Array of the ng of mRNA or cDNA used 
%           CP - Crossing points data. These data can be inserted as an
%           array or as a matrix (2 or more replicates for each ng data).
%           In the last case the mean will be calculated
%     Outputs:
%           - Summary of MYREGR function
%           - PCR efficiency and error.
% Example:
%       ng=[30 15 3 0.5 0.12];
%       cp=[19.120 20.230 20.070; 20.770 20.780 20.860; 23.890 23.900
%       23.770; 26.650 26.540 26.680; 29.510 29.590 29.680];
%
%   Calling on Matlab the function: 
%             pcreff(ng,cp)
%
%   Answer is:
%
% (...) All the outputs of MYREGR function + calibration plot
%                    PCR Efficiency
% ------------------------------------------------------------
%      Value (%)               S.E.
% ------------------------------------------------------------
%    1.76494 (76.49421)        0.02726
% ------------------------------------------------------------
%
%           Created by Giuseppe Cardillo
%           [e-mail]
%
% To cite this file, this would be an appropriate format:
% Cardillo G. (2008) PCREfficiency: set the Efficiency of a RT-PCR to use
% in the relative quantification of transcripts.
% http://www.mathworks.com/matlabcentral/fileexchange/20887

%Input errors handling
if nargin~=2
    error('Two input arguments are required')
end

if isscalar(ng) || isscalar(cp)
    error('ng must be a Cx1 vector. Cp must be a Cx1 vector or a CxN matrix.')
end

if ~isvector(ng)
    error('ng must be a vector.')
else
    ng=ng(:);
end

if isvector(cp)
    cp=cp(:);
else
    cp=mean(cp,2);
end

if size(ng,1)~=size(cp,1)    
    error('ng and Cp must have the same number of rows')
end

try
    slope=myregr(log(ng),cp);
    A.value=exp(-1/slope.value); %PCR efficiency (1<=A<=2)
    A.err=1/slope.value^2*exp(-1/slope.value)*slope.se; %Error propagation

%Display the results
    tr=repmat('-',1,60);
    disp(' ')
    disp('                   PCR Efficiency')
    disp(tr)
    disp('     Value (%)               S.E.')
    disp(tr)
    fprintf('%10.5f (%0.5f)     %10.5f\n',A.value,(A.value-1)*100,A.err)
    disp(tr)

    if nargout
        E=A;
    end
catch ME
    disp(ME)
    disp('I am trying to download the myregr function by Giuseppe Cardillo from FEX')
    [F,Status]=urlwrite('http://www.mathworks.com/matlabcentral/fileexchange/15473-myregression?controller=file_infos&download=true','myregr.zip')
    if Status
        slope=myregr(log(ng),cp);
        A.value=exp(-1/slope.value); %PCR efficiency (1<=A<=2)
        A.err=1/slope.value^2*exp(-1/slope.value)*slope.se; %Error propagation

        %Display the results
        tr=repmat('-',1,60);
        disp(' ')
        disp('                   PCR Efficiency')
        disp(tr)
        disp('     Value (%)               S.E.')
        disp(tr)
        fprintf('%10.5f (%0.5f)     %10.5f\n',A.value,(A.value-1)*100,A.err)
        disp(tr)

        if nargout
            E=A;
        end
    end
    clear F Status
end