diff --git a/N803_collector_2.m b/N803_collector_2.m index 629a510..dcc14e8 100644 --- a/N803_collector_2.m +++ b/N803_collector_2.m @@ -4,9 +4,9 @@ % This script will load the MCMC parameter samples, thin them by selecting % every 100th sample, solve the ODE model, and store the results. -%% ------------------------------------------------------------------------ +%% ======================================================================== % INPUTS -% ------------------------------------------------------------------------ +% ======================================================================== name = 'N803_MCMC' ;% file name for storing output @@ -19,30 +19,24 @@ Scenario{2} = 'Fitting to both (C2)' ;% Simultaneous fitting - Cohort 2 Scenario{3} = 'Fitting to Cohort 1 only' ;% Seperate fitting - Cohort 1 Scenario{4} = 'Fitting to Cohort 2 only' ;% Seperate fitting - Cohort 2 -Scenario{5} = 'Fitting no N to both (C1)' ;% Model w/o non-SIV-spec CD8 T -Scenario{6} = 'Fitting no N to both (C2)' ;% Model w/o non-SIV-spec CD8 T -Scenario{7} = 'Webb18 Validation (naive)' ;% Validation of shared fitting -Scenario{8} = 'Webb18 Validation (SIV)' ;% Validation of shared fitting +Scenario{5} = 'Webb18 Validation (naive)' ;% Validation of shared fitting +Scenario{6} = 'Webb18 Validation (SIV)' ;% Validation of shared fitting % x-values (days) corresponding to outputs X{1} = [ (- 4:0)*7 (.1:.1:6.9) (7:.5:44*7) ] ; X{2} = [ (-22:0)*7 (.1:.1:6.9) (7:.5:10*7) ] ; X{3} = [ (- 4:0)*7 (1:.5:44*7) ] ; X{4} = [ (-22:0)*7 (1:.5:10*7) ] ; -X{5} = [ (- 4:0)*7 (1:.5:44*7) ] ; -X{6} = [ (-22:0)*7 (1:.5:10*7) ] ; -X{7} = -7:.1: 4*7 ; -X{8} = -7:.2:18*7 ; +X{5} = -7:.1: 4*7 ; +X{6} = -7:.2:18*7 ; % x-values (days) for doses given Dose{1} = 7*[0:3 5:8 37:40] ; Dose{2} = 7*[0 2 4] ; Dose{3} = 7*[0:3 5:8 37:40] ; Dose{4} = 7*[0 2 4] ; -Dose{5} = 7*[0:3 5:8 37:40] ; -Dose{6} = 7*[0 2 4] ; -Dose{7} = 0 ; -Dose{8} = 7*[0 1 2 7] ; +Dose{5} = 0 ; +Dose{6} = 7*[0 1 2 7] ; % Descriptions of additional model outputs (see 'N803_model_2.m') Output{01} = 'virions [log fold change]' ; @@ -111,9 +105,9 @@ Output{64} = 'R generation by aN [log wrt C1 pre]' ; Output{65} = 'R gen by aN fraction' ; -%% ------------------------------------------------------------------------ +%% ======================================================================== % BODY -% ------------------------------------------------------------------------ +% ======================================================================== % Percolate results structure N = length(Dose) ; @@ -128,7 +122,7 @@ disp('Collecting...') ; % For each model... -for n = 1:8 +for n = 1:6 disp(Scenario{n}) ; % Declare model function handle and results file name @@ -140,15 +134,11 @@ elseif (n==3) || (n==4) saveFile = ['N803_cohort_' num2str(n-2) '_Results'] ; modelFun = @(P) N803_single(X{n}, Dose{n}, P, []) ; - elseif (n==5) || (n==6) - saveFile = 'N803_shared_noN_Results' ; - Doses = { Dose{5} , Dose{6} } ; - modelFun = @(P) N803_shared_noN(X{n}, Doses, P, {[],[]}, n-4) ; - elseif (n==7) + elseif (n==5) saveFile = 'N803_shared_Results' ; modelFun = @(P) N803_shared(X{n}, {Dose{n},[]}, P, {[],[]}, 1,... 'fullOut',true,'ivDosing',880*[.06 0]/.03,'isNaive',true) ; - elseif (n==8) + elseif (n==6) saveFile = 'N803_shared_Results' ; modelFun = @(P) N803_shared(X{n}, {Dose{n},[]}, P, {[],[]}, 1,... 'fullOut',true,'ivDosing',880*[.06 .06 .06 1 0]/.03) ; @@ -221,9 +211,9 @@ C(n).ParSample = PCOH ;% sample of calculated parameters end -%% ------------------------------------------------------------------------ +%% ======================================================================== % OUTPUTS -% ------------------------------------------------------------------------ +% ======================================================================== save([name '.mat'],'C') ; diff --git a/N803_plotter_2.m b/N803_plotter_2.m new file mode 100644 index 0000000..ed8b2ea --- /dev/null +++ b/N803_plotter_2.m @@ -0,0 +1,291 @@ +% N803_plotter_2.m - plot outputs from 'N803_model_2.m' +addpath Analyzer +addpath Data +load('N803_MCMC','C') + +% This script will generate all figure panels using the model outputs saved +% by 'collector.m' and using custom plotting function 'gridplot.m'. + +%% ======================================================================== +% BODY +% ======================================================================== + +% close all figures and undock new figures +close all +set(0,'DefaultFigureWindowStyle','normal') + +% scenario colors +Color{1} = [ 2 1 0 ]/2.4 ;% 'Fitting to both (C1)' ; +Color{2} = [ 0 1 2 ]/2.4 ;% 'Fitting to both (C2)' ; +Color{3} = [ 1 1 1 ]/5 ;% 'Fitting to Cohort 1 only' ; +Color{4} = [ 1 1 1 ]/5 ;% 'Fitting to Cohort 2 only' ; +Color{5} = [ 1 1 1 ]/5 ;% 'Webb18 Validation (naive)' ; +Color{6} = [ 1 1 1 ]/5 ;% 'Webb18 Validation (SIV)' ; + +% Loop through 'Fig' and create each figure +% NOTE: 'for' indent is removed +for Fig = 2:.5:14.5 + +clear P L + +% Declare axis specifications (see 'gridplot.m') +FigSizes = [0 0] ;% figure [x,y] size (not enforced) +Gaps = [.25 .25] ;% axes [x,y] gaps +isBroken = [1 0] ;% x-axis is broken +Margins = [.875 .5 .1 .1] ;% figure [left,bottom,right,top] margins +AxisProps.FontSize = 10 ; +AxisProps.LineWidth = 1 ; +AxisProps.TickDir = 'out' ; +AxisProps.XTickLabelRotation = 45 ; + +% Initialize some variables used for setting up 'gridplot.m' inputs +DatID = zeros(1,10) ;% used for mapping experimental data to axes +OutID = zeros(1,10) ;% used for mapping model outputs to axes +plotMed = false(1) ;% plot median line +plotCI = true(1) ;% plot 95% CI region ( OR plot whole sample ) +xSpace = 7 ;% converts days to weeks +XTickSpace = .2 ;% inches between X Tick Marks +DataColor = [1 1 1]*.55 ;% color for data points + +% Update figure settings for current figure +% NOTE: 'switch' indent is removed +switch Fig + +case {2,2.5} ; c = 1+2*rem(Fig,1) ;% Fitting (fold) ----------------------- + + S = c ;% scenario + figName = [ 'Fig2 Fitting c' num2str(c) ] ; + + YLengths = [.1 1.5 1.5] ; + P{1,3}.YLabel.String = {'Virus','fold change'} ; + P{1,2}.YLabel.String = {'{CD8+ T cells}','fold change',' '} ; + [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( -3:1 , [], 1 ) ; + [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( 0:6 ) ; + OutID(3) = 01 ; DatID(3) = 01 ; + OutID(2) = 02 ; DatID(2) = 02 ; + + load(['N803_DataSet_' num2str(c)]) + X = {DataSet(1).X_data, DataSet(3).X_data} ; + Y = {DataSet(1).Y_data, DataSet(3).Y_data} ; + [X,Y,Sid,Cid] = prepper(X,Y,{2,[]},{2,[]},{0,0},[1 0]) ; + + if c == 1 + [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -4:2:16 ) ; + [P{2,1}.XTick,P{2,1}.XTickLabel] = tickmaker( 32:2:44 ) ; + else + [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -22:6:-4 ) ; + [P{2,1}.XTick,P{2,1}.XTickLabel] = tickmaker( 0:2:8 ) ; + end + +case {3,3.5} ; c = 1+2*rem(Fig,1) ;% GZB Validation ----------------------- + + S = c ;% scenario + figName = [ 'Fig3 Validation c' num2str(c) ] ; + Margins([1 3]) = .5 ; + AxisProps.XTickLabelRotation = 0 ; +% AxisProps.YAxisLocation = 'right' ; + + YLengths = [.1 1.5] ; + P{1,2}.YTick = 0:.2:1 ; + P{1,2}.YTickLabel = {'0%','20%','40%','60%','80%','100%'} ; + OutID(2) = 27 ; DatID(2) = 01 ; + + xSpace = 1 ; + XTickSpace = .15 ; + [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -2:8 , [-2,-1] ) ; + + load('N803_DataSet_2') + if c == 1 + X = {DataSet(16).X_data(1:4) } ; + Y = {DataSet(16).Y_data(1:4) } ; + else + X = {DataSet(15).X_data } ; + Y = {DataSet(15).Y_data } ; + end + [X,Y,Sid,Cid] = prepper(X,Y,{[],[]},{[],[]},{[],[]},[0 0]) ; + Y{1} = Y{1}/100 ; + +case 4 ; c = 1 ;% Killing Weeds (alternate) ------------------------------- + + S = [1 2] ;% scenario + figName = 'Fig4 Killing Weeds' ; + Margins(1) = .5 ; + + YLengths = [0.2 .75 .75 1 .75] ; + [P{1,5}.YTick,P{1,5}.YTickLabel] = tickmaker( -1:2 , [], 1 ) ; + [P{1,4}.YTick,P{1,4}.YTickLabel] = tickmaker( -1:3 , [], 1 ) ; + [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( -3:0 , [], 1 ) ; + [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( -2:1 , [], 1 ) ; + OutID(5) = 58 ;% S killing rate (gS*) + OutID(4) = 29 ;% Sa + OutID(3) = 60 ;% R multiplier for gS* + OutID(2) = 4 ;% S0 + + [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -1:10 ) ; + XTickSpace = .15 ; + +case 4.5 ; c = 1 ;% Cell Weeds (alternate) -------------------------------- + + S = [1 2] ;% scenario + figName = 'Fig4 Cell Weeds' ; + Margins(1) = .5 ; + + YLengths = [0.2 [3 3 4 5]*3.25/15 ] ; + [P{1,5}.YTick,P{1,5}.YTickLabel] = tickmaker( -1:4 , [], 1 ) ; + [P{1,4}.YTick,P{1,4}.YTickLabel] = tickmaker( -2:2 , [], 1 ) ; + [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( -3:0 , [], 1 ) ; + [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( 0:3 , [], 1 ) ; + OutID(5) = 39 ;% S0 activation + OutID(4) = 53 ;% V multiplier for aS* + OutID(3) = 51 ;% R multiplier for aS* + OutID(2) = 49 ;% C multiplier for aS* + + [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -1:10 ) ; + XTickSpace = .15 ; + +case {13,13.5} ; c = 1+2*rem(Fig,1) ;% Single Fitting --------------------- + + S = [0 2]+c ;% scenario + figName = [ 'FigS3 Single Fitting c' num2str(c) ] ; + + YLengths = [.1 1.5 1.5] ; + P{1,3}.YLabel.String = {'Virus','fold change'} ; + P{1,2}.YLabel.String = {'{CD8+ T cells}','fold change',' '} ; + [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( -3:1 , [], 1 ) ; + [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( 0:6 ) ; + OutID(3) = 01 ; DatID(3) = 01 ; + OutID(2) = 02 ; DatID(2) = 02 ; + + load(['N803_DataSet_' num2str(c)]) + X = {DataSet(1).X_data, DataSet(3).X_data} ; + Y = {DataSet(1).Y_data, DataSet(3).Y_data} ; + [X,Y,Sid,Cid] = prepper(X,Y,{2,[]},{2,[]},{0,0},[1 0]) ; + + if c == 1 + [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -4:2:16 ) ; + [P{2,1}.XTick,P{2,1}.XTickLabel] = tickmaker( 32:2:44 ) ; + else + [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -22:6:-4 ) ; + [P{2,1}.XTick,P{2,1}.XTickLabel] = tickmaker( 0:2:8 ) ; + end + +case {14,14.5} ; c = 1+2*rem(Fig,1) ;% Webb Validation -------------------- + + S = 4+c ;% scenario + + figName = ['FigS4 Webb ' num2str(c) ] ; + + YLengths = [.1 1.5 1.25] ; + P{1,2}.YLabel.String = {'CD8+ T cells','per µL blood'} ; + [P{1,2}.YTick,P{1,2}.YTickLabel] = tickmaker( 0:500:3000 ) ; + OutID(2) = 21 ; DatID(2) = 3 ; + + XTickSpace = .35 ; + [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -1:2 ) ; + + load('Webb_DataSet') + X = {DataSet(1).X_data, DataSet(2).X_data, DataSet(3).X_data} ; + Y = {DataSet(1).Y_data, DataSet(2).Y_data, DataSet(3).Y_data} ; + [X,Y,Sid,Cid] = prepper(X,Y,[],[],[],[]) ; + + if c == 2 + P{1,3}.YLabel.String = {'SIV RNA copies','per mL blood'} ; + [P{1,3}.YTick,P{1,3}.YTickLabel] = tickmaker( 0:5 , [], 1 ) ; + OutID(3) = 3 ; DatID(3) = 1 ; + DatID(2) = 2 ; + XTickSpace = .2 ; + [P{1,1}.XTick,P{1,1}.XTickLabel] = tickmaker( -1:18 ) ; + end + +otherwise + continue +end + +[nXaxes,nYaxes] = size(P) ;% number of [x,y] axes + +% Plot doses +if OutID(1) == 0 + [P{1,1}.YTick,P{1,1}.YTickLabel] = tickmaker( [0 .5] , [0 .5] ) ; + L{1,1}{1}.X = C(S(1)).Dose / xSpace ; + L{1,1}{1}.Y = 0.1*ones(size(C(S(1)).Dose)) ; + L{1,1}{1}.LineStyle = 'none' ; + L{1,1}{1}.Marker = '^' ; + L{1,1}{1}.MarkerSize = 5 ; + L{1,1}{1}.MarkerEdgeColor = 'none' ; + L{1,1}{1}.MarkerFaceColor = 'k' ; + if any(Fig == [4,4.5]) % for Weeds plots + L{1,1}{1}.MarkerFaceColor = Color{S(1)} ; + L{1,1}{2} = L{1,1}{1} ; + L{1,1}{2}.X = C(S(2)).Dose / xSpace ; + L{1,1}{2}.Y = 0.3*ones(size(C(S(2)).Dose)) ; + L{1,1}{2}.MarkerFaceColor = Color{S(2)} ; + end +end + +% Plot model outputs +for y = 1:length(YLengths) + o = OutID(y) ; + if o == 0 ; continue ; end + lid = 0 ; + for s = S + lid = lid + 1 ; + L{1,y}{lid}.X = C(s).X / xSpace ; + L{1,y}{lid}.Color = Color{s} ; + L{1,y}{lid}.Y = [ C(s).YSummary(o).lo95 C(s).YSummary(o).hi95 ] ; + L{1,y}{lid}.LineWidth = .5 ; + L{1,y}{lid}.FillColor = [ Color{s} .4 ] ; + end +end + +% Plot data +dataLine.MarkerEdgeColor = DataColor ; +dataLine.MarkerFaceColor = DataColor ; +dataLine.Color = DataColor ; +dataLine.LineWidth = 1 ; +dataLine.LineStyle = ':' ; +DataMarkers = repmat ( {'o','^','s'} , 1 , 5 ) ; +MarkerSize = repmat ( [ 4 , 4 , 5 ]/2 , 1 , 5 ) ; +for y = 1:length(YLengths) + d = DatID(y) ; + if d == 0 ; continue ; end + if OutID(y) == 0 ; lid = 0 ; end + for r = 1:Sid{d}(end) + subMods.X = X{d}(Sid{d}==r)/xSpace ; + subMods.Y = Y{d}(Sid{d}==r) ; + subMods.CensorID = Cid{d}(Sid{d}==r) ; + subMods.Marker = DataMarkers{r} ; + subMods.MarkerSize = MarkerSize(r) ; + L{1,y}{lid+r} = copyfields(subMods,dataLine) ; + if any(Fig == [14,14.5]) + if d == 2 ; L{1,y}{lid+r}.ErrorBars = ... + [104,205,0,0,130,157,124,0,170,580] ; + elseif d == 3 ; L{1,y}{lid+r}.ErrorBars = ... + [191,183,156,109,69,178,94,59] ; + end + L{1,y}{lid+r}.Marker = 's' ; + L{1,y}{lid+r}.MarkerSize = 5 ; + L{1,y}{lid+r}.CensorID = [] ; + break + end + end +end + +% Set [x,y] axes lengths bases on number of ticks +for x = 1:nXaxes + XLengths(x) = ( length(P{x,1}.XTick) - 1 ) * XTickSpace ; %#ok +end + +% Make and save plots +gridplot(P,L ... + ,'AxisProps', AxisProps ... + ,'figName' , figName ... + ,'FigSizes' , FigSizes ... + ,'Gaps' , Gaps ... + ,'imageRes' , '-r300' ... + ,'isBroken' , isBroken ... + ,'Margins' , Margins ... + ,'XLengths' , XLengths ... + ,'YLengths' , YLengths ... + ) ; + +end \ No newline at end of file