BINGO_example.m

% This file contains an example of the method with some instructions. For
% more information and for more functionalities, please refer to the
% readme_BINGO file.

%Load data
addpath('./BINGO_files/')
load('./BINGO_files/Example_Data.mat')

%% === EXAMPLE 1: Basic use ===

%The data consist of two time series X1 and X2 with dimension 5. The first
%time series consists of 21 measurements at sampling rate dt1 = 0.5. The 
%second time series consists of 14 measurements at measurement times stored
%in the variable 'times'.

clear('data')
%The time series data should be put to field ts:
data.ts={X1,X2};

%Either the sampling rate or measurement times are put to field Tsam:
data.Tsam={dt1,times};

%% Run BINGO

%Initialization
[data,state,parameters]=BINGO_init(data);

% MCMC Burn-in
[~,chain,~,state,stats]=BINGO(data,state,parameters);
disp_stats(' BURN-IN COMPLETE',stats,chain,parameters.its)

% Actual sampling
parameters.its=10000;
[Plink,chain,xstore,state,stats]=BINGO(data,state,parameters);
disp_stats(' SAMPLING COMPLETE',stats,chain,parameters.its)

%The end result:
confidence_matrix=Plink/chain;


%% Collect more samples

%It is possible to continue collecting MCMC samples to improve the accuracy
%of the method.

chain_old=chain;
Plink_old=Plink;

parameters.its=10000;
xstore_old=xstore;
[Plink,chain,xstore,state,stats]=BINGO(data,state,parameters);

xstore=chain_old/(chain+chain_old)*xstore_old+chain/(chain+chain_old)*xstore;
Plink=Plink_old+Plink;
chain=chain_old+chain;
disp_stats(' SAMPLING COMPLETE',stats,chain,parameters.its)

%The end result:
confidence_matrix=Plink/chain;

%% --- Visualization: Link confidence histogram ---

%It is a good idea to plot a histogram of the confidence matrix to help
%decide on a threshold. Here the ground truth is known, and the links
%corresponding to true links are highlighted in the histogram.

%The ground truth adjacency matrix of the system that generated the data
GroundTruth=[0 0 0 1 1; 1 0 0 0 0; 0 1 0 0 1; 0 0 1 0 0; 0 1 0 0 0];

%Exclude diagonal values
conf_to_plot=confidence_matrix(1:5,1:5)-2*eye(5);

figure; hold on; grid on;
histogram(conf_to_plot,0:.1:1,'FaceColor',[.7 .7 .7])
histogram(conf_to_plot(find(GroundTruth)),0:.1:1,'FaceColor',[.4 .4 .4])


%% --- Visualization: The trajectory estimate ---

%The posterior mean of the continuous expression trajectory is stored in
%the 'xstore' variable, which contains all trajectories concatenated. The
%variable data.plot_index{j} contains the indices of the j^th time series,
%and the variable data.fine_times{j} contains the corresponding time points
%where this trajectory is computed. These are generated by BINGO_init.

%NOTE: The example data is from a white noise driven linear system, and 
%therefore the trajectory estimate mainly resembles a piecewise linear 
%function instead of a smooth curve.

figure('Position',[60 360 1290 420])
for jx=1:5
    subplot(2,5,jx); grid on; hold on;
    plot(data.fine_times{1},xstore(jx,data.plot_index{1}),'LineWidth',1)
    plot(data.Tsam{1},data.ts{1}(jx,:),'ok','MarkerFaceColor','k','MarkerSize',3)
    axis([0 10 0 inf])
    title(['Time series 1, variable ' num2str(jx)'])
    
    subplot(2,5,5+jx); grid on; hold on;
    plot(data.fine_times{2},xstore(jx,data.plot_index{2}),'LineWidth',1)
    plot(data.Tsam{2},data.ts{2}(jx,:),'ok','MarkerFaceColor','k','MarkerSize',3)
    axis([0 10 -inf inf])    
    title(['Time series 2, variable ' num2str(jx)'])
end




%% === EXAMPLE 2: perturbation data series ===

%In the third time series X3, a static perturbation has been applied on the
%system. This can be taken into account in the input-field. The confidence 
%matrix is now 5-by-6 where the sixth column corresponds to the 
%perturbation targets. The sampling rate is the same as in time series 1.

%Form the data structure
clear('data')
data.ts={X1,X3};
data.Tsam={dt1};
data.input={zeros(1,size(X1,2)),ones(1,size(X3,2))};

%Initialization
[data,state,parameters]=BINGO_init(data);

% MCMC Burn-in
[~,chain,~,state,stats]=BINGO(data,state,parameters);
disp_stats(' BURN-IN COMPLETE',stats,chain,parameters.its)

% Actual sampling
parameters.its=20000;
[Plink,chain,xstore,state,stats]=BINGO(data,state,parameters);
disp_stats(' SAMPLING COMPLETE',stats,chain,parameters.its)

%The end result:
confidence_matrix=Plink/chain;



%% === EXAMPLE 3: missing measurements ===

Xmiss=X2;

%Throw away a couple of time points
Xmiss(1,12)=0;
Xmiss(3,4)=0;

clear('data')
data.ts={X1,Xmiss};
data.Tsam={dt1,times};

%Missing measurements:
data.missing={[],[1 12; 3 4]};

%NOTE: the initialization file replaces the missing values by interpolated
%values. The whole procedure of handling missing measurements is described
%in the readme-file.


%Then run normally
[data,state,parameters]=BINGO_init(data);
[~,chain,~,state,stats]=BINGO(data,state,parameters);
disp_stats(' BURN-IN COMPLETE',stats,chain,parameters.its)
parameters.its=20000;
[Plink,chain,xstore,state,stats]=BINGO(data,state,parameters);
disp_stats(' SAMPLING COMPLETE',stats,chain,parameters.its)
confidence_matrix=Plink/chain;



%% === EXAMPLE 4: knockout data, prior information ===

%Time series X4 (sampling times in variable 'times') is formed by
%artificially setting the fourth variable to zero during simulation,
%corresponding to a gene knockout experiment.

clear('data')
data.ts={X1,X4};
data.Tsam={dt1,times};

%Give the indices of the knocked-out genes in each time series
data.ko={[],[4]};

%If it is known a priori that gene 1 is regulated by gene 5 and that gene 2
%is certainly not regulated by gene 3, it can be given like this:
%data.sure=zeros(5,5);
%data.sure(1,5)=1;
%data.sure(2,3)=-1;


%Then run normally
[data,state,parameters]=BINGO_init(data);
[~,chain,~,state,stats]=BINGO(data,state,parameters);
disp_stats(' BURN-IN COMPLETE',stats,chain,parameters.its)
parameters.its=20000;
[Plink,chain,xstore,state,stats]=BINGO(data,state,parameters);
disp_stats(' SAMPLING COMPLETE',stats,chain,parameters.its)
confidence_matrix=Plink/chain;