%This script is associated with Bisarad et al.
%It reads in the data from Averaged_Array_Data.xlsx

clear
close all

%choose which mAb you want to model
%choices are Ab1, Ab8, 4C1, DM1A, Lnkb2, TF3B5, C3 and 9E10.3
current_mAb='DM1A';
Data_file='Averaged_Array_Data.xlsx'; %This contains all the mAb averaged data

num_rep=12;%number of times it repeats the training and prediction to average for each run
numrand=1000000;%number of random sequences used in the ranking test


%Main Parameter set*******************************************************
%Network parameters
num_layers=2;
num_nodes=250;
miniBatchSize = 1000;
MaxEpochs=20;
InitialLearnRate=0.002; %originally 0.002
LearnRateDropFactor=0.9; %originally 0.9
%
%Perameters associated with shifting sequences
shift=5; %0 means no shifting

%Parameters associated with weighting sequences
weight_high_binders=true(1); %Options to weight the high binders
maxdensity=300;%this is the number of copies desired in each bin for the max version
maxcopy=100;%this is the maximun number of copies for any particular peptide in either version
conc2compare=8;%this is the concentration it looks at to determine how much to weight the high binders (in nM)

%endmarkers -- parentheses that mark the N and C terminus
EndMarkers=true(1);

%extrapolate saturated points based on lower concentrations
adjust_saturated_points=true(1);

%subtract compositional binding
Composition_fit=true(1);

%exclude the really long sequences
ExcludeLong=true(1);
maxlen=11;%if we exclude the long sequences, this is the max length
%************************************************************************

%Parameters associated with sequences
aminos='ADEFGHKLNPQRSVWYX';%these are the amino acids used in array sequences

ifig=0;

%The following are the cognate epitopes that are on the array and the
%number of times they are repeated.  We will average them before doing
%analysis so there is only one value and we will remove that value when
%doing neural network training.  Note that there are no cognate epitope sequences
%on the array for TF3B5, C3 or 9E10.3.

% %
%    mAb    epitope_seq      num_repeat on array
%    ___    ______________    __________
%
%     Ab1   {'RHSVV'     }       210
%     Ab1   {'RHSVVV'    }       100
%     Ab8   {'SDLWKL'    }       245      SVEPPLSQETFSDLWKLLPENNV
%     Ab8   {'SDLWKLL'   }       100
%     4C1   {'FDSHYD'    }       200
%     4C1   {'LQAFDS'    }       100
%     4C1   {'QAFDSH'    }       100
%     DM1A  {'ALEKDYE'   }       280
%     DM1A  {'AALEKDY'   }       100
%     DM1A  {'LEKDYEE'   }       100
%     Lnkb2 {'PLEEVLN'   }       100
%     TF3B5 {'EYLGLD'    }         0
%     C3    {'LPNPEGRYT' }         0
%     9E10.3{'EQKLISEED' }         0

%These are just the epitope sequences actually on the array that were
%removed for training
Array_epitope_sequences=[...
    {'RHSVV'},{'RHSVVV'},{''};...
    {'SDLWKL'},{'SDLWKLL'},{''};...
    {'FDSHYD'},{'LQAFDS'},{'QAFDSH'};...
    {'ALEKDYE'},{'AALEKDY'},{'LEKDYEE'};...
    {'PLEEVLN'},{''},{''};...
    {''},{''},{''};...
    {''},{''},{''};...
    {''},{''},{''}...
    ];
Array_epitope_index=~cellfun( @isempty, Array_epitope_sequences);
Array_episeq=Array_epitope_sequences(Array_epitope_index);%this is a list of epitope sequences on the array

%These are the epitope sequences used in the ranking (in the same format as the ones
%present on the array for consistency)
mAb_epitope_sequences=[...
    {'RHSVVVP'};...
    {'SDLWKLL'};...
    {'QAFDSHY'};...
    {'LEKDYEE'};...
    {'PLEEVLN'};...
    {'PEYLGLD'};...
    {'SLPNPEG'};...
    {'KLVSEED'}...
    ];

num_mAb=length(mAb_epitope_sequences);%num_mAb is the number of different mAbs in the dataset

%******************Read in the Data****************************************
%Read the data from the data file (assuming it is in the same directory)
% load('C:\Users\laserweb\Dropbox (ASU)\Current working folders\Research\papers\mAb1\Build the Paper 3\Step 1 Create the dataset/Data_wo_mAb_Sequences');

%definitions of variables below:
%conc_values - the four values of the concentration in increasing order
%Concentration - the concentration of each column in the Data matrix
%Data - 32 column matrix with the 8 mAb samples in order and each with 4
%   concentration values in order
%Description - a sentence describing the load file
%epitope_data - the 32 column data associated with the epitope cognate sequences
%epitope_sequence - the corresponding sequences of the cognate epitopes on
%   the array
%mAb - the vector containing the mAb name associated with each column of
%   Data
%mAb_list - the list of mAb names in order 'Ab1','Ab8','4C1','DM1A','Lnkb2','TF3B5','C3','9E10.3'

[numbers,text,~]=xlsread(Data_file); %read in numbers and text from the file
%Note, the way that Matlab works, the "numbers" skips the header row but the
%"text" does not.
All_Data=numbers(2:end,2:end);%extract the part that has binding values
All_Sequence=text(3:end,1); %extract the first column that has the Sequences jumping below the header and concentration rows

%Below are the key variables we need going forward from the data file
Concentration=numbers(1,2:end);%The concentrations are in the first row associated with each column
Cognate_index=numbers(2:end,1)>0;%The second column of the numbers is 1 for cognate sequences, 0 if not
mAb=text(1,3:end); %This row has the mAb names associated with each column
mAb_list=[{'Ab1'},{'Ab8'},{'4C1'},{'DM1A'},{'Lnkb2'},{'TF3B5'},{'C3'},{'9E10.3'}]; %this is the list of mAb names
conc_values=unique(Concentration); %this is the list of concentration values
conc_values=conc_values';
Data=All_Data(~Cognate_index,:);%This removes the cognate sequences from the data matrix
Sequence=All_Sequence(~Cognate_index,:);%This removes the cognate sequences from the sequence cell vector

[N,S]=size(Data);

%*************************************************************************

Data=log10(Data);%We will use the log10 of the data for the model
orig_data=Data; %make sure we keep a copy of the data before preparing it for training

num_conc=length(conc_values); %the number of different concentrations

%Set up an index that can be used to grab data from different mAbs
Sample_index=false(num_mAb,S); 
for imAb=1:num_mAb
    Sample_index(imAb,:)=strcmp(mAb,mAb_list(imAb));
end

%decide which mAb we are modeling
ix=1:num_mAb;
imAb=ix(strcmp(mAb_list,current_mAb)); 

%**********Option adjust the saturating points*************************
%This goes through the data and sets saturated points at a value based on
%extrapolating from the next lower concentration and assuming a linear
%dependence on concentration.  Note that the data is set up in order of
%concentration and then mAb, so the first 4 are the four concentrations of
%Ab1 in increasing order, etc.
if adjust_saturated_points
    num_sat=0;
    for c=2:num_conc
        A=orig_data(:,(imAb-1)*4+c);%need to measure saturation in the original data
        B=Data(:,(imAb-1)*4+(c-1));
        %note, 4.8 is approximately saturated (10^4.8 counts)
        Data(A>4.8,(imAb-1)*4+c)=B(A>4.8)+log10(4); %increase saturated data by 4x over previous concentration
        num_sat=num_sat+sum(A>4.8);
    end
    fprintf('for %s the number of saturated data points over all four concentrations was %d\n',mAb_list{imAb},num_sat)
end
%************************************************************************

%Determine the lengths of the sequences
plen=zeros(N,1); %array of lengths of the array sequences
for n=1:N
    plen(n)=length(Sequence{n});
end

%*************Option exclude very long sequences************************
if ExcludeLong
    Data=Data(plen<=maxlen,:); %this removes all data from sequences > maxlen
    Sequence=Sequence(plen<=maxlen); %this removes all sequences > maxlen
    [N,S]=size(Data);  %N and S are the number of sequence-data pairs and the number of different samples (num_conc*num_mAb), respectively
else
    maxlen=max(plen);
end
%*************************************************************************

%Plot the distribution of lengths of sequences used
ifig=ifig+1;
figure(ifig)
plen_range=min(plen):maxlen;
for i=1:maxlen-min(plen)+1 %create the distribution
    plen_dist(i)=sum(plen==plen_range(i));
end
bar(plen_range,plen_dist) %generate a simple bar graph
title('Distribution of the Original Sequence Lengths')
xlabel('Sequence Lengths')
ylabel('Frequency')
drawnow

%*************Option subtract compositional fit*************************
%Calculate a linear composition fit and subtract that binding energy.  A
%vector is calculated with the number of each type of amino acid for each
%sequence.  A set of coefficients and a bias term is calculated relating
%this composition to the binding value for the sequences on the array
if Composition_fit
    Composition_binding=Fit2composition(Data,Sequence,aminos);
    Data=Data-Composition_binding;
end

%*********Option add end markers to the sequences ********************
%This inserts a "(" in front of the N-term residue and a ")" behind the
%C-terminal residue
if EndMarkers
    plen=zeros(N,1);
    for n=1:N
        Sequence(n)=cellstr(['(',Sequence{n},')']);
        plen(n)=length(Sequence{n});
    end
    maxseqlen=max(plen);
    fprintf('Maximum sequence length after adding end markers is %d\n',maxseqlen)
    aminos=[aminos,'(',')']; %add end markers as separate amino acids
end
%************************************************************************

current_mAb(current_mAb<'0')='_'; %there is a '.' in 9E10
fprintf('Starting run for %s. \n',current_mAb);

%Define the data to be used
TrainData=Data(:,Sample_index(imAb,:));%just use the data for that sample
TrainSeq=Sequence;
YTrain=TrainData;

%***************Option Weight the high binders*******************
if weight_high_binders
    ix=1:num_conc;
    iconc=ix(conc_values==conc2compare);%index of concencentration used in determining the weighting
    [Index,Dist]=Weight_constant_density_extended(TrainData,iconc,maxdensity,maxcopy); %Weight function
    TrainData=TrainData(Index,:);%Index specifies which peptides to copy and how many
    TrainSeq=TrainSeq(Index);%repeat for the sequences
    YTrain=TrainData;
end
%******************************************************************

[NN_N,NN_S]=size(TrainData);

%************Option add shifted versions************************
%The function Sequence_Data_Shift returns a 2D matrix of sequences and a 3D
%matrix of data, where the extra dimension is the shifted versions.  It
%then reshapes this into a Sequence vector containing all sequence versions
%and the correspond data 2D matrix
if shift>0
    [ShiftData,ShiftSeq]=Sequence_Data_shift(TrainData,TrainSeq,shift);
    TrainDataShift=reshape(ShiftData,[NN_N*(shift+1),NN_S]);%back to 2D
    TrainSeqShift=reshape(ShiftSeq,[NN_N*(shift+1),1]);%back to a vector
    YTrain=TrainDataShift;
    TrainSeq=TrainSeqShift;
end
%*****************************************************************

%Convert the sequences to a binary representation
XTrain=seq2onehot(TrainSeq,aminos);

%set up a neural network regression using the one hot values as the input
%values and the binding values as the target
numFeatures = size(XTrain,2);
numtargets = size(YTrain,2);

%create the structure of the neural network based on numbers of layers and
%features in the layers.  This is a fully connected, regressor
layers=featureInputLayer(numFeatures,'Normalization', 'zscore'); %set up the input layer
for i_layers=1:num_layers %add the right number of hidden layers and activation functions
    layers=[
        layers
        fullyConnectedLayer(num_nodes)
        reluLayer];
end
layers=[layers
    fullyConnectedLayer(numtargets)
    regressionLayer]; %add the regression layer to the end

%Depending on the machine you are using, you might do this differently.  We
%are running this on a machine with a large number of cpus available but
%only one GPU.  In this case, we run many modeling runs in parallel and use
%the cpus.  However, if you run them one at a time, it is faster on a GPU.
%It defaults to GPU if you remove the "ExecutionEvironment" statement
options = trainingOptions('adam', ...
    'MiniBatchSize',miniBatchSize, ...
    'Shuffle','every-epoch', ...
    'MaxEpochs',MaxEpochs, ...
    'InitialLearnRate',InitialLearnRate,...
    'LearnRateDropFactor',LearnRateDropFactor,...
    'LearnRateSchedule','piecewise',...
    'LearnRateDropPeriod',2,....
    'ExecutionEnvironment','cpu',...
    'Verbose',false);


%Options can also include these, but it will run slower
%     'Plots','training-progress', ...
%     'ValidationData',{XTest,YTest}, ...
%     'ValidationFrequency',30, ...


%fit all the reps and save the models
NN_model_save=cell(1,num_rep); %this is where we store the models

fprintf('Starting fit for %s\n',current_mAb)
fprintf('%s model run with size of XTrain %d x %d and size of YTrain %d x %d\n',...
    current_mAb,size(XTrain,1),size(XTrain,2),size(YTrain,1),size(YTrain,2))
tic
%This loop does the actual training and creates num_rep models.  The reason
%for doing it multiple times is that the neural network starts at random
%places and sometimes that leads it to a non optimal solution.  We normally
%do about 12 of these reps and do them in parallel on a machine with
%sufficent cores and memory to handle that.  Of course you can change
%"parfor" to "for" below and run them one at a time
parfor irep=1:num_rep %this can be made a parfor loop
    fprintf('Running rep #%d\n',irep)
    NN_model = trainNetwork(XTrain,YTrain,layers,options);
    NN_model_save{irep}=NN_model;
end
toc
%at this point we have a set of models

%**********************************************************************
%**********************************************************************
%Now we evaluate the cognate epitope among 1 million random sequences using
%the models developed above

%**************************************************************************
 %Set up the epitope sequences.  Need to pad the sequences add end
    %markers and shift
    padded_epitope_char=repmat('X',[1,maxseqlen]);%pad with X everywhere that there is no sequence
    epitope_cell=mAb_epitope_sequences(imAb);%grab the current epitope sequence
    current_epitope_char=epitope_cell{1};
    if EndMarkers
        padded_epitope_char(1)='(';%put "(" on N-term
        padded_epitope_char(2:length(current_epitope_char)+1)=current_epitope_char;
        padded_epitope_char(length(current_epitope_char)+2)=')';%put ")" on C-term
    else
        padded_epitope_char(1:length(current_epitope_char))=current_epitope_char;
    end
    padded_epitope_cell=cellstr(padded_epitope_char);%convert to cell
    
    %*********************************************************************
    %Now add the shifted versions of the cognate epitope sequence if shift>0
    if shift>0
        %This adds shifted versions to the sequences to be predicted.
        epitopeSeqShift=Sequence_shift(padded_epitope_cell,shift);%function to add shifted versions
        epitopeSeq=epitopeSeqShift';
    else
        epitopeSeq=padded_epitope_cell;
    end
            
    %*********************************************************************
    %Set up the random sequences, pad and shift
    aminos16='ADEFGHKLNPQRSVWY';
    Randseq=repmat('X',[numrand,maxseqlen]);%set up the sequence matrix so it has X where no sequence
    %This next section creates the random sequences with endmarkers if
    %desired.  The random sequences are the same length as the cognate
    if EndMarkers
        Randseq(:,1)='(';%put "(" on N-term
        Randseq(:,2:length(current_epitope_char)+1)=aminos16(randi(16,[numrand,length(current_epitope_char)]));%create random sequences
        Randseq(:,length(current_epitope_char)+2)=')';%put ")" on C-term
    else
        Randseq(:,1:length(current_epitope_char))=aminos16(randi(16,[numrand,length(current_epitope_char)]));%create random sequences
    end
    Randseq=cellstr(Randseq);
    Randseq_no_shift=Randseq;

    %*********************************************************************
    %Next create the shifted versions of the random sequences if desired
    if shift>0
        %This adds shifted versions to the sequences to be predicted.  It
        %then puts them into one long vector
        RandSeqShift=Sequence_shift(Randseq,shift);
        Randseq=reshape(RandSeqShift,[numrand*(shift+1),1]);
    end

    current_mAb(current_mAb<'0')='_';%have to get rid of the '.' in 9E10.3 name

    %*********************************************************************
    %here we start doing the predictions for the num_rep models determined
    %above
    tic
    Pred_Rand=zeros(num_rep,numrand,num_conc);
    Pred_epitope=zeros(num_rep,num_conc);
    
    
    fprintf('Starting prediction for %s\n',current_mAb)
    parfor irep=1:num_rep %this can be made a parfor loop


        NN_model = NN_model_save{irep};%grab a model
        num_residues=NN_model.Layers(1).InputSize/length(aminos); %look at the model and determine the number of residues used
        
        fprintf('Evaluate the random sequence and epitope sequence datasets, rep %d\n',irep)

        %predicte the random sequence data
        %The next function returns a 2D dataset that contains all the
        %shifted sequences
        Pred=NN_Predict(NN_model,Randseq,num_residues,numtargets,miniBatchSize,aminos);
        %the next line then reshapes this so that the shifted predictions
        %of each sequence are in the first dimension
        Pred=reshape(Pred,[shift+1,numrand,num_conc]);
        Pred_Rand(irep,:,:)=squeeze(max(Pred,[],1));%this grabs the max binding prediction for each shifted sequence
        %now work on the epitope sequence
        Pred=NN_Predict(NN_model,epitopeSeq,num_residues,numtargets,miniBatchSize,aminos);
        Pred_epitope(irep,:)=squeeze(max(Pred,[],1)); %Note this has only one epitope so 1 x 4 (as above, grabs the max binder fo the shifted sequences)
    end
    toc
    
    fprintf('calculating means, slopes and ranks for %s\n',current_mAb)
    
    %Calculate mean values for this cycle (Nx4)
    mean_Pred_Rand=squeeze(mean(Pred_Rand,1));%take the mean of the num_rep predictions
    mean_Pred_epitope=squeeze(mean(Pred_epitope,1));%take the mean of the num_rep predictions
    
    
    %Calculate ranks based on max mean binding (Nx1)
    P=[mean_Pred_epitope;mean_Pred_Rand];
    [~,ind1]=sort(P,'descend');
    [~,ind2]=sort(ind1,'ascend');
    ind3=min(ind2,[],2);%this time I take the best ranking concentration for all of them
    epi_rand_maxmean_rank=ind3(1);
    
    %Set up and display an output Rank
    mAb_name=mAb_list(imAb);
    Epitope_sequence=epitope_cell;
    %     Rank=epi_rand_slope_rank;
    Rank=epi_rand_maxmean_rank;
    Trand=table(mAb_name,Epitope_sequence, Rank)
    
    %Display top ten random binding sequences
    [sort_seq,ix]=sort(mean_Pred_Rand(:,4),'descend'); %use highest concentration
    Top_sequences=Randseq_no_shift;
    Top_sequences=(Top_sequences(ix(1:10)));
    for i=1:10
        a=Top_sequences{i};
        a=a(a>'(');
        a=a(a~='X');
        Top_sequences(i)=cellstr(a);
    end
    Tseq=table(Top_sequences)
    
    
    

%*************************************************************************
%*************************************************************************
%                  FUNCTIONS USED
%*************************************************************************


function [m,b]=find_slopes(x,y)

%find the slope and intercept of a line fit through values from each row of
%the x and y variables

xmean=mean(x,2);
x_xmean=x-repmat(xmean,[1,4]);
ymean=mean(y,2);
y_ymean=y-repmat(ymean,[1,4]);
m=sum(x_xmean.*y_ymean,2)./sum(x_xmean.^2,2);
b=ymean-(m.*xmean);
end 

%*************************************************************************

function Composition_binding=Fit2composition(Data,Sequence,aminos)

%This determines the composition and performs a fit
N=length(Sequence);
plen=zeros(N,1);
for n=1:N %determine the length of each sequence
    a=Sequence{n};
    plen(n)=length(a);
end
maxseqlen=max(plen);

%set up a matrix to put the sequences in padded with X
seq=repmat('X',[N,maxseqlen]);
for n=1:N
    a=Sequence{n};
    seq(n,1:length(a))=a;
end

%This determines the number of each amino acid in each sequence
for m=1:length(aminos)
    X(:,m)=sum(seq==aminos(m),2);
end

%now do a linear regression with this
%Lets normalize these so we can compare them 
Xmean=mean(X); %this takes the mean of each column
Xmean_all=repmat(Xmean,[N,1]);
Xnorm=X./Xmean_all;%normalizes prior to the linear fit
Xnorm=[Xnorm,ones(N,1)]; %adds a bias to the fit
Y=Data;
XX=Xnorm'*Xnorm;
XY=Xnorm'*Y;
A=XX\XY;%A is the set of coefficients in the linear fit
Composition_binding=Xnorm*A;
end

%*************************************************************************

function YPred=NN_Predict(NNmodel,Seq,num_residues,numtargets,miniBatchSize,aminos)
%This performs a prediction of the binding given a sequence set and a model
%to predict with

%Need to know how much to pad or cut the sequences
PadSeq=repmat('X',[length(Seq),num_residues]);
Truncate=0;
for i=1:length(Seq)
    a=Seq{i};
    lena=length(a);
    if lena <= num_residues %looks to make sure the input sequences were not too long
    PadSeq(i,1:lena)=a;
    else %if they are, truncate and flag
        PadSeq(i,:)=a(1:num_residues);
        Truncate=Truncate+1;
    end
end
if Truncate>0
    fprintf('%d of the input sequences had to be truncated to run this model\n',Truncate)
end

cellPadSeq=cellstr(PadSeq);

%It turns out that you want to predict in chunks of about 10k for speed and
%memory.  Note that the seq2onehot conversion is also done this way because
%the memory use when you are predicting very large numbers of sequences
%gets very large otherwise.
chunk=10000;
N=length(Seq);
YPred=zeros(N,numtargets);
num_chunk=ceil(N/chunk);
k=1;
for i=1:num_chunk
    if i<num_chunk
        %first one-hot convert the sequences and then do the prediction of
        %binding
        X=seq2onehot(cellPadSeq(k:k+chunk-1),aminos);%this keeps memory use low particularly when running parallel
        YPred(k:k+chunk-1,:)=predict(NNmodel,X,'MiniBatchSize',miniBatchSize);
    else %have to deal with the last chunk differently
        X=seq2onehot(cellPadSeq(k:end),aminos);
        YPred(k:end,:)=predict(NNmodel,X,'MiniBatchSize',miniBatchSize);
    end
    k=k+chunk;
end
end

%**************************************************************************

function OnehotSeq=seq2onehot(Sequence,varargin)
%this function converts a sequence to a one hot form.  You can include a
%list of amino acids if you want.  This is important if there is a specific
%requirement from a previously created network but the sequence you are
%converting does not have all the amino acids.  Note that in this version I
%pad the sequence with 'X', not with char(0) as previously, and I assign X
%a value
N=length(Sequence);

%Is there a common linker at the C-term?  (sometimes we use sequences with
%a common linker sequence on the C-term connecting to the surface)
k=0;
if N>=10
    for i=1:10
        if iscell(Sequence)
            a=Sequence{i};
        else
            a=Sequence(i,:);
        end
        if length(a)>2
            if strcmp(a(end-2:end),'GSG')
                k=k+1;
            end
        end
    end
end
if k==10
    linker=3;
else
    linker=0;
end

if iscell(Sequence) %this is for if the sequences are in cells
    %Need to create a character array of the sequence
    Seq=repmat('X',[N,20]);%a character array with up to 20 characters per line padded with X
    Plen=zeros(N,1);%this is where we put the lengths of the peptide sequences
    for i=1:N
        a=Sequence{i};%convert each sequence to a string of characters
        Plen(i)=length(a)-linker;
        Seq(i,1:Plen(i))=a(1:end-linker);
    end
    R=max(Plen); %need to know how long the max sequence is
else %this is for if the sequences are a character matrix
    Seq=Sequence;
    R=size(Sequence,2);
end

%Now determine which amino acids are being used for the one-hot conversion
Seq=Seq(:,1:R); %get rid of extra spaces past the longest sequence
aa=Seq(1:end);
aa=aa(int32(aa)>32);%get rid of blanks

if isempty(varargin)
    aminos=unique(aa);%find the unique characters
else
    aminos=varargin{1};
end
m=length(aminos);
numseq=zeros(N,R);

%convert all the amino acids to numerical representation
for i=1:m
    numseq(Seq==aminos(i))=i;
end

v=ones(1,m);
U=diag(v); %This is a matrix with ones on the diagonal and zeros elsewehre

%now assign "1" in the appropriate places
OnehotSeq=zeros(N,R*m);
for n=1:N
    pep2D=zeros(m,R);
    for r=1:R
        if numseq(n,r)>0
            pep2D(:,r)=U(numseq(n,r),:);
        end
    end
    OnehotSeq(n,:)=pep2D(1:end);
end
end

%*************************************************************************

function [ShiftData,ShiftSeq]=Sequence_Data_shift(Data,Sequence,shift)
%this function shifts the sequence either starting at the left and shifing
%in 
%num_residues is the total number of residues in the peptide including the
%blanks which will be X
%shift is the number of different shift positions

[N,S]=size(Data);

%calculate max length of Seq
plen=zeros(N,1);
for n=1:N
    plen(n)=length(Sequence{n});
end
maxlen=max(plen);

%Create character array X padded
seqlen=maxlen+shift;
CharSeq=char(zeros(N,seqlen)+double('X'));
for n=1:N
    CharSeq(n,1:plen(n))=Sequence{n};
end

if shift>0
    ShiftData=zeros(shift+1,N,S);
    ShiftSeq=cell(shift+1,N);
    ShiftSeq(1,:)=cellstr(CharSeq);%no shift
    ShiftData(1,:,:)=Data;%no shift
    %this does the shifting:
    for ishift=1:shift
        ShiftSeq(ishift+1,:)=cellstr(circshift(CharSeq,ishift,2)); %shifted sequences in a second dimension
        ShiftData(ishift+1,:,:)=Data;%makes a copy of the data in dimensions 2 and 3 to go with the shifted sequences
    end
else
    ShiftSeq=cellstr(CharSeq);%no shift
    ShiftData=Data;%no shift
end
fprintf('After adding shifted sequences, there are %d Sequences and Data rows\n',N*(shift+1))


%Note that if you set up ShiftSeq(Shift,N) or ShiftData(Shift,N,S) then
%you can use reshape(ShiftData,[Shift*N,S]) and
%reshape(ShiftSeq,[Shift*N,1]) to put them in the right form for training
%and then use reshape(ShiftData,[Shift,N,S]) and reshape(ShiftSeq,[Shift,N])
%to bring them back to 2D or 3D.

end

%**************************************************************************

function ShiftSeq=Sequence_shift(Sequence,shift)
%this function shifts the sequence either starting at the left and shifing

%num_residues is the total number of residues in the peptide including the
%blanks which will be X
%shift is the number of different shift positions

N=length(Sequence);

%calculate max length of Seq
plen=zeros(N,1);
for n=1:N
    plen(n)=length(Sequence{n});
end
maxlen=max(plen);

%Create character array X padded
seqlen=maxlen+shift;
CharSeq=char(zeros(N,seqlen)+double('X'));
for n=1:N
        CharSeq(n,1:plen(n))=Sequence{n};
end

if shift>0
ShiftSeq=cell(shift+1,N);
ShiftSeq(1,:)=cellstr(CharSeq);%no shift
%this does the shifting:
    for ishift=1:shift
        ShiftSeq(ishift+1,:)=cellstr(circshift(CharSeq,ishift,2));
    end
else
    ShiftSeq=cellstr(CharSeq);%no shift
end
fprintf('After adding shifted sequences, there are %d Sequences\n',N*(shift+1))


%Note that if you set up ShiftSeq(Shift,N) and ShiftData(Shift,N,S) then
%you can use reshape(ShiftData,[Shift*N,S]) and
%reshape(ShiftSeq,[Shift*N,1]) to put them in the right form for training
%and then use reshape(ShiftData,[Shift,N,S]) and reshape(ShiftSeq,[Shift,N])
%to bring them back to 2D or 3D.

end

%**************************************************************************

function [Index,Dist]=Weight_constant_density_extended(Data_in,Conc_level,maxdensity,maxcopy)
%This is a function that weights the high binders by creating more
%replicates of them

%First we divide the binding values into bins
X_inc=0.2;
meddata=median(Data_in(:,Conc_level));
data_c=Data_in(:,Conc_level)-meddata;
maxdata=max(data_c);
num_bins=ceil(maxdata/X_inc); %only work on things > median

[N,~]=size(Data_in);

weight=ones(N,1);
Dist=zeros(num_bins,1);
Dist_orig=zeros(num_bins,1);


for i=1:num_bins

    index=(data_c>(i-1)*X_inc & (data_c<=i*X_inc));
    if sum(index)>0 && sum(index)<maxdensity
        weight(index)=maxdensity/sum(index); %create weighting if number in bin too low
    else
        weight(index)=1;
    end
    if weight(index)>maxcopy %however, the number of any one sequence cannot be greater than maxcopy
        weight(index)=maxcopy;
    end

end

weight(weight>maxcopy)=maxcopy;
weight(weight<1)=1;
k=0;

%now make copies of weighted sequences
R=rand([N,1])-1;%random number from-0.5 to 0.5; we use this to randomly pick sequences to weigth
for n=1:N
    weight(n)=round(weight(n)+R(n));% this will randomly push it up or down such that 2.6 will be 3 60% of the time
    if weight(n)<1
        weight(n)=1;
    end
    for i=1:weight(n)
        k=k+1;
        Index(k)=n;
    end
end
fprintf('\nAfter weighting there are now %d Sequences and Data rows\n',k)




end







