-
Notifications
You must be signed in to change notification settings - Fork 0
/
Plot_Fig_S1.m
109 lines (103 loc) · 5.1 KB
/
Plot_Fig_S1.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
%% Pannala et al., (2018)
% Script to plot Figure S1: model predictions against metabolic data using heatmaps
% iRno model predictions under no MFA and MFA conditions.
% Gene expression change induced metabolite prediciton can be generated by following these steps listed below:
% 1) First, raw RNA-sequencing files were analyzed using the Kallisto/Sleuth platform and the results of
% differential expressed genes (DEGs) for both time points
% "Sleuth_table_liver_5h.csv" and "Sleuth_table_liver_10h.csv" stored in the folder Weights_apap
% 2) Using these DEGs, an R script "Generate_timbr_weights.R", with
% helper file "ncomm_helper.R", reactions weights can be generated
% and results are stored in Weights_apap folder as timbr_weights_rno_liver5h.txt and timbr_weights_rno_liver5h.txt
% 3) The generated reaction weights based on gene expression changes, can
% be integrated now to indentify metabolite alterations. This is done using the m-files for
% two different cases mentioned in the manuscript: "Timbr_5h_noMFA.m" and Timbr_5h_MFA.m"
% These files use the timbr.m file
% Output of these files are stored in Table format for each case MFA_5h and NoMFA_5h
% 4) Figure S1 is generated based on the stored timbr predicitons for each case and compared
% against the metabolite alterations from the global metabolic profiling
% data mapped to the model TableS4 (For simplicity we stored mapped
% metabolites as text file in this folder as
% Mets_increased/decreased_5h_data)
%% Install COBRA toolbox
% intialize the toolbox using the below command
% initCobraToolbox()
% set path to the folder
Spath = mfilename('fullpath');
Loc = regexp(Spath,filesep);
Fpath = Spath(1:Loc(end));
addpath (Fpath)
% Change the solver
changeCobraSolver('glpk'); % or glpk
% % Load COBRA models
load iRno_v2.mat;
% no specific objective functions (biomass or ATP)
obj_value_required = 0;
atp_value_required = 0;
% Set unconstrained reaction bounds to +/- 10^5 because the largest
% consumption rate is fairly close to the default 10^3 reaction bound.
default_bound = 10000; % 10^5
rno_irrev_base = ncomm_blais_model2irrev(rno_cobra);
default_rxns_rno = rno_irrev_base.rxns(rno_irrev_base.ub == 1000);
rno_irrev_default = changeRxnBounds(rno_irrev_base, ...
default_rxns_rno, default_bound, 'u');
rno_irrev_atp = changeRxnBounds(rno_irrev_default,...
'RCR11017_f',atp_value_required,'l');
rno_irrev_biomass = changeRxnBounds(rno_irrev_atp,...
'RCR99999_f',obj_value_required,'l');
% Note: as mentioned above, TIMBR predictions were not performed while
% requiring biomass production.
rno_irrev = rno_irrev_biomass;
%% list of increased metabolites from data for 10h exposure
In_Data5h=readtable('Mets_increased_5h_data.txt','Delimiter','\t');
De_Data5h=readtable('Mets_decreased_5h_data.txt','Delimiter','\t');
% load the model results
FilePath5 = [Fpath,'MFA_5h',filesep,'Table5h'];
load(FilePath5)
V = cell2mat(Table_5h(:,2));
[g,h]=ismember(In_Data5h{:,1},Table_5h(:,1));
if g==1
V_m = V(h);
end
N_Im = V_m(V_m>0.1); % Number of increased metabolites
I_m= [In_Data5h{:,1} num2cell(V_m)]; % list of all metabolites in the increased direction with model predictions
% De_D = De_Data10h{:,1};
% De_Data10h = [De_D(1:19);De_D(21:32)];
[g1,h1]=ismember(De_Data5h{:,1},Table_5h(:,1));
if g1==1
Vd = V(h1);
end
N_dm = Vd(Vd<-0.1); % Number of decreased metabolites
d_m = [De_Data5h{:,1} num2cell(Vd)]; % list of all metabolites in the decreased direction with model predictions
%% get model results for no mfa data from the folder
nFilePath5 = [Fpath,'NoMFA_5h',filesep,'Table5h'];
load (nFilePath5)
Table_5h_n = Table_5h;
Vn = cell2mat(Table_5h_n(:,2));
[gn,hn]=ismember(In_Data5h{:,1},Table_5h_n(:,1));
if gn==1
Vn_m = Vn(hn);
end
Nn_Im = Vn_m(Vn_m>0.1); % Number of increased metabolites
In_m= [In_Data5h{:,1} num2cell(Vn_m)]; % list of all metabolites in the increased direction with model predictions
% De_D = De_Data10h{:,1};
% De_Data10h = [De_D(1:19);De_D(21:32)];
[gn1,hn1]=ismember(De_Data5h{:,1},Table_5h_n(:,1));
if gn1==1
Vdn = Vn(hn1);
end
Nn_dm = Vdn(Vdn<-0.1); % Number of decreased metabolites
dn_m = [De_Data5h{:,1} num2cell(Vdn)]; % list of all metabolites in the decreased direction with model predictions
%%
Mets_changed = [In_Data5h{:,1};De_Data5h{:,1}];
Mets=regexprep(Mets_changed,' \exchange \(fwd)','');
model5h = [cell2mat(I_m(:,2));cell2mat(d_m(:,2))];
In_mets = Mets(1:52);De_mets = Mets(53:end);
Data5h = [In_Data5h{:,2};De_Data5h{:,2}];
In_va = [log2(Data5h(1:52)) cell2mat(In_m(:,2)) cell2mat(I_m(:,2)) ];
De_va = [log2(Data5h(53:end)) cell2mat(dn_m(:,2)) cell2mat(d_m(:,2))];
figure(1); set(figure(1),'Units','inches','Position',[0.5 0.5 5.6 9.4]);
h1 = heatmap({'Data','No MFA','MFA'},In_mets,round(In_va,2),'Colormap',jet,'FontSize',11,'FontName','arial');
h1.CellLabelFormat = '%0.2f';h1.ColorLimits = [-4.0 3];
figure(2);set(figure(2),'Units','inches','Position',[0.5 0.5 5.6 8.4]);
h2 = heatmap({'Data','No MFA','MFA'},De_mets,round(De_va,2),'Colormap',jet,'FontSize',11,'FontName','arial');
h2.CellLabelFormat = '%0.2f'; h2.ColorLimits = [-4.0 3];