-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathPlot_Fig_7.m
99 lines (95 loc) · 4.87 KB
/
Plot_Fig_7.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
%% Pannala et al., (2018)
% Script to plot Figure 7: 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_10h_noMFA.m" and Timbr_10h_MFA.m"
% These files use the timbr.m file
% Output of these files are stored in Table format for each case MFA_10h and NoMFA_10h
% 4) Figure 7 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_10h_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)
changeCobraSolver('glpk'); % or glpk
% % Load COBRA model
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');
rno_irrev = rno_irrev_biomass;
%% list of increased/decreased metabolites from data for 10h exposure based on FDR <0.1 from Supplementary TableS4
In_Data10h=readtable('Mets_increased_10h_data.txt','Delimiter','\t');
De_Data10h=readtable('Mets_decreased_10h_data.txt','Delimiter','\t');
% Load the model results stored in MFA_10h folder for conditions under MFA
FilePath10 = [Fpath,'MFA_10h',filesep,'Table10h'];
load(FilePath10)
V = cell2mat(Table_10h(:,2));
[g,h]=ismember(In_Data10h{:,1},Table_10h(:,1));
if g==1
V_m = V(h);
end
N_Im = V_m(V_m>0.1); % Number of increased metabolites
I_m= [In_Data10h{:,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_Data10h{:,1},Table_10h(:,1));
if g1==1
Vd = V(h1);
end
N_dm = Vd(Vd<-0.1); % Number of decreased metabolites
d_m = [De_Data10h{:,1} num2cell(Vd)]; % list of all metabolites in the decreased direction with model predictions
%% Load the model results stored in NoMFA_10h folder for conditions under NoMFA
nFilePath10 = [Fpath,'NoMFA_10h',filesep,'Table10h'];
load (nFilePath10)
Table_10h_n = Table_10h;
Vn = cell2mat(Table_10h_n(:,2));
[gn,hn]=ismember(In_Data10h{:,1},Table_10h_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_Data10h{:,1} num2cell(Vn_m)]; % list of all metabolites in the increased direction with model predictions
[gn1,hn1]=ismember(De_Data10h{:,1},Table_10h_n(:,1));
if gn1==1
Vdn = Vn(hn1);
end
Nn_dm = Vdn(Vdn<-0.1); % Number of decreased metabolites
dn_m = [De_Data10h{:,1} num2cell(Vdn)]; % list of all metabolites in the decreased direction with model predictions
%%
Mets_changed = [In_Data10h{:,1};De_Data10h{:,1}];
Mets=regexprep(Mets_changed,' \exchange \(fwd)','');
model10h = [cell2mat(I_m(:,2));cell2mat(d_m(:,2))];
In_mets = Mets(1:25);De_mets = Mets(26:end);
Data10h = [In_Data10h{:,2};De_Data10h{:,2}];
In_va = [log2(Data10h(1:25)) cell2mat(In_m(:,2)) cell2mat(I_m(:,2)) ];
De_va = [log2(Data10h(26:end)) cell2mat(dn_m(:,2)) cell2mat(d_m(:,2))];
figure(7); set(figure(7),'Units','inches','Position',[0.5 0.5 5.6 10.4]);
h1 = heatmap({'Data','No MFA','MFA'},[In_mets;De_mets],[round(In_va,2,'significant');round(De_va,2,'significant')],'Colormap',jet,'FontSize',11,'FontName','arial');
h1.CellLabelFormat = '%0.2f';