Skip to content

Commit

Permalink
Merge pull request #39 from sysbiolux/v1.3_Refactor
Browse files Browse the repository at this point in the history
V1.3 refactor
  • Loading branch information
sdelandtsheer authored Sep 8, 2020
2 parents c966d53 + 9680e7d commit ad1af70
Show file tree
Hide file tree
Showing 26 changed files with 1,681 additions and 2,595 deletions.
159 changes: 44 additions & 115 deletions FALCON/DriverFalcon.m
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,9 @@
% |/ |/ \|(_______/(_______/(_______)|/ )_) %
% %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Current version: v1.2 %
% (July 2019) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Current version: v1.3 %
% (July 2020) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CONTEXTUALIZATION OF DBN MODELS OF BIOLOGICAL NETWORKS %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sébastien De Landtsheer: [email protected] %
Expand All @@ -22,161 +22,90 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% cite: De Landtsheer et al., 2017, Bioinformatics, 33, 3431–6 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc, clear all % clear screen and workspace


% ===========================
% Define optimization options
% ===========================
optRound = 4; % Number of optimisation round
MaxFunEvals = 3000; % Number of maximal function being evaluated (3000 = default)
MaxIter = 3000; % Number of maximal iteration being evaluated (3000 = default)
Parallelisation = 1; % Use multiple cores for optimisation? (0=no, 1=yes)
HLbound = 0.5; % Qualitative threshold between high and low inputs
InitIC = 2; % Initialise parameters' distribution (1=uniform, 2=normal, 3=zeros)
ObjFunction = 'weighted'; % Either 'unweighted' or 'weighted'

% Define plotting and saving (0=no, 1=yes)
PlotFitEvolution = 1; % Graph of optimise fitting cost over iteration
PlotFitSummary = 1; % Graph of state values at steady-state versus measurements (all in 1)
PlotFitIndividual = 0; % Graph of state values at steady-state versus measurements (individual)
PlotHeatmapCost = 1; % Heatmaps of optimal costs for each output for each condition absolute cost
PlotStateSummary = 1; % Graph of only state values at steady-sate (all in 1)
PlotStateEvolution = 1; % Graph of state values evolution over the course of the simulation (two graphs)
PlotBiograph = 0; % Graph of network topology, nodes activities, and optimised parameters
PlotAllBiographs = 0; % (Only for machines with strong GPUs) Plot all Biographs above

% Additional analyses after the optimisation with the default setting (0=no, 1=yes)
Resampling_Analysis = 1; % Resampling of experimental data and re-optimise
NDatasets = 50;% Number of artificial datasets from which to resample.

LPSA_Analysis = 1; % Local parameter sensitivity analysis
Fast_Option = 0; % Performing faster LPSA by stopping if fitting costs go over a set threshold value
optRound_LPSA = 5; % Number of optimizations for each perturbed datapoint
LPSA_Increments = 3; % Number of increments for LPSA. Increase for finer resolution

KO_Analysis = 1; % Parameter knock-out analysis
KO_Nodes_Analysis = 1; % Node knock-out analysis
optRound_KO = 5; % Number of optimizations for each KO datapoint.
KO_Nodes_fast = 1; %

KO_Nodes_Analysis_eff = 1; % test different KO efficencies on each node and analyse the entire network based on this information
efficiency_range = 0:0.1:1; % indicate the different KO efficencies you want to test (vector from 0 to 1)

% Model and measurement files
InputFile = 'PDGF_model.xlsx'; %for xls
% InputFile='PDGF_model.sif'; %for sif
MeasFile = 'PDGF_meas.xlsx'; %for xls
% MeasFile={'PDGF_meas_in.csv';'PDGF_meas_out.csv';'PDGF_meas_err.csv'}; %for csv


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Creates a save folder
SaveFolderName = ['Results_' datestr(now)]; FinalFolderName=strrep(SaveFolderName, ':', '.');
mkdir(FinalFolderName) % Automatically generate a folder for saving
% InputFile='PDGF_model.sif'; %for sif
% MeasFile={'PDGF_meas_in.csv';'PDGF_meas_out.csv';'PDGF_meas_err.csv'}; %for csv

% Builds a FALCON model for optimisation
estim = FalconMakeModel(InputFile, MeasFile, HLbound); %make the model
estim = FalconMakeModel(InputFile, MeasFile); %make the model

% Defines optimisation options
estim.options = optimoptions('fmincon', 'TolCon', 1e-6, 'TolFun', 1e-6, 'TolX', 1e-10, 'MaxFunEvals', MaxFunEvals, 'MaxIter', MaxIter); % Default setting
estim.SSthresh = eps*100 ; estim.ObjFunction = ObjFunction;
if InitIC == 1, IC_Dist = 'uniform'; elseif InitIC == 2, IC_Dist = 'normal'; elseif InitIC == 3, IC_Dist = 'scratch'; else, error('Please choose the initial parameter distribution'), end
%%%%
% modify optimization parameters here
%%%%

% Optimization
toc_all = []; x_all = []; fval_all = [];

if Parallelisation
parfor counter = 1:optRound %'parfor' will run the parallel computing toolbox
tic, k = FalconIC(estim, IC_Dist); %initial conditions
[xval, fval] = FalconObjFun(estim, k); %objective function
toc_all = [toc_all; toc]; x_all = [x_all; xval]; fval_all = [fval_all; fval];
end

else
for counter = 1:optRound %'parfor' will run the parallel computing toolbox
tic, k = FalconIC(estim, IC_Dist); %initial conditions
[xval, fval] = FalconObjFun(estim, k); %objective function
toc_all = [toc_all; toc]; x_all = [x_all; xval]; fval_all = [fval_all; fval];
end
close(h);
end

fxt_all = [fval_all x_all toc_all];
beep; pause(0.5); beep;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Retrieving the results

[bestx, meanx, stdx, estim] = FalconResults(estim, fxt_all, estim.param_vector, FinalFolderName);
estim = FalconOptimize(estim);

%%% Re-simulate results based on the best optimised parameter set
[MeanStateValueAll, StdStateValueAll, MeanCostAll, StdCostAll, estim] = FalconSimul(estim, bestx, [PlotFitSummary PlotFitIndividual PlotHeatmapCost PlotStateSummary PlotStateEvolution], FinalFolderName);
[estim, StateValues, MSE] = FalconSimul(estim);

save([FinalFolderName filesep 'OptimizedModel'])
save([estim.FinalFolderName filesep 'OptimizedModel'])

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Analyzing the evolution of fitting cost
if PlotFitEvolution
estim = FalconFitEvol(estim, IC_Dist, FinalFolderName);
if estim.PlotFitEvolution
estim = FalconFitEvol(estim);
end

%%% Displayed optimised network with weights
if PlotBiograph
FalconShowNetwork(estim, PlotAllBiographs, FinalFolderName)
if estim.PlotBiograph
FalconShowNetwork(estim)
end

%%% Obtaining robust estimates for the parameters by resampling
if Resampling_Analysis
optRound_Sampling = optRound;
CV_cutoff = 10; % Percent cut-off for large coefficient of variation
[~,~,estim] = FalconResample(estim, bestx, NDatasets, 3, CV_cutoff, FinalFolderName);
if estim.Resampling_Analysis
estim = FalconResample(estim);
end

%%% Sensitivity analysis
if LPSA_Analysis
if Fast_Option, IsFast = 'fast'; else, IsFast = 'slow'; end
Estimated_Time_LPSA = mean(fxt_all(:, end)) * optRound_LPSA * LPSA_Increments * 3 * length(estim.param_vector);
disp(['Estimated Time for LPSA analysis (fast): ' num2str(Estimated_Time_LPSA) ' seconds']); beep; pause(3); beep;
[~, estim] = FalconLPSA(estim, bestx, MeasFile, HLbound, optRound_LPSA, LPSA_Increments, IsFast, Parallelisation, FinalFolderName);
if estim.LPSA_Analysis
estim = FalconLPSA(estim);
end

%%% Knock-out analysis
if KO_Analysis
Estimated_Time_KO = mean(fxt_all(:, end)) * optRound_KO * length(estim.param_vector);
disp(['Estimated Time for KO analysis: ' num2str(Estimated_Time_KO) ' seconds']); beep; pause(3); beep;
estim = FalconKO(estim, fxt_all, HLbound, optRound_KO, Parallelisation, FinalFolderName);
if estim.KO_Analysis
estim = FalconKO(estim);
end

%%% Nodes Knock-out analysis
if KO_Nodes_Analysis
Estimated_Time_KO = mean(fxt_all(:, end)) * optRound_KO * (length(estim.state_names) - length(estim.Input_idx(1,:)));
disp(['Estimated Time for KO analysis: ' num2str(Estimated_Time_KO) ' seconds']); beep; pause(3); beep;
estim = FalconKONodes(estim, fxt_all, HLbound, optRound_KO, Parallelisation, FinalFolderName);
end

%%% Nodes Knock-out efficency analysis
if KO_Nodes_Analysis_eff
Estimated_Time_KO = mean(fxt_all(:,end)) * numel(efficiency_range) * optRound_KO * (length(estim.state_names)-length(estim.Input_idx(1,:)));
disp(['Estimated Time for KO analysis: ' num2str(Estimated_Time_KO) ' seconds']); beep; pause(3); beep;
estim=FalconKONodes_eff(estim, fxt_all, efficiency_range, HLbound, optRound_KO, FinalFolderName);

if estim.KO_Nodes_Analysis
estim = FalconKONodes(estim);
end

%%% Fast KO
if KO_Nodes_fast
estim = FalconKONodes_fast(estim, fxt_all, HLbound, optRound_KO, Parallelisation, FinalFolderName);
if estim.KO_Nodes_fast
estim = FalconKONodes_fast(estim);
end
%
% %%% Nodes Knock-out efficency analysis
% if estim.KO_Nodes_Analysis_eff
% estim=FalconKONodes_eff(estim, fxt_all, efficiency_range, HLbound, optRound_KO, FinalFolderName);
%
% end

%%% Discriminate fast effects from long-term effects
X = estim.Results.KnockOutNodes.BIC_values;
Y = estim.Results.KnockOutNodesFast.BIC_values;
M1 = min(X,Y); M2 = min(X,Y);
figure, plot(X,Y,'.k', 'MarkerSize', 15), hold on
M1 = min(X,Y); M2 = min(X,Y);
plot([M1, M2], [M1, M2], '-k');
Nodes = estim.state_names;
Nodes(estim.Input_idx(1, :)) = [];
text(X(1), Y(1), 'base')
for tx = 1:length(X)-1
text(X(tx+1), Y(tx+1), Nodes(tx))
end
xlabel('BIC for re-optimized model')
ylabel('BIC for non-adapted model')




% === End of the script === %
Loading

0 comments on commit ad1af70

Please sign in to comment.