First commit

CODE_OF_CONDUCT.md

@@ -0,0 +1,34 @@
+Contributor Covenant Code of Conduct
+====================================
+
+## Our Pledge
+In the interest of fostering an open and welcoming environment, we as contributors and maintainers pledge to making participation in our project and our community a harassment-free experience for everyone, regardless of age, body size, disability, ethnicity, gender identity and expression, level of experience, nationality, personal appearance, race, religion, or sexual identity and orientation.
+
+## Our Standards
+Examples of behavior that contributes to creating a positive environment include:
+* Using welcoming and inclusive language
+* Being respectful of differing viewpoints and experiences
+* Gracefully accepting constructive criticism
+* Focusing on what is best for the community
+* Showing empathy towards other community members
+
+Examples of unacceptable behavior by participants include:
+* The use of sexualized language or imagery and unwelcome sexual attention or advances
+* Trolling, insulting/derogatory comments, and personal or political attacks
+* Public or private harassment
+* Publishing others' private information, such as a physical or electronic address, without explicit permission
+* Other conduct which could reasonably be considered inappropriate in a professional setting
+
+## Our Responsibilities
+Project contributors are responsible for clarifying the standards of acceptable behavior and are expected to take appropriate and fair corrective action in response to any instances of unacceptable behavior.
+Project contributors have the right and responsibility to remove, edit, or reject comments, commits, code, wiki edits, issues, and other contributions that are not aligned to this Code of Conduct, or to ban temporarily or permanently any contributor for other behaviors that they deem inappropriate, threatening, offensive, or harmful.
+
+## Scope
+This Code of Conduct applies both within project spaces and in public spaces when an individual is representing the project or its community. Examples of representing a project or community include using an official project e-mail address, posting via an official social media account, or acting as an appointed representative at an online or offline event. Representation of a project may be further defined and clarified by project contributors.
+
+## Enforcement
+If a contributor engages in harassing behaviour, the project organizers may take any action they deem appropriate, including warning the offender or expelling them from online forums, online project resources, face-to-face meetings, or any other project-related activity or resource.
+If you are being harassed, notice that someone else is being harassed, or have any other concerns, please contact a member of the project team immediately.  Instances of abusive, harassing, or otherwise unacceptable behavior may be reported by contacting the project team.  All complaints will be reviewed and investigated and will result in a response that is deemed necessary and appropriate to the circumstances.  The project team is obligated to maintain confidentiality with regard to the reporter of an incident.  Further details of specific enforcement policies may be posted separately.
+
+## Attribution
+Portions of this Code of Conduct were adapted from Electron's [Contributor Covenant Code of Conduct](https://github.com/electron/electron/blob/master/CODE_OF_CONDUCT.md), which itself was adapted from the [Contributor Covenant](http://contributor-covenant.org/version/1/4), version 1.4.

CONTRIBUTING.md

@@ -0,0 +1,11 @@
+# Guidelines for contributing to this project
+Any constructive contributions – bug reports, pull requests (code or documentation), suggestions for improvements, and more – are welcome.
+
+## Conduct
+Everyone is asked to read and respect the [code of conduct](CODE_OF_CONDUCT.md) before participating in this project.
+
+## Coordinating work
+A quick way to find out what is currently in the near-term plans for this project is to look at the [GitHub issue tracker](https://github.com/mhucka/readmine/issues), but the possibilities are not limited to what you see there – if you have ideas for new features and enhancements, please feel free to write them up as a new issue or contact the developers directly!
+
+## Submitting contributions
+Please feel free to contact the author directly, or even better, jump right in and use the standard GitHub approach of forking the repo and creating a pull request.  When committing code changes and submitting pull requests, please write a clear log message for your commits.

README.md

@@ -0,0 +1,52 @@
+# Inferring metabolic fluxes from time-course metabolomics data using Dynamic Flux Analysis (DFA) 
+Metabolism is highly dynamic, and is tied to essential cellular processes such as cellular differentiation and proliferation. However, it is challenging to measure the activity of every reaction in a cell. Time-course metabolomics is a useful measure for tracing metabolite dynamics, but most approaches that utilize this data set examines a small part of the metabolic network. 
+
+**Dynamic flux analysis (DFA)** addresses this issue by fitting time-course metabolomics with genome-scale metabolic models to infer metabolic activity throughout the entire metabolic network. 
+
+*Advantages of DFA include*:
+  * **Flexibility of data type:** DFA can use both absolute and relative metabolomics values. 
+  * **Handling large and small data:** DFA can use as little as two time points, while other models require several time points to fit statistical or genome-scale model. 
+  * **Ease of parameter tuning:** DFA models are easily tuned by a single parameter **kappa**, which is a weight proportional to the optimization strength of the objective function.
+
+## Installation
+DFA was implemented in [MATLAB](https://www.mathworks.com/products/matlab.html) (recommended version: 2018+), and requires the [Gurobi Mathematical Programming Solver](https://www.gurobi.com/) (recommended version: 2018+). The file `dfa.m` or the livescript `dfa.mlx` in this repository are used to run DFA. 
+
+## Usage
+The function `dfa` has the following syntax:
+```MATLAB
+[dfaModel, solution] = DFA(model, metabolomics, params)
+```
+
+**INPUTS:**
+  * `model`:             A `structure` of the genome-scale metabolic model in COBRA format. 
+  * `metabolomics`:      A `structure` of the metabolomics data. It requires the following `fields`:
+    * `data`:            A numerical array of the time course metabolomics data. Rows correspond to metabolites, columns correspond to individual time points. 
+    * `position`:        A numerical array mapping the metabolite to compartments in the metabolic model. Rows correspond to metabolites, columns correspond to individual compartments.
+  * `params (optional)`: A `structure` containing parameters and hyperparameters. It takes the following `fields`:
+    * `kappa`:           A scalar denoting the weight to penalize the metabolomics data. The default value is 1.
+    * `kappa2`:          A scalar denoting the weights to minimize the sum of fluxes using parsimonious Flux Balance Analysis (pFBA). The default value is 1E-6.
+    * `norm`:            A string denoting which function to use for normalizing the flux activity coefficients. There are three possible arguments. The default argument is `None`:
+      * `None`:          No normalization for the flux activity coefficients
+      * `MAV`:           Normalizing the flux activity coefficients from the maximum absolute value
+      * `Quantile`:      Quantile normalization of the flux activity coefficients
+    * `pfba`:            A Boolean that determines whether to run pFBA on or off.
+
+**OUTPUTS:**
+  * `dfaModel`:          A `structure` of the genome-scale metabolic model in COBRA format fitted with the flux activity coefficients computed from the time-course metabolomics.
+  * `solution`:          A `structure` containing solutions from the Gurobi solver, including the metabolic fluxes.
+
+## Contributing
+Contributions are welcome! Please read the contributions guide to get started. Also feel free to submit bugs, feature requests, and pull requests.
+
+Additionally, you can support development for DFA by citing the original manuscript: 
+
+Chandrasekaran, Sriram, Jin Zhang, Zhen Sun, Li Zhang, Christian A. Ross, Yu-Chung Huang, John M. Asara, Hu Li, George Q. Daley, and James J. Collins. ["Comprehensive mapping of pluripotent stem cell metabolism using dynamic genome-scale network modeling."](https://www.cell.com/cell-reports/fulltext/S2211-1247(17)31027-6) Cell reports 21, no. 10 (2017): 2965-2977.
+
+## Publications using Dynamic Flux Analysis
+**Stem cell metabolism**
+1. Chandrasekaran, Sriram, Jin Zhang, Zhen Sun, Li Zhang, Christian A. Ross, Yu-Chung Huang, John M. Asara, Hu Li, George Q. Daley, and James J. Collins. ["Comprehensive mapping of pluripotent stem cell metabolism using dynamic genome-scale network modeling."](https://www.cell.com/cell-reports/fulltext/S2211-1247(17)31027-6) Cell reports 21, no. 10 (2017): 2965-2977.
+2. Shen, Fangzhou, Camden Cheek, and Sriram Chandrasekaran. ["Dynamic network modeling of stem cell metabolism."](https://link.springer.com/protocol/10.1007%2F978-1-4939-9224-9_14) In Computational stem cell biology, pp. 305-320. Humana, New York, NY, 2019.
+
+**Cancer metabolism**
+1. Nelson, Barbara S., Lin Lin, Daniel M. Kremer, Cristovão M. Sousa, Cecilia Cotta-Ramusino, Amy Myers, Johanna Ramos et al. ["Tissue of origin dictates GOT1 dependence and confers synthetic lethality to radiotherapy."](https://cancerandmetabolism.biomedcentral.com/articles/10.1186/s40170-019-0202-2) Cancer & Metabolism 8, no. 1 (2020): 1-16.
+2. Campit, Scott, and Sriram Chandrasekaran. ["Inferring Metabolic Flux from Time-Course Metabolomics."](https://link.springer.com/protocol/10.1007%2F978-1-0716-0159-4_13) In Metabolic Flux Analysis in Eukaryotic Cells, pp. 299-313. Humana, New York, NY, 2020.

src/DFA.m

@@ -0,0 +1,155 @@
+function [dynamicModel, solution] = DFA(model, metabolomics, params)
+%% DFA Dynamic flux analysis for time-course metabolomics data
+% |dfa| computes metabolic fluxes using time-course metabolomics data.
+% 
+% *INPUTS*
+%% 
+% # |model|:                A geneome-scale metabolic model in COBRA format
+% # |file|:                  An excel spreadsheet containing the metabolomics 
+% dataset. 
+% # |sheetname|:       The sheet name in the Excel spreadsheet. The default 
+% sheet is 'Sheet1'
+% # |kappa|:                The relative weight for the consistency of the flux 
+% distribution compared to the weight for maximizing the objective function. The 
+% default value is 1
+% # |kappa2|:              The relative weight for minimizing the sum of absolute 
+% flux values of other reactions. The default value is 1E-3
+%% 
+% *OUTPUTS*
+%% 
+% # |dynamicModel|:  A metabolic model with pseudoreactions that correspond 
+% to the time-course metabolomics data
+% # |solution|:          The growth rate objective value from single gene KO
+
+    % Default parameters
+    if (~exist('params.kappa', 'var')) || (isempty(params.kappa))
+        params.kappa = 1;
+    end
+    if (~exist('params.kappa2', 'var')) || (isempty(params.kappa2))
+        params.kappa2 = 1E-3;
+    end
+    if (~exist('params.norm', 'var')) || (isempty(params.norm))
+        params.norm = 'None';
+    end
+
+    % Initialize Gurobi model parameters
+    tmpModel               = model;
+    tmpModel.A             = model.S;
+    tmpModel.obj           = model.c;
+    tmpModel.rhs           = model.b;
+    tmpModel.sense         = repmat('=', [size(model.S, 1), 1]);
+    tmpModel.lb            = model.lb;
+    tmpModel.ub            = model.ub;
+    tmpModel.vtype         = repmat('C', [size(model.S, 2), 1]);
+    tmpModel.modelsense    = 'max';
+    tmpModel.original_size = length(tmpModel.rxns);
+
+    % Compute the linear flux activity coefficient
+    slope = zeros(size(metabolomics.data, 1), 1);
+    for m = 1:size(metabolomics.data, 1)
+        weight   = polyfit(1:size(metabolomics.data, 2), metabolomics.data(m, :), 1);
+        slope(m) = weight(1) / weight(2);
+    end
+
+    % Normalize the flux activity coefficients
+    switch params.norm
+        case 'None'
+            slope                      = slope;
+        case 'MAV'
+            slope(:, 1)                = slope(:, 1) ./ max(abs(slope(:, 1)));
+        case 'Quantile'
+            slope(:, 1)                = slope(:, 1) ./ quantilenorm(abs(slope(:, 1)), 'MEDIAN', true);
+            slope(slope(:, 1) > 1, 1)  = log10(slope(slope(:, 1) > 1, 1)) + 1;
+            slope(slope(:, 1) < -1, 1) = -log10(-1 * slope(slope(:, 1) < -1, 1)) - 1;
+    end
+
+    % Create the right hand side and constrain the model using the flux
+    % activity coefficients
+    constrainedModel = tmpModel;
+    epsilon          = zeros(size(slope(:, 1)));
+    for m = 1:size(metabolomics.positions, 1)
+        metabolite = (metabolomics.positions ~= 0);
+        if ~isempty(metabolite)
+            positions                  = metabolomics.positions(m, :);
+            positions(positions == 0)  = [];
+            epsilon(m)                 = slope(m); 
+
+            alpha = size(constrainedModel.A, 2) + 1;
+            beta  = alpha + 1;
+            gamma = size(constrainedModel.A, 1) + 1;
+            constrainedModel.rhs(positions) = epsilon(m);
+            constrainedModel.rhs(isnan(constrainedModel.rhs)) = 0;
+
+            % metabolite constraint = tolerance - alpha + beta
+            constrainedModel.A(positions, alpha) = 1;
+            constrainedModel.A(positions, beta)  = -1;
+            constrainedModel.vtype(alpha)        = 'C';
+            constrainedModel.lb(beta)            = 0;
+            constrainedModel.lb(alpha)           = 0;
+            constrainedModel.ub(beta)            = 1000;
+            constrainedModel.ub(alpha)           = 1000;
+            constrainedModel.obj(beta)           = -1 * params.kappa; 
+            constrainedModel.obj(alpha)          = -1 * params.kappa;
+            constrainedModel.c(alpha)            = -1 * params.kappa;
+            constrainedModel.c(beta)             = -1 * params.kappa;
+            constrainedModel.rxns{alpha}         = 'alpha';
+            constrainedModel.rxns{beta}          = 'beta';
+            constrainedModel.mets{gamma}         = 'pmet';
+            constrainedModel.sense(alpha)        = '=';
+            constrainedModel.sense(beta)         = '=';
+
+        end
+    end
+
+    % Minimize the sum of absolute flux values for all other reactions
+    % using kappa2 (pFBA)
+    for r = 1:length(constrainedModel.rxns)
+        if constrainedModel.c(r) == 0
+
+            % xi + si >= -eps2
+            %      si >= 0
+            % rho(ri + si)
+            pmet                           = size(constrainedModel.A, 1) + 1;
+            prxn                           = size(constrainedModel.A, 2) + 1;
+            constrainedModel.A(pmet, r)    = 1;
+            constrainedModel.A(pmet, prxn) = 1;
+            constrainedModel.rhs(pmet)     = -0;
+            constrainedModel.sense(pmet)   = '>';
+            constrainedModel.lb(prxn)      = 0;
+            constrainedModel.ub(prxn)      = 1000;
+            constrainedModel.obj(prxn)     = -1 * params.kappa2;
+            constrainedModel.c(prxn)       = -1 * params.kappa2;
+            constrainedModel.rxns{prxn}    = 'pseudorxn';
+            % xi - ri <= eps2
+            %      ri >= 0
+            pmet                           = size(constrainedModel.A, 1) + 1;
+            prxn                           = size(constrainedModel.A, 2) + 1;
+            constrainedModel.A(pmet, r)    = 1;
+            constrainedModel.A(pmet, prxn) = -1;
+            constrainedModel.rhs(pmet)     = 0;
+            constrainedModel.sense(pmet)   = '<';
+            constrainedModel.lb(prxn)      = 0;
+            constrainedModel.ub(prxn)      = 1000;
+            constrainedModel.obj(prxn)     = -1 * params.kappa2;
+            constrainedModel.rxns{prxn}    = 'pseudorxn';
+        end
+    end
+
+    constrainedModel.vtype     = repmat('C', size(constrainedModel.A, 2), 1);
+    dynamicModel               = constrainedModel;
+
+    % Set params for Gurobi
+    params2.outputflag          = 0;
+    params2.Threads             = 2;
+    params2.Seed                = 314;
+    params2.NumericFocus        = 3;
+    solution                    = gurobi(dynamicModel, params2);
+
+    try
+        solution.v              = solution.x(1:length(model.rxns));
+        solution.solverObj      = solution.objval;
+        solution.rhs            = slope;
+    catch ME
+        warning('Infeasible solution obtained for a metabolic model')
+    end
+end