% NNS primer design
%
% Designs primers for whole-protein saturation mutagenesis. For each amino 
% acid position designs a pair of complementary mutagenesis primers 
% (forward/sense and reverse/antisense) with the following guidelines: the 
% codon corresponding to the amino acid to be mutagenized is replaced by 
% NNS and centered in the primer, the 5’ and 3’ ends should end in C or G, 
% and the melting temperature (Tm) should be approximately 70 °C. Tm is
% determined by the QuikChange formula (Papworth C, Bauer J, Braman J and 
% Wright DA. (1996) Strategies. 9:3-4). 
%
% INPUTS
% 1. The nucleotide sequence of the gene to mutate, plus additional
%    flanking nucleotides to accomodate design of primers for first and 
%    last residue
% 2. number of additional flanking nucleotides (buff_nuc)
% 3. desired minimum Tm
%
% OUTPUTS
% 1. NNSprimer: primer sequences with rows corresponding to amino acid 
%    positions; columns one and two are forward and reverse primers,
%    respectively
% 2. NNSprimer_Tm: Tm for each primer
% 3. NNSprimer_length: three columns containing length of forward primer
%    upstream of NNS, downstream of NNS, and full length, respectively
%
% Michael A. Stiffler  mstiffler@post.harvard.edu
% Victor Salinas       victorhernan.salinas@gmail.com
% 09Sept2015

%% Initialization. Enter sequence information specific to experiment.
% Illustrative default values are for TEM-1 beta-lactamase from pBR322; 

% load TEM-1 sequence having additional nucleotides (e.g. 30) on each side to accomodate primers for first and last residue
seq='tttgcggcattttgccttcctgtttttgctcacccagaaacgctggtgaaagtaaaagatgctgaagatcagttgggtgcacgagtgggttacatcgaactggatctcaacagcggtaagatccttgagagttttcgccccgaagaacgttttccaatgatgagcacttttaaagttctgctatgtggcgcggtattatcccgtgttgacgccgggcaagagcaactcggtcgccgcatacactattctcagaatgacttggttgagtactcaccagtcacagaaaagcatcttacggatggcatgacagtaagagaattatgcagtgctgccataaccatgagtgataacactgcggccaacttacttctgacaacgatcggaggaccgaaggagctaaccgcttttttgcacaacatgggggatcatgtaactcgccttgatcgttgggaaccggagctgaatgaagccataccaaacgacgagcgtgacaccacgatgcctgcagcaatggcaacaacgttgcgcaaactattaactggcgaactacttactctagcttcccggcaacaattaatagactggatggaggcggataaagttgcaggaccacttctgcgctcggcccttccggctggctggtttattgctgataaatctggagccggtgagcgtgggtctcgcggtatcattgcagcactggggccagatggtaagccctcccgtatcgtagttatctacacgacggggagtcaggcaactatggatgaacgaaatagacagatcgctgagataggtgcctcactgattaagcattggtaactgtcagaccaagtttactcatatata';

buff_nuc=30; % set number of additional nucleotides added on each side of sequence
Tm_min=70; % set desired Tm

%% design primers
for i=1:(length(seq)-(2*buff_nuc))/3;

    min_length=5; % set a minimum number of nucleotides primers must have on each side of mutation
    
    Tm=0; % initialize Tm

    while Tm<Tm_min;
        
        % design 5' end
        f_start=buff_nuc-min_length+((i-1)*3);
        f_end=buff_nuc+((i-1)*3);
        f_seq=seq(f_start:f_end);
        f_clamp=0;
        f_extra_nuc=1;
        % 5' end must be G or C, otherwise extend until criteria met
        while f_clamp==0 
            if f_seq(1,1)~='g'&&f_seq(1,1)~='c';
                f_seq=seq((f_start-f_extra_nuc):f_end);
                f_extra_nuc=f_extra_nuc+1;
            else f_clamp=1;
            end
        end
        
        % design 3' end
        r_start=buff_nuc+(i*3)+1;
        r_end=buff_nuc+(i*3)+min_length;
        r_seq=seq(r_start:r_end);
        r_clamp=0;
        r_extra_nuc=1;
        % 3' end must be G or C, otherwise extend until criteria met
        while r_clamp==0
            if r_seq(end)~='g'&&r_seq(end)~='c';
                r_seq=seq(r_start:(r_end+r_extra_nuc));
                r_extra_nuc=r_extra_nuc+1;
            else r_clamp=1;
            end
        end
        
        % determine Tm according to QuikChange formula 
        N=length(f_seq)+length(r_seq)+3; % primer length
        GC=100*(sum(f_seq=='g')+sum(f_seq=='c')+sum(r_seq=='g')+sum(r_seq=='c'))/N; % percent GC
        Tm=81.5+0.41*GC-675/N-3*100/N;
        
        % if Tm is below the desired Tm then repeat design, extending by one nucleotide at both 5' and 3' ends
        if Tm<Tm_min
            min_length=min_length+1;
        end
    end
    
    NNSprimer_Tm(i,1)=Tm; % create variable with Tm
    
    % make full-length sense (forward primer) and antisense (reverse primer) sequences
    sense=f_seq;
    sense(1,(length(f_seq)+1):(length(f_seq)+3))='NNS';
    sense(1,(length(f_seq)+4):(length(f_seq)+3+length(r_seq)))=r_seq;
    antisense=seqrcomplement(sense);
    % create variable containing full-length primer sequences
    NNSprimer(i,1)=cellstr(sense);
    NNSprimer(i,2)=cellstr(antisense);
    
    % create variable containing length of full-length sequence and 5' and 3' to NNS 
    NNSprimer_length(i,1)=length(f_seq);
    NNSprimer_length(i,2)=length(r_seq);
    NNSprimer_length(i,3)=length(sense);
    
    clearvars -except seq buff_nuc Tm_min NNSprimer NNSprimer_Tm NNSprimer_length
end       