function makecurves(NUC_A, NUC_B, NUC_C)


pep = 'ACTG';
cod_whole = [];
for i = 1:4
    cod_tmp(1) = pep(i);
    for j = 1:4
        cod_tmp(2) = pep(j);
        for k = 1:4
            cod_tmp(3) = pep(k);
            cod_whole = [cod_whole;cod_tmp];
        end
    end
end


cod_TriNuc = ['TAC';'ATC';'GTT';'TTC';'ACT';...
    'GAC';'TCT';'CAT';'GCT';'CAG';'GAA';'ATG';...
    'AAA';'AAC';'CGT';'CCG';'CTG';'GGT';'TGG'];
% 19 codons in the NT-TriNuc library

id_ex = [];
for i = 1:64
    for j = 1:19
        if(strcmp(cod_whole(i,:),cod_TriNuc(j,:)))
            id_ex = [id_ex;i];
            break;
        end
    end
end
cod_noTriNuc = cod_whole;
cod_noTriNuc(id_ex,:) = [];

%%%
pep_rate_B = Calculatecodons(NUC_B,[cod_TriNuc;cod_noTriNuc]);
pep_rate_C = Calculatecodons(NUC_C,[cod_TriNuc;cod_noTriNuc]);
pep_rate_A = Calculatecodons(NUC_A,[cod_TriNuc;cod_noTriNuc]);



pep_rate(1,:) = pep_rate_A./max(pep_rate_A)*100;
pep_rate(2,:) = pep_rate_B./max(pep_rate_B)*100;
pep_rate(3,:) = pep_rate_C./max(pep_rate_C)*100;


%%% sort the fraction values
[~,IN] = sort (pep_rate(1,:), 'descend');

pep_rate(1,:) = pep_rate(1,IN);
pep_rate(2,:) = pep_rate(2,IN);
pep_rate(3,:) = pep_rate(3,IN);
pep_rate(pep_rate == 0) = 0.001;
codons = [cod_TriNuc;cod_noTriNuc];

codons = codons(IN,:);

fraction = pep_rate';


% calculate total fraction of forbidden codons
pep_rate_A = pep_rate_A(IN');
pep_rate_B = pep_rate_B(IN');
pep_rate_C = pep_rate_C(IN');

FC_A =  pep_rate_A./sum(pep_rate_A);
FC_B =  pep_rate_B./sum(pep_rate_B);
FC_C =  pep_rate_C./sum(pep_rate_C);

TFC_A = round(sum(FC_A(20:end,1))*100,2);
TFC_B = round(sum(FC_B(20:end,1))*100,2);
TFC_C = round(sum(FC_C(20:end,1))*100,2);



% plot curves
x= 1:size(fraction,1);

figure

plot(x,fraction(:,1),'LineStyle','-','Color',[0.2235 0.3255 0.6392],'Marker','o');
hold on
plot(x,fraction(:,2),'LineStyle','-','Color',[0.0431 0.5020 0.2510],'Marker','o');
hold on
plot(x,fraction(:,3),'-ro');
hold on
xlim([0,65]);
text(40,50, 'total contribution of these 45 codons ', 'horiz', 'right', 'color', 'k');
text(45,50, [num2str(TFC_A) '%'], 'horiz', 'right', 'color', [0.2235 0.3255 0.6392]);
text(48,50,[num2str(TFC_B) '%'], 'horiz', 'right', 'color', [0.0431 0.5020 0.2510]);
text(51,50,[num2str(TFC_C) '%'], 'horiz', 'right', 'color', 'r');
set(gca,'YScale','log');
cod_cell = mat2cell(codons,ones(size(cod_whole,1),1),size(cod_whole,2));
set(gca,'xTick',1:size(cod_whole,1));
set(gca,'XTickLabel',cod_cell);
set(gca,'yTick', [0.001,0.01,0.1,1,10,100]);
set(gca, 'YTickLabel', {'0','10^{-2}','10^{-1}','1','10','100'});
set(gca,'YScale','log');
ylim([0.001,100]);
legend('perfect match','Reads with A18T or G59T',...
    'Reads with FR-mismatches in library region',...
    'Location','southwest');
xtl=get(gca,'XTickLabel'); 
xt=get(gca,'XTick');                  
xtextp=xt;                       
text(xtextp(1:19),0.001*ones(1,19),xtl(1:19),'FontWeight','Bold','HorizontalAlignment','right','rotation',90,'fontsize',9); 
text(xtextp(20:end),0.001*ones(1,64-19),xtl(20:end),'Color',[0.2,0.2,0.2],'HorizontalAlignment','right','rotation',90,'fontsize',9); 
set(gca,'xticklabel','');
set(gcf,'Units', 'centimeters', 'Position', [2, 2, 46.5, 11.1]);

hold on
plot([19.5,19.5],[0.001,100],'k--');

end

function Numcodons = Calculatecodons(Nuc,cod)

Numcodons = zeros(size(cod,1),1);
for i = 1:size(Nuc,2)
    if(i == 1)
        Nuc_new = Nuc{i};
    else
        Nuc_new = [Nuc_new;Nuc{i}];
    end
end

%%%%%%%% get the random region %%%%%%%%%%%%%%%%%%%%
for i=1:size(Nuc_new,1)
    Nuc_new{i,1} = Nuc_new{i,1}([4:6 10:18]);

end

% look back inside the variable region
Nuc_new = char(Nuc_new);

IXJ = find(Nuc_new =='J');
IXB = find(Nuc_new == 'B');
IXD = find(Nuc_new == 'D');
IXK = find(Nuc_new == 'K');
IXM = find(Nuc_new == 'M');
IXE = find(Nuc_new == 'E');

Nuc_new(IXJ) = 'A';
Nuc_new(IXB) = 'C';
Nuc_new(IXD) = 'G';
Nuc_new(IXK) = 'T';
Nuc_new(IXM) = 'T';
Nuc_new(IXE) = 'C';

Nuc_new = cellstr(Nuc_new);
        
block_size = 10000;
nob = ceil(size(Nuc_new,1)./block_size);
for i = 1:nob
    if i == nob
        Nuc_block = char(Nuc_new((1+i*block_size-block_size):end));
    else
        Nuc_block = char(Nuc_new((1+i*block_size-block_size):i*block_size));
    end
    Nuc_block = [Nuc_block,repmat(' ',size(Nuc_block,1),size(cod,2)-mod(size(Nuc_block,2),size(cod,2)))];
    Nuc_block = reshape(Nuc_block',size(cod,2),[])';
    for j = 1:size(cod,1)
        Numcodons(j) = Numcodons(j) + sum(sum(Nuc_block == repmat(cod(j,:),size(Nuc_block,1),1),2) == size(cod,2));
    end
end
end