forked from BIMIB-DISCo/scFBA
-
Notifications
You must be signed in to change notification settings - Fork 0
/
scFBA.asv
137 lines (103 loc) · 6.03 KB
/
scFBA.asv
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
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
% Script replicate scFBA paper
%% Load HMRcore and dictionary
load('data/DictCORE_ENS2HGNC.mat');
HMRcore = readCbModel('HMRcore.xml');
HMRcore = buildRxnGeneMat(HMRcore); %generate rxnGeneMat field
HMRcore.rev = zeros(length(HMRcore.rxns), 1);
HMRcore.rev(HMRcore.lb < 0 & HMRcore.ub > 0) = 1; %generate rev field for reversible reaction
HMRcore = setScFBABound(HMRcore, 1); %set up boundaries of scFBA paper
%Remove mithocondrial complex I-IV rules
GPR2Remove = {'Complex1','Complex1ROS','Complex2','Complex3','Complex4'};
HMRcore.rules(findIdString(HMRcore.rxns, GPR2Remove)) = {''};
clear GPR2Remove
[~, Ex_id] = EditBoundaries(HMRcore, 'Ex_'); %take index of exchange reactions
IdxExRxns = Ex_id.ID; %cambiare se cambiano la posizione delle ex reaction
[~, Coop_id] = EditBoundaries(HMRcore, '_COOP'); %take index of coop reactions
[~, Biomass_id] = EditBoundaries(HMRcore, 'biomass_synthesis'); %take index of biomass reaction
IdxCoopRxn = [Coop_id.ID; Biomass_id.ID];
ThrSigGen = 0.1; %Threshold for significative gene computing.
%% LUAD
%% Load dataset
load('data/LUADdataset.mat')
% filter out genes not in HMRcore model
LUAD_filt = extractGeneRow(HMRcore, LUADdataset, 1, DictCORE);
% Build scStructures
% H358 cell line
H358 = makeSCdataset(LUAD_filt.H358_Pooled, LUAD_filt{:,7:56}, LUAD_filt.Properties.VariableNames(7:56), LUAD_filt.HGNC_ID, 10^-4);
H358 = Genes_Sign(H358, ThrSigGen, 0, 0);
H358 = RepairNegFalse(H358);
% LCPT45 xenograft from primary tumour
LCPT45 = makeSCdataset(LUAD_filt.LCPT45_Pooled, LUAD_filt{:,63:96}, LUAD_filt.Properties.VariableNames(63:96), LUAD_filt.HGNC_ID, 10^-4);
LCPT45 = Geni_Sign(LCPT45, ThrSigGen, 0, 0);
LCPT45 = RepairNegFalse(LCPT45);
% LCMBT15 zenograft from brain metastasis
LCMBT15 = makeSCdataset(LUAD_filt.LCMBT15_Pooled, LUAD_filt{:,153:202}, LUAD_filt.Properties.VariableNames(154:202), LUAD_filt.HGNC_ID, 10^-4);
LCMBT15 = Geni_Sign(LCMBT15, ThrSigGen, 0, 0);
LCMBT15 = RepairNegFalse(LCMBT15);
%% Performe analisys
changeCobraSolver('gurobi');
%H358
% Compute RAS and integrate them in modelPop
H358 = single2IntPopModel(H358, HMRcore, IdxExRxns, IdxCoopRxn, 's'); % Compute RAS and integrate them in modelPop
fluxIntH358 = optimizeCbModel(H358.modelFVAInt); % optimize model with RAS integrated
% Split single optimize solution in nPop solutions, one for each single
% cells for clustering analisys
[FluxPopH358, RxnPop] = splitScFluxes(H358.modelFVAInt, fluxIntH358, length(H358.CellType));
FluxPopNormH358 = FluxPopH358; % normalize fluxes between 0 and 1.
for i=1:length(RxnPop)
FluxPopNormH358(i,:) = ((FluxPopH358(i,:) - min(FluxPopH358(i,:)))./(max(FluxPopH358(i,:)) - min(FluxPopH358(i,:))));
end
clustergram(FluxPopNormH358( ~isnan(FluxPopNormH358(:,1)) ,:),'Cluster',2,'Colormap',redgreencmap,'symmetric',false)
%LCPT45
% Compute RAS and integrate them in modelPop
LCPT45 = single2IntPopModel(LCPT45, HMRcore, IdxExRxns, IdxCoopRxn, 's'); % Compute RAS and integrate them in modelPop
fluxIntLCPT45 = optimizeCbModel(LCPT45.modelFVAInt); % optimize model with RAS integrated
% Split single optimize solution in nPop solutions, one for each single
% cells for clustering analisys
[FluxPopLCPT45, RxnPop] = splitScFluxes(LCPT45.modelFVAInt, fluxIntLCPT45, length(LCPT45.CellType));
FluxPopNormLCPT45 = FluxPopLCPT45; % normalize fluxes between 0 and 1.
for i=1:length(RxnPop)
FluxPopNormLCPT45(i,:) = ((FluxPopLCPT45(i,:) - min(FluxPopLCPT45(i,:)))./(max(FluxPopLCPT45(i,:)) - min(FluxPopLCPT45(i,:))));
end
clustergram(FluxPopNormLCPT45( ~isnan(FluxPopNormLCPT45(:,1)) ,:),'Cluster',2,'Colormap',redgreencmap,'symmetric',false)
%LCMBT15
% Compute RAS and integrate them in modelPop
LCMBT15 = single2IntPopModel(LCMBT15, HMRcore, IdxExRxns, IdxCoopRxn, 's'); % Compute RAS and integrate them in modelPop
fluxIntLCMBT15 = optimizeCbModel(LCMBT15.modelFVAInt); % optimize model with RAS integrated
% Split single optimize solution in nPop solutions, one for each single
% cells for clustering analisys
[FluxPopLCMBT15, RxnPop] = splitScFluxes(LCMBT15.modelFVAInt, fluxIntLCMBT15, length(LCMBT15.CellType));
FluxPopNormLCMBT15 = FluxPopLCMBT15; % normalize fluxes between 0 and 1.
for i=1:length(RxnPop)
FluxPopNormLCMBT15(i,:) = ((FluxPopLCMBT15(i,:) - min(FluxPopLCMBT15(i,:)))./(max(FluxPopLCMBT15(i,:)) - min(FluxPopLCMBT15(i,:))));
end
clustergram(FluxPopNormLCMBT15( ~isnan(FluxPopNormLCMBT15(:,1)) ,:),'Cluster',2,'Colormap',redgreencmap,'symmetric',false)
%% Breast
%% Load dataset
load('data/BC04dataset.mat')
load('data/BCLN03dataset.mat')
BC04_filt = extractGeneRow(HMRcore, BC04dataset, 1, DictCORE);
BC03LN_filt = extractGeneRow(HMRcore, BC03LNdataset, 1, DictCORE);
% Build scStructures
% BC04 primary breast cancer
BC04 = makeSCdataset(BC04_filt.BC04_Pooled, BC04_filt{:,7:end}, BC04_filt.Properties.VariableNames(7:end), BC04_filt.HGNC_ID, 10^-4);
BC04 = Genes_Sign(BC04, ThrSigGen, 0, 0);
BC04 = RepairNegFalse(BC04);
% BC03LN lymph node metastasis
BC03LN = makeSCdataset(BC03LN_filt.BC03LN_Pooled, BC03LN_filt{:,7:end}, BC03LN_filt.Properties.VariableNames(7:end), BC03LN_filt.HGNC_ID, 10^-4);
BC03LN = Genes_Sign(BC03LN, ThrSigGen, 0, 0);
BC03LN = RepairNegFalse(BC03LN);
%% Performe analisys
changeCobraSolver('gurobi');
%BC04
% Compute RAS and integrate them in modelPop
BC04 = single2IntPopModel(BC04, HMRcore, IdxExRxns, IdxCoopRxn, 's'); % Compute RAS and integrate them in modelPop
fluxIntBC04 = optimizeCbModel(BC04.modelFVAInt); % optimize model with RAS integrated
% Split single optimize solution in nPop solutions, one for each single
% cells for clustering analisys
[FluxPopBC04, RxnPop] = splitScFluxes(BC04.modelFVAInt, fluxIntBC04, length(BC04.CellType));
FluxPopNormBC04 = FluxPopBC04; % normalize fluxes between 0 and 1.
for i=1:length(RxnPop)
FluxPopNormBC04(i,:) = ((FluxPopBC04(i,:) - min(FluxPopBC04(i,:)))./(max(FluxPopBC04(i,:)) - min(FluxPopBC04(i,:))));
end
clustergram(FluxPopNormBC04( ~isnan(FluxPopNormBC04(:,1)) ,:),'Cluster',2,'Colormap',redgreencmap,'symmetric',false)