% [sIPdata tdata] = serp_ORF(data, pdata startstr, endstr, snt)
%
% Plots ORF of interest from the bulk of data.
%
% data     - output of serp_read
% pdata    - combined replicate output of serp_analyzer
% startstr - sequence tag unambiguously identifying the beginning of the ORF 
% endstr   - sequence tag unambiguously identifying the end of the ORF
% snt      - number of nucleotides to smooth over 
%
% MB 30/04/21
% MS 04/05/21 Figure 3 added
% MS 20/01/22 Figure 4 added to script


function [sIPdata tdata] = serp_ORF(data, pdata, startstr, endstr, snt);

% determine strand
sp = 0;ep = 0;sm = 0; em = 0;
start_idx = strfind(char(data(1).IP{1,5}(:))',startstr);
if size (start_idx , 2) == 1
    disp (['Start matches to ' num2str(start_idx) ' on + strand']),
    sp = 1;
elseif size (start_idx , 2) > 1,
    disp (['AMBIGIOUS! - No. of start matches is ' num2str(size (start_idx , 2)) ' on + strand']),
end
% look for start on - strand
start_idxm = strfind(char(data(1).IP{1,5}(:))',serp_revcomp(startstr));
if size (start_idxm , 2) == 1
    disp (['Start matches to ' num2str(start_idxm) ' on - strand']),
    sm = 1;
elseif size (start_idxm , 2) > 1,
    disp (['AMBIGIOUS! - No. of start matches is ' num2str(size (start_idx , 2)) ' on - strand']),
end    

% look for end on + strand
end_idx = strfind(char(data(1).IP{1,5}(:))',endstr);
end_idx = end_idx + size (endstr,2);
if size (end_idx , 2) == 1
    disp (['End matches to ' num2str(end_idx) ' on + strand']),
    ep = 1;
elseif size (end_idx , 2) > 1,
    disp (['AMBIGIOUS! - No. of end matches is ' num2str(size (end_idx , 2)) ' on + strand']),
end

% look for end on - strand
end_idxm = strfind(char(data(1).IP{1,5}(:))',serp_revcomp(endstr));
end_idxm = end_idxm + size (endstr,2);
if size (end_idxm , 2) == 1
    disp (['End matches to ' num2str(end_idxm) ' on - strand']),
    em = 1;
elseif size (end_idxm , 2) > 1,
    disp (['AMBIGIOUS! - No. of end matches is ' num2str(size (end_idxm , 2)) ' on - strand']),
end

% decide on which strand to work
if sp == 1 & ep == 1,
    disp ( [' Working on + strand, Orf length is ' num2str(end_idx-start_idx) ] );
    tdata = pdata.plustotal;
    IPdata = pdata.plusIP;
    ORF_vector = 1:1:end_idx-start_idx+1;
elseif sm == 1 & em == 1,
    start_idx = end_idxm;
    end_idx = start_idxm;
    disp ( [' Working on - strand, Orf length is ' num2str(end_idx-start_idx) ] );
    tdata = pdata.minustotal;
    IPdata = pdata.minusIP;
    ORF_vector = end_idx-start_idx+1:-1:1;
else disp ('WARNING - ORF assignment is AMBIGIOUS !!!');
    return;
end



close all
%figure (1); title ('ORF Combined - IP in blue - Total in red - Division in cyan'); xlabel ('ORF position'); ylabel ('reads'); hold on;
%plot (ORF_vector , IPdata(start_idx:end_idx),'.b' );
%plot (ORF_vector , tdata(start_idx:end_idx),'.r');
%plot (ORF_vector , IPdata(start_idx:end_idx)./tdata(start_idx:end_idx) ,'.c');
%hold off

%figure (2); title ('Replicate 1 - IP in blue - Total in red - Division in cyan'); xlabel ('ORF position'); ylabel ('reads'); hold on;
%plot (ORF_vector , data(1).IP{1,6}(start_idx:end_idx),'.b' );
%plot (ORF_vector , data(1).total{1,6}(start_idx:end_idx),'.r');
%plot (ORF_vector , data(1).IP{1,6}(start_idx:end_idx) ./ data(1).total{1,6}(start_idx:end_idx) ,'.c');
%hold off



figure (1); title ('ORF Smoothed - IP in blue - Total in red - Division in cyan'); xlabel ('ORF position'); ylabel ('reads'); hold on;
sIPdata = smooth (IPdata(start_idx:end_idx), snt);
stdata = smooth (tdata(start_idx:end_idx), snt );
plot (ORF_vector , sIPdata,'.b' );
plot (ORF_vector , stdata,'.r');
plot (ORF_vector , sIPdata./stdata ,'.c');
hold off

figure (2); title ('IP-red - Total-blue');
plot(ORF_vector/3, sIPdata, '-r');hold on
plot(ORF_vector/3, stdata, '-b'); 
hold off


figure (3); title ('Title_Name'); xlabel ('Codon/ AA'); ylabel ('enrichment [IP/total]'); hold on;
sIPdata = smooth (IPdata(start_idx:end_idx), snt);
stdata = smooth (tdata(start_idx:end_idx), snt );
plot (ORF_vector/3 , sIPdata./stdata, '-b', 'Linewidth', 1.5);
yline([2], '--', 'Color', [0.5,0.5,0.5]);
xlim ([]);
ylim ([]);
box on;
set(gca,'XMinorTick','off','YMinorTick','off');
set(gca,'TickDir','out');
hold off

figure (4); title ('Total Translatome'); xlabel ('Codon/ AA'); ylabel ('Normalized Footprints]'); hold on;
stdata = smooth (tdata(start_idx:end_idx), snt );
bar(ORF_vector/3, stdata, 2, 'histc');
xlim ([665 695]);
ylim ([0 10]);
box on;
set(gca,'XMinorTick','off','YMinorTick','off');
set(gca,'TickDir','out');
hold off

