SCRIPT_EWC_NHDBN.m

% This is an example Matlab script on how to apply the EWC NH-DBN model 
% from Shafiee Kamalabad and Grzegorczyk (2019)
% to the yeast gene expression data from Cantone et al. (2009)

% Since the yeast data consist of two separate time series, 
% the yeast data require a special data pre-processing
% this is done in the first 6 commands of the function: 'PROC_YEAST.m'

%%% global DATA_ALL;
%%%
%%% [DATA_1] = SHIFT_DATA(data_1);
%%% [DATA_2] = SHIFT_DATA(data_2);
%%%
%%% for i=1:length(DATA_1)
%%%   DATA_ALL{i} = [DATA_1{i},DATA_2{i}];
%%% end


% At the end of this file, we also describe how to proceed 
% when only one single time series has to be analysed.
% Then the function 'PROC_ONLY_ONE.m' has to be called.
% In 'PROC_ONLY_ONE.m' the 6 inital command lines of 'PROC_YEAST.m' 
% are replaced by two command lines.

%%% global DATA_ALL;
%%% [DATA_ALL] = SHIFT_DATA(data);


% We note that there are no other differences between the PROC functions.
% 'PROC_YEAST.m' and 'PROC_ONLY_ONE.m' call the functions:
% INITIALISE.m
% START.m
% UPDATE.m

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

% For the mathematical details we refer to our paper:
% 'Non-homogeneous dynamic Bayesian networks 
% with edge-wise sequentially coupled parameters'
% accepted by Bioinformatics on 27 August 2019 

% Please note that the text contained in this file can be copied&pasted 
% (line-by-line, piece-wise or fully) into the Matlab Command window.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% There are two data matrices available:
% 'data_on_original.mat'
% and 
% 'data_off_original.mat'

% The yeast gene expression data stem from Cantone et al. (2009)
% and below we will pre-process the data, as described in the paper.

% HOW TO INTERPRET THE DATA?

% Both data matrices have 5 rows, corresponding to the five yeast genes:
% row 1 - CBF1
% row 2 - GAL4
% row 3 - SWI5
% row 4 - GAL80
% row 5 - ASH1

% And the columns refer to equidistant time points.
% t=1,...,16 (on)
% t=1,...,21 (off)

% The true network, as reported in Cantone et al. (2009), has 8 edges: 

% CBF1  -> GAL4  (1->2)
% GAL4  -> GAL80 (2->4)
% GAL4  -> SWI5  (2->3)
% SWI5  -> ASH1  (3->5)
% SWI5  -> CBF1  (3->1)
% SWI5  -> GAL80 (3->4)
% GAL80 -> GAL4  (4->2)
% ASH1  -> CBF1  (5->1)

% Load the true network ('TRUE.mat') from the current working directory:
load('TRUE');
TRUE
% TRUE(i,j) = 1 means that there is an edge from node i to node j
% TRUE(i,j) = 0 means that there is no edge from node i to node j

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

% PREPARATION OF THE MCMC-BASED MODEL INFERENCE

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SELECT THE NUMBER OF MCMC ITERATIONS
% NOTE:
% In total: 'steps*step_iterations' 
% MCMC iterations will be performed 
% where 'step_iterations' is the thin-out factor

% E.g. for performing 5000 MCMC iterations
% and keeping only every 5-th MCMC sample (thining out by the factor 5),
% set:
steps = 1000;        
step_iterations = 5; 
% Total number of MCMC iterations: steps*step_iterations = 5000
% And every 5-th sample will be kept.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SELECT THE HYPERPARAMETERS OF THE EWC NH-DBN MODEL (CF. PAPER) %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% WE HAVE DEFINED SOME HYPERPARAMETERS AS GLOBAL VARIABLES, SO THAT
% THEY ARE ACCESSIBLE FROM ALL FUNCTIONS AND DO NOT ALWAYS HAVE 
% TO BE PROVIDED AS INPUT ARGUMENTS FOR THEM.

% THE HYPERPARAMETERS FOR THE TWO LAMBDA PARAMETERS 
% (LAMBDA_U AND LAMBDA_C) ARE GLOBAL
global alpha_snr;
global beta_snr;
global alpha_coup;
global beta_coup;

% SET HYPERPARAMETERS OF LAMBDA_U 
alpha_snr  = 2;
beta_snr   = 0.2;
% SET HYPERPARAMETERS OF LAMBDA_C
alpha_coup = 2;
beta_coup  = 0.2;

% HYPERPARAMATER FOR THE NOISE VARIANCE PARAMETER (NOT GLOBAL)
nue_var      = 0.01;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LOAD THE TWO YEAST DATA FROM THE CURRENT WORKING DIRECTORY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

load('data_on_original');
load('data_off_original');

% In both data sets 
% - the 5 rows refer to the 5 genes, as explained above
% - the columns refer to the temporal observations

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% STANDARDIZE THE YEAST DATA (CF. SUPPLEMENTARY PAPER, PART C)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% REMOVE WASHING PERIOD POINTS (= REMOVE FIRST TIME POINT FROM BOTH)
data_on  = data_on_original(:,2:end);
data_off = data_off_original(:,2:end);

% DETERMINE NUMBER OF NODES 'n' 
% AND DETERMINE THE NUMBERS OF REMAINING TIME POINTS 'm_1' and 'm_2'
[n,m_1] = size(data_on);  % n=5, m_1 = 16-1 = 15
[n,m_2] = size(data_off); % n=5, m_2 = 21-1 = 20

% APPLY ZSCORE STANDARDIZATION ON MERGED DATA
data = zscore([data_on,data_off]')';

% THEN SEPARATE THE DATA AGAIN
data_on  = data(:,1:m_1);
data_off = data(:,(m_1+1):end);

% THIS IS REQUIRED, AS WE WANT TO SHIFT BOTH DATA SETS SEPARATELY.
% AFTER THE SHIFTS THE SHIFTED DATA WILL BE MERGED AGAIN
% THIS HAPPENS IN THE FUNCTION 'PROC_YEAST.m'
% FOR DETAILS PLEASE SEE SUPPLEMENTARY PAPER, PART C

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SPECIFY THE INITIAL DATA SEGMENTATION THROUGH AN ALLOCATION MATRIX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% For all 'n' genes we initialize with H=1 (no changepoints).

MATRIX = ones(n,m_1+m_2-2);
% Note that we will lose two more observations. 

% For both data sets we lose one data point, 
% as we will shift both data sets separately.
% (Each edge interaction is subject to a time lag of 1.)

% HOW TO INTERPRET MATRIX?
% MATRIX(i,j) = k
% MEANS THAT THE 'j'-TH DATA POINT OF TARGET GENE 'i' BELONGS TO SEGMENT 'h'


% E.G. IF WE WANTED TO INITIALISE WITH THE TRUE DATA SEGMENTATION

% MATRIX = ones(n,m_1+m_2-2);
% MATRIX(:,m_1:end) = 2;

% Then for all genes (=all rows)
% - the first (m_1-1) data points belong to segment h=1
% - the last  (m_2-1) data points belong to segment h=2

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

% BEFORE RUNNING AN MCMC SIMULATION 
% SPECIFY SOME TUNING PARAMETERS
% AND INITIALIZE THE TWO LAMBDA PARAMETERS 

% Specify the maximal number of data segments per gene:
H_max = 10;

% Because of the changepoint location prior, neighbouring changepoints 
% must have at least the distance (\tau_{h+1}-tau_{h}= 2). 
% Give this piece of information to the MCMC algorithm: 
k_transition = 2;

% Initial lambda parameters
lambda_snr  = 1; % corresponds to $\lambda_{u}$ in the paper
lambda_coup = 1; % corresponds to $\lambda_{c}$ in the paper

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% START THE MCMC SIMULATION 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Depending on the selected number of MCMC iterations 
% (selected via 'steps' and 'step_iterations', as explained above)
% the MCMC simulation might take some time.

[Run] = PROC_YEAST(data_on, data_off, steps, step_iterations, H_max, k_transition, lambda_snr, lambda_coup, nue_var, MATRIX);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THE OUTPUT FILE 'Run.mat'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Most important are the sampled network structures.
% From the sampled network structures, the marginal edge posterior 
% probabilities (=edge scores) can be computed.

% Number of sampled network structures:
n_samples = length(Run.dag); % n_samples = steps + 1;

% The MCMC sampled networks are
% Run.dag{2},...,Run.dag{n_samples}

% Note that Run.dag{1} only contains the initial network structure.
% Default initialisation: A network without any edges.
    
% Average the sampled network structures to get estimates of the 
% marginal edge posterior probabilities (= edge scores)

% Initialize:
NET = zeros(n,n);

% When ignoring the first 100 network structures 
% to take a burn-in phase into account: 

burn_in = 100; % 'burn_in' must be lower than 'steps'

% Add all networks matrices (sampled after burn-in) up:
for i_sample=(burn_in+1):n_samples
    NET = NET + Run.dag{i_sample};
end

NET = NET/(n_samples-burn_in);

% Consider the matrix of edge scores:
NET

% DAG(i,j) is the score of the edge i->j

% E.g. the score of the edge 1->2 is
NET(1,2)
% In this application, it is the score for the edge CBF1 -> GAL4 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Network reconstruction accuracy (AUC value)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load the true network:
load('TRUE');

% And compute the AUC value:
[AUC] = COMPUTE_AUC(Run,TRUE,burn_in);

% The precision-recall AUC for this MCMC simulations is:
AUC


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Note that the 'Run' file also contains other interesting output 
% Here we give a few example diagnostics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% [1] Make a plot of the logarithmic scores (log_likelihood + log_priors)
% along the MCMC iterations:

figure(1)
clf
plot(Run.Log_Scores)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [2] The 10-th network structure in Run is

j = 10;

Run.dag{j}
% The corresponding data point segmentation (allocation) is:
Run.matrix{j}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [3] Identifying coupled/uncoupled edges:
% To get insight, which edges of the j-th sampled network 
% were coupled/uncoupled, extract the information as follows:

% e.g. if j=10
j = 10;

OUT = [];

for i_node = 1:n
    OUT = [OUT,Run.VECTORS{j}{i_node}];
end

% Remore the first row (as the first row corresponds to the intercept)
OUT_EDGES_ONLY = OUT(2:end,:);
OUT_EDGES_ONLY

% Compare with the j-th DAG: Run.dag{j}

Run.dag{j}


% OUT_EDGES(i,j) = -1 means that there is           no edge from i to j
% OUT_EDGES(i,j) =  0 means that there is an uncoupled edge from i to j
% OUT_EDGES(i,j) = +1 means that there is a coupled    edge from i to j

% For DAG(i,j) = 0 we have: OUT_EDGES(i,j) = -1 (no edge Z_i->Z_j)
% For DAG(i,j) = 1 we have: either OUT_EDGES(i,j) = 1 or OUT_EDGES(i,j) = 0
% There is an edge i->j and it is either coupled (1) or uncoupled (0).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [4] To get the fractions of coupled/uncoupled edges, proceed as follows:

% The edge score for the edge 1 -> 2 is
parent_node = 1;
child_node  = 2;
edge_score = NET(parent_node,child_node);

edge_score

% Count the number of coupled and uncoupled edges

n_coupled   = 0;
n_uncoupled = 0;

for i_sample=(burn_in+1):n_samples
    
    edge_status = Run.VECTORS{i_sample}{child_node}(1+parent_node); 
    % add 1 to parent_node to jump over the first element for the intercept  

    if(edge_status==0) % edge was uncoupled
        n_uncoupled = n_uncoupled + 1; 
    elseif(edge_status==1) % edge was coupled
         n_coupled = n_coupled + 1;
    else % edge was not even present
         % so do nothing.
    end
end

% Compute the fractions:
p_coupled   = n_coupled/(n_samples-burn_in);
p_uncoupled = n_uncoupled/(n_samples-burn_in);

% Note that
p_coupled + p_uncoupled
% is equal to
edge_score

% With which fractions was the edge 'parent_node -> child_node' 
% 'coupled' and 'uncoupled'? 

% Simply compute the ratios:

[p_coupled,p_uncoupled]/edge_score

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [5] Make trace plots of the sampled lambda parameters

% For example, let's make the trace plots for the third target node (=SWI5)

% Set:
node = 3;

% Collect the sampled parameters in:
lambda_u = []; % uncoupled
lambda_c = []; % coupled

% Let's here include the burn-in samples
for i_sample=1:n_samples % without burn_in:  for i_sample=(burn_in+1):n_samples
    lambda_u = [lambda_u,Run.lambda_snr_vec{i_sample}(node)];
    lambda_c = [lambda_c,Run.lambda_coup_vec{i_sample}(node)];
end

% Make trace plots
figure(2)
clf
subplot(2,1,1)
plot(lambda_u)
subplot(2,1,2)
plot(lambda_c)

% END OF EXAMPLE DIAGNOSTICS



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HOW TO USE THE CODE WHEN THERE IS ONLY ONE SINGLE TIME SERIES?
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If there is only one single time series 'data' available, 
% then the function 'PROC_ONLY_ONE.m' 
% must be used instead of 'PROC_YEAST.m'

% For example:
% Generate a random data set with 'n=5' nodes and 'm=21' time points:
data = randn(5,21);

% For demonstrating purposes, let's introduce 
% a very strong uncoupled edge: '1->2'

data(2,2:11)  = (+2) * data(1,1:10)  + 0.01 * randn(1,10);
data(2,12:21) = (-2) * data(1,11:20) + 0.01 * randn(1,10);

% Alternatively, you can here use your own 'data' matrix instead

% E.g. via:
% data = ...
% load('data')

% From now on 'data' is a Matlab matrix.
% Every row corresponds to a network variable (e.g. a gene).
% Every column corresponds to a time point.


[n,m] = size(data);
% There are: 'n' network variables and 'm' time points.

% Initialize the data point segmentation matrix so that there are
% no changepoints (i.e. H=1 for all variables).
MATRIX = ones(n,m-1); 
% because of the time lag, we have to shift the data
% so that we lose one data point (m -> m-1)

% Start an MCMC simulation with 'steps*step_iterations' MCMC iterations
steps = 1000;
step_iterations = 5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% And as before, we again set:
H_max = 10;
k_transition = 2;
lambda_snr = 1;
lambda_coup = 1;
nue_var = 0.01;
%global alpha_snr;
%global beta_snr;
%global alpha_coup;
%global beta_coup;
alpha_snr  = 2;
beta_snr   = 0.2;
alpha_coup = 2;
beta_coup  = 0.2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Start the MCMC simulation"
% The 'data' will be shifted within the function 'PROC_ONLY_ONE.m':
[Run] = PROC_ONLY_ONE(data, steps, step_iterations, H_max, k_transition, lambda_snr, lambda_coup, nue_var, MATRIX);

% As before, compute the edge scores from the 'Run' file:

n_samples = length(Run.dag);

burn_in = 100;

NET = zeros(n,n);

for i_sample=(burn_in+1):n_samples
    NET = NET + Run.dag{i_sample};
end

% Compute the edge scores:
NET = NET/(n_samples-burn_in);

NET

% Has EWC NH-DBN  'found' the true edge (1->2)?
% The edge score is hopefully near 1:
NET(1,2)

% Also the other diagnostics [1-5] can be applied, as described above.
% For example:

% Make a trace plot of logarithmic scores:
figure(3)
clf
plot(Run.Log_Scores)

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

% Has EWC NH-DBN inferred the changepoint?

% Let's have a look at the sampled allocation vectors for target node 2.

% Collect all allocation vectors for target node 2

node = 2;

ALLOCATION = [];

for i_sample=1:n_samples
    ALLOCATION = [ALLOCATION;Run.matrix{i_sample}(node,:)];
end

% Every row of 'ALLOCATION' is a sampled allocation vector


% What are the probabilities of a changepoint a the different locations? 

% Hopefully there is a high probability at location 10 
% Note that data point 10 refers to the original time point 11, 
% because of the time lag (data shift).

% Posterior probabilities per location:
sum(ALLOCATION(:,2:end) - ALLOCATION(:,1:end-1))/n_samples

% Make a plot of those probabilities:
figure(4)
clf
plot(sum(ALLOCATION(:,2:end) - ALLOCATION(:,1:end-1))/n_samples)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% END OF SCRIPT FILE 
% BY MARCO GRZEGORCZYK 
% BERNOULLI INSTITUTE, GRONINGEN UNIVERSITY, NL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%