Skip to content

Commit

Permalink
Add files via upload
Browse files Browse the repository at this point in the history
  • Loading branch information
cody2 authored Jan 8, 2023
1 parent 9b04ce8 commit a38b9ab
Show file tree
Hide file tree
Showing 2 changed files with 306 additions and 25 deletions.
40 changes: 15 additions & 25 deletions N803_collector_2.m
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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]' ;
Expand Down Expand Up @@ -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) ;
Expand All @@ -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
Expand All @@ -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) ;
Expand Down Expand Up @@ -221,9 +211,9 @@
C(n).ParSample = PCOH ;% sample of calculated parameters
end

%% ------------------------------------------------------------------------
%% ========================================================================
% OUTPUTS
% ------------------------------------------------------------------------
% ========================================================================

save([name '.mat'],'C') ;

Expand Down
291 changes: 291 additions & 0 deletions N803_plotter_2.m
Original file line number Diff line number Diff line change
@@ -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<SAGROW>
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

0 comments on commit a38b9ab

Please sign in to comment.