Using gene expression data to improve flux predictions in CobraPy

[Ryan Ames]

Hi all,

I want to use measurements of gene expression to constrain/improve flux predictions in a model. This approach defines some simple rules to map gene expression data to reactions. The authors give reactions a weight with an associated standard deviation. I can calculate the weights easy enough but I'm not sure how I can associate these weights to reactions in cobrapy such that they will be used to constrain flux predictions during optimisation or FBA.

Does anyone know how I might achieve this?

Thanks

Ryan

[Christian Diener]

If you only need to change the objective function of the optimization you can use 

> mod.objective = rdict

where rdict is a dictionary of {reaction_id: coefficient}. If your objective function is not linear anymore that would be more complicated and you would need to set up your own solver. For  instance  if  it is quadratic you would need to use the quadratic solvers as shown in https://cobrapy.readthedocs.io/en/stable/qp.html. As a note, even though there are many methods to integrate continuous gene expression values into metabolic modeling, systematic tests of those usually find no improvement when using those as shown in http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003580, so make sure to validate your results. 

[Ryan Ames]

Thanks for your response. As I'm not that familiar with this type of analysis a bit more information would be useful. What would be my new protocol? Do I then optimise twice like:

Set the objective function to a dict of reaction coefficients (gene expression) 

Optimise

Set the objective function to biomass

Optimise

Perform deletion analysis

Also thanks for the PLoS Comp reference. I have done this analysis without the integration of gene expression and have some experiments lined up for validation.

Ryan

[Christian Diener]

Hi Ryan. The method you mentioned does not optimize biomass, but rather looks for the flux distribution that coincides (correlates) most with the gene expression data. It does not require

a biomass function at all (see equation 3 in the paper you mentioned). Looking at the paper in detail I think you would have to implement you own method since the method

they describe also defines some new variables and an iterative method to solve the non-convexity of the absolute value problem. Unfortunately, there is nothing already implemented in

cobrapy that does that for you. 

An alternative approach would be to not use continuous gene expression values, but rather make a data-specific reconstruction based on your gene expression data. Here you first

use a barcode approach to identify enzyme genes that are present or absent in your data and than reconstruct a new model that is consistent with these calls. This is usually done 

with also enforcing some required metabolic targets such as the ability to grow or to produce essential metabolites (this is often called functional validation and is for instance ignored

in the paper you mentioned). So an example: the gene for PFK is not expressed in your data, so as long as its deletion does not prohibit one of the required metabolic functions, the most

consistent thing would be to remove the corresponding reaction from the model (since it is very likely that the enzyme protein is missing). There are several methods that perform

those reconstructions and they generally have better perfromance than the one using continuous gene expression values, but only few are available in Python. 

I have recently implemented one of the newer ones (http://journals.plos.org/ploscompbiol/article?id=10.1371%2Fjournal.pcbi.1004808) and it can be found at https://github.com/cdiener/corda. 

After you have the reconstruction you can predict fluxes using pFBA or FVA and do deletion analysis.  

Cheers

Christian