function [class, normpvalue, predlr] = DrugPredictionScore(training, test, groups, numbern)
% This function takes as an argument a (training) set that has
% observations(drugs)as rows and variables as columns, (test) is a test vector to classify using
% training and in this version is limited to n=1. (groups) is a column vector
% of group memberships(that are sequential and numeric starting at n=1).  Note if you have 7 groups they must be ordered 1-7 for 
% the observations in (training). Numbern is the
% number of nearest neighbors to include, the default is majority rule.  The class is the predicted group
% mebership in test, and pvalue is the predicted significance of the
% inclusion of the test compound in the predicted class
% Currently a KNN classifier is hardcoded into the algorithm, but the
% algorithm can be used with any matlab multi-class classifier if the code
% is changed.  In our shRNA dataset we have noted that Naive-Bayes
% Classifiers perform similarly well.  In either case a euclidean metric of
% cluster size is used to assess significance.

%Step1 creating the reference distances and null distributions
%__________________________________________________________________________
%__________________________________________________________________________
DistanceMatrix=squareform(pdist(training,'euclidean'));
M=max(groups)
i=1
nulldist={}
ALtot=[]
while i< M+1
    j=1%makes an index that terminates when the max group number is reached
    AL=[]
    gpos=[]
    nongpos=[]
    while j< length(groups)+1  %looks through the groups vector sequentially to determine which entries match the group label i
        if groups(j,1)==i %sequentially checks the equality
            gpos=horzcat(gpos,j)%horizontally concatenates the reference row for a given group.  
        elseif groups(j,1)~=i
            nongpos=horzcat(nongpos,j)
        end
        j=j+1
    end
    S=length(gpos)
    Sn=length(nongpos)
    submati=[DistanceMatrix(gpos(1,1):gpos(1,S),gpos(1,1):gpos(1,S))]%builds submatrix of group specific distances and then calculates the averages on the basis of this
    AL=sum(sum(submati))/(S*(S-1))
    ALtot=horzcat(ALtot,AL)
    k=1
    newalvec=[]
    while k < Sn+1 %this step build null distributions based on the empircal prediction using false negatives
        alldistances=DistanceMatrix(nongpos(1,k),:)
        l=1
        nullpredictionsdist=[]
        while l < S+1
           nullpredictionsdist=vertcat(nullpredictionsdist,alldistances(1,gpos(1,l)))
           l=l+1
        end
        newal=(sum(sum(submati))+sum(nullpredictionsdist))/((S*(S-1))+S)
        newalvec=vertcat(newalvec,newal)
        k=k+1
    end
    storednull{1,i}=newalvec     
    i=i+1
end

%__________________________________________________________________________
%step2 is the prediction of drug class and updating a new training set to
%take into account that prediction, and then calculating the average
%linkage following the inclusion of the test drug
%________________________________________________________________________
training2=training
test2=test
groups2=groups
class=knnclassify(test, training, groups, numbern, 'euclidean')
m=1
while m < length(class)+1
    n=1
    while n < length(groups)+1
        if groups(n,1)==class(m,1)
            if n<length(groups)
                trainingnew=vertcat(training2(1:n,:),test2(m,:),training2(n+1,:))
                groupsnew=vertcat(groups2(1:n,:),class(m,:),groups2(n+1,:))
            elseif n==length(groups)
                trainingnew=vertcat(training2(1:n,:),test2(m,:))
                groupsnew=vertcat(groups2(1:n,:),class(m,:))
            end  
        end
        n=n+1
    end
    m=m+1
end
DistanceMatrixnew=squareform(pdist(trainingnew,'euclidean'))
Mnew=max(groupsnew)
i=1
ALnew=[]
while i< Mnew+1
    j=1%makes an index that terminates when the max group number is reached
    AL2=[]
    gposnew=[]
    nongposnew=[]
    while j< length(groupsnew)+1  %looks through the groups vector sequentially to determine which entries match the group label i
        if groupsnew(j,1)==i %sequentially checks the equality
            gposnew=horzcat(gposnew,j)%horizontally concatenates the reference row for a given group.  
        elseif groupsnew(j,1)~=i
            nongposnew=horzcat(nongposnew,j)
        end
        j=j+1
    end
    Snew=length(gposnew)
    Snnew=length(nongposnew)
    submatinew=[DistanceMatrixnew(gposnew(1,1):gposnew(1,Snew),gposnew(1,1):gposnew(1,Snew))]%builds submatrix of group specific distances and then calculates the averages on the basis of this
    AL2=sum(sum(submatinew))/(Snew*(Snew-1))
    ALnew=horzcat(ALnew,AL2)
    i=i+1
end

%__________________________________________________________________________
%_Step3 is comparing to the background distribution
%1st calculate for the group predicted in class, what is the distribution of
%linkage ratios for the given group, calculate from altot(yourclass) and
%the classes entry in storednull what the distribution of increases for
%null(false) predictions are.  Then calculate the linkage ratio alnew/altot
%then get pvalue
%_________________________________________________________________________
lrdist=[]
p=1
while p < length(storednull{1,class})+1 
    lrdist=vertcat(lrdist,(storednull{1,class}(p))/ALtot(1,class))
    p=p+1
end
predlr=ALnew(1,class)/(ALtot(1,class))
lr=lrdist
lr = lr(:);

% Set up figure to receive datasets and fits
f_ = clf;
figure(f_);
set(f_,'Units','Pixels','Position',[547.333 181 680 475.45]);
legh_ = []; legt_ = {};   % handles and text for legend
ax_ = newplot;
set(ax_,'Box','on');
hold on;

% --- Plot data originally in dataset "lr data"
t_ = ~isnan(lr);
Data_ = lr(t_);
[Y_,X_] = ecdf(Data_,'Function','cdf'...
              );  % compute empirical function
h_ = stairs(X_,Y_);
set(h_,'Color',[0.333333 0 0.666667],'LineStyle','-', 'LineWidth',1);
xlabel('Data');
ylabel('Cumulative probability')
legh_(end+1) = h_;
legt_{end+1} = 'lr data';

% Nudge axis limits beyond data limits
xlim_ = get(ax_,'XLim');
if all(isfinite(xlim_))
   xlim_ = xlim_ + [-1 1] * 0.01 * diff(xlim_);
   set(ax_,'XLim',xlim_)
end

x_ = linspace(xlim_(1),xlim_(2),100);
%______________________________

% --- Create fit "fit 1"
%____________________________

% Fit this distribution to get parameter values
t_ = ~isnan(lr);
Data_ = lr(t_);
% To use parameter estimates from the original fit:
%     p_ = [ 5.5, 3.027650354097];
pargs_ = cell(1,2);
[pargs_{:}] = normfit(Data_, 0.05);
p_ = [pargs_{:}];
y_ = normcdf(x_,p_(1), p_(2)); % compute cdf
h_ = plot(x_,y_,'Color',[1 0 0],...
          'LineStyle','-', 'LineWidth',2,...
          'Marker','none', 'MarkerSize',6);
legh_(end+1) = h_;
legt_{end+1} = 'fit 1';

hold off;
leginfo_ = {'Orientation', 'vertical', 'Location', 'NorthWest'}; 
h_ = legend(ax_,legh_,legt_,leginfo_{:});  % create legend
set(h_,'Interpreter','none');
normpvalue=normcdf(predlr,p_(1), p_(2))
%this works