
****************************************************************************
*                       QUALITATIVE FLUX CONTROL                           *
*                                (QFC)                                     *
* A Matlab package to calculate the driver reactions of metabolic networks *
*                      Last update: 14 May 2015                            *
****************************************************************************

Georg Basler (1), Zoran Nikoloski (2), Abdelhalim Larhlimi (3), Albert-László Barabási (4-6), Yang-Yu Liu (5,6)

1 Department of Environmental Protection, Estación Experimental del Zaidín CSIC, Granada, Spain
2 Systems Biology and Mathematical Modeling, Max Planck Institute of Molecular Plant Physiology, Potsdam, Germany
3 Laboratoire d’Informatique de Nantes Atlantique, Université de Nantes, France
4 Center for Complex Network Research and Departments of Physics, Computer Science and Biology, Northeastern University,
Boston, Massachusetts, USA
5 Channing Division of Network Medicine, Department of Medicine, Brigham and Women’s Hospital, Harvard Medical School,
Boston, Massachusetts, USA.
6 Center for Cancer Systems Biology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, Massachusetts, USA

********************
* RUNNING THE CODE *
********************

REQUIREMENTS:

To run this code, you need Matlab (version 2014a or newer recommended) and Tomlab (http://tomopt.com/). A demo version of
Tomlab is available at http://tomopt.com/scripts/register.php.

RUNNING THE CODE:

1. Extract the contents of the enclosed zip file to a local directory.

2. Start Matlab and include the directory containing the unzipped files, including its subdirectories, in the Matlab
search path (see http://www.mathworks.com/help/matlab/ref/path.html).

3. The input files of the metabolic networks analyzed in the paper are found in the exctracted subdirectory 'networks'.
To run the analyes on a different network, the following files must be created and named accordingly:
	- <name>.S
	  Stoichiometric matrix as tab-separated table.
	- <name>.rev
	  List of reversible reactions, where each line in the file specifies either '0' for an irreversible reaction,
	  or '1' for a reversible reaction.
	- <name>.Reactions
	  List of reaction names, where each line in the file contains the name of the corresponding reaction.
	- <name>.Metabolites
	  List of metabolite names, where each line in the file contains the name of the corresponding metabolite.
<name> is an arbitrary network name. See the 'networks' directory for examples of input files.

4. To run all analysis on a given network, execute 'runAll(<name>)', where <name> refers to the network name specifying
the input files. This will run the following calculations:
	- calculation of blocked reactions, full, partial, and directional coupling using the F2C2 tool (Larhlimi et al.
	  (2012). BMC Bioinformatics 13, 57)
	- calculation of anti-coupling and inhibitive coupling
	- creation of flux coupling profiles and flux coupling graph
	- sampling of 100 reaction activity patterns
	- calculation of driver reactions
	
5. The default Tomlab solver which will be used is 'cplex'. If you want to use a different solver, specify the name as
second string argument to runAll, e.g. runAll('Ecoli_Core', 'glpk') for using the glpk solver.

EXAMPLE:

The following command runs all analyses on the E. coli central metabolic network 'Ecoli_Core' contained in the extracted
directory 'networks/Ecoli_Core':
	runAll('Ecoli_Core');
A description of the generated output files is given in the next section.


****************
* OUTPUT FILES *
****************

The following output files are generated by executing the 'runAll' command:

1. <name>.blocked
Tab-separated list of blocked reactions, where the j-th entry indicates '1', if reaction j is blocked, and '0', if reaction j
is non-blocked.

2. <name>.fctable
Coupling table generated by the F2C2 tool, where an entry (i,j) specifies the following relation for the reactions i and j:
	0: the reactions are not fully, partially, nor directionally coupled
	1: the reactions are fully coupled
	2: the reactions are partially coupled
	3: reaction i is directionally coupled to reaction j
	4: reaction j is directionally coupled to reaction i

3. <name>.anticoupling
Coupling table including full, partial, directional, and anti-coupling. This is a temporary table created before calculating
inhibitive coupling. The final table including inhibitive coupling is described next.

4. <name>.coupling
Coupling table including full, partial, directional, anti, and inhibitive coupling. An entry (i,j) specifies the following
relation for the reactions i and j:
	0: the reactions are uncoupled
	1: the reactions are fully coupled
	2: the reactions are partially coupled
	3: reaction i is directionally coupled to reaction j
	4: reaction j is directionally coupled to reaction i
	5: the reactions i and j are anti-coupled
	8: the forward direction of i is inhibitively coupled to the forward direction of j
	16: the forward direction of i is inhibitively coupled to the reversed direction of j
	32: the reversed direction of i is inhibitively coupled to the forward direction of j
	64: the reversed direction of i is inhibitively coupled to the reversed direction of j
For other values, the reactions are coupled by inhibitive coupling and another coupling type, which is
specified by the sum of their values. For example, the value '12' indicates that the forward direction of i is inhibitively
coupled to the forward direction of j (as with '8'), and reaction j is directionally coupled to reaction i (as with '4').

5. <name>.profile
Absolute and relative frequencies of coupling types in the network, and the number of essential reactions. The flux coupling
profile is obtained from the frequencies of coupling types divided by the total number of coupled reaction pairs.

6. <name>.profile.mat
Same as <name>.profile, but as Matlab struct (binary file).

7. <name>.fcgraph
Flux coupling graph as adjacency list. Every line contains source reaction, target reaction, and coupling type, separated by
tabs. The coupling type is:
	1: the reactions are fully coupled
	2: the reactions are partially coupled
	3: reaction i is directionally coupled to reaction j
	5: the reactions i and j are anti-coupled
	6: reaction i is is inhibitively coupled to reaction j

8. <name>.split.fcgraph
As <name>.fcgraph, but here reversible reactions are represented separately by their name, followed by [+] for the forward
direction, or [-] for the reversed direction.

9. <name>.png
Image of the plotted flux coupling profile.

10. <name>.samples
Sampled reaction activity patterns, where each row corresponds to a sample, and each column corresponds to a reaction.
An entry (i,j) = 1 indicates that the reaction j is active in sample i, (i,j) = 0 indicates that the reaction is inactive.

11. <name>.drivers
Driver reactions over sampled reaction activity patterns. An entry (i,j) of '1' indicates that reaction j is a driver for controlling
the activity pattern given by sample i, while '0' indicates that reaction j is not a driver for the corresponding activity pattern.

*************************
* COMMENTS, SUGGESTIONS *
*************************

If you have comments on the tool, or if you receive an error message, please send the error output and the file
'ErrorXYZ_<date,time>', which is created with the date and time the error occurred, to: georgbasler@yahoo.de.
