a question regarding Shlomi's method implementation within CobraToolbox2

[mg]

Dear all,

I have a question about the implementation of Shlomi's method for
tissue specific flux prediction,
it seems that in cobra2 implementation, after creating the entries for
the equations on RH and RL, a unique run for solving one (global) MILP
problem is done, and one can get a result pretty fast.
On the original paper from Shlomi et al. (2008) there's a step
describing an individual run when setting each reaction (at a time) to
either active or inactive, and then checking its influence on the
overall similarity of the predicted fluxes for all reactions to the
expression data -- meaning sequentially solving the MILP problem for
all reactions when setting a specific one.
Is this the case? and if so, the implementation on cobra2 is doing
somewhat a different thing? Do you think the results are not affected?


We have been using an in-house version of Shlomi's method, based on
the original paper (2008) description, using cplex solver. Our
implementation takes some hours for Recon1 (as well as the original
webservices provided by Shlomi's lab, initially MEAT and then iMAt,
which are not working atm).
Since the version on cobra2 provides a result in some seconds, we were
wondering if the Shlomi's version on the toolbox is actually a 'full'
version of the original..or a 'lighter' version?
Any help on this would be very appreciated!


Thank you very much in advance.

Best regards,
Mafalda.

[Aarash]

Hi Mafalda,

You are correct, this is a 'lighter' version of iMAT as it represents
only one run. A full iMAT would require running multiple times as
described in the original paper.

Aarash

[mg]

Hello Aarash,

thank you for the reply.
So in terms of the results, to what extent can one use them for
additional simulations/analysis?

And just to re-check, the gimme implementation (default) corresponds
to a 'full' one right? which is also available under the TIGER
project.

Thank you very much,
Mafalda.

[Aarash]

The GIMME implementation is the "full" one. It is done similarly to
TIGER where the boolean rules are parsed and present/absence calls are
used.

A quick run of iMAT with one iteration from COBRA 2 gives a solution
similar to an entire full run (sorry I don't have more robust results
to show). However, I caution that one solution is not unique and thus
the original algorithm does many iteration to determine whether a
reaction should be included or not depending on the number of times
the reaction is included versus not.

Aarash

[mg]

Ok, thank you for clarifying.

Best regards,
Mafalda.

[mg]

Dear Aarash,

i'm having again a doubt concerning shlomi's results. It got a bit
confused on my mind.
Shlomi's original method is predicting fluxes based in expression data
and network topology, allowing for the possibility that a reaction
(correspondent gene(s)) is highly expressed but carries zero flux,
e.g. if the surrounding reactions are 'off', and this would the
account for post-transcription downregulation (mRNA high but enzyme
low/inactive or no substrate).

Is this also what one gets when running Shlomi's implementation on
Cobra2? I mean, the output contains ExpressedRxns, UnExpressedRxns,
unknown, UpRegulated, DownRegulated and UnknownIncluded.. and the
description says it's based in mRNA data.. is this then just the
mapping of expression to reactions (through gpr's) and giving results
based on the 'and' & 'or' rulles? and therefore any actual flux
prediction?
or does ExpressedRxns stands for reactions predicted to be active
(flux>0) and UnExpressed to reactions predicted to be inactive
(flux=0)?

Thank you very much in advance.

Best regards,
Mafalda

[Aarash]

The iMat implementation in COBRA 2 is the same as the original Nat
Biotech implementation. The MILP problem has been formulated in the
COBRA syntax. The Rxns output tells you a couple of things. The
Expressed/Unexpressed/Unknown results are based solely off of the
expression data. HOWEVER, the Upregulated, Downregulated, and
UnknownIncluded are determined by the algorithms (either GIMME or
iMat) depending on the one used (post-transcription modification as
the original text calls it). The final model (tissueModel) that is
outputed uses the expression data and the algorithms to build the
model with "post-transcription modification" included.

The only difference between the COBRA 2 implementation of iMAT and the
original is the repeated iterations to determine alternate optimal
solutions to make a call whether a reaction should be included or not.
Similar to FBA which has alternate optimal solutions, iMAT solutions
do as well. The COBRA 2 version can be easily modified with a for loop
around the createTissueSpecificModel module to implement iMAT "full".
This will however take much much longer to run.

Thanks
Aarash

[mg]

Dear Aarash,

thank you very much for your reply.
I am still a bit confused with Shlomi's implementation output.
So following what you said, Upregulated, Downregulates and
UnknownIncuded results are the ones on flux prediction, and they come
from the algorithms specifications (either gimme or shlomi).
Using the data provided on cobra toolbox 2, testTissueModel.mat, i ran
the command "createTissueSpecificModel" with gimme (default) and then
shlomi:

>[testModel,testRxns] = createTissueSpecificModel(model,expressionData)

testRxns =

Expressed: {1769x1 cell}
UnExpressed: {497x1 cell}
unknown: {41x1 cell}
UpRegulated: {52x1 cell}
DownRegulated: {0x1 cell}
UnknownIncluded: {1476x1 cell}

>[testModelS,testRxnsS] = createTissueSpecificModel(model,expressionData,[],[],[],'Shlomi')

testRxnsS =

solution: [1x1 struct]
Expressed: {1769x1 cell}
UnExpressed: {497x1 cell}
unknown: {41x1 cell}
UpRegulated: {497x1 cell}
DownRegulated: {0x1 cell}
UnknownIncluded: {1476x1 cell}


--> what I still don't get is:
1) shlomi puts all UnExpressed reactions to UpRegulated. What does it
mean? Does it mean these 497 rxns are predicted to carry a flux
(active)? or is this just because on shlomi algorithm rxns are not
tentatively removed?
2) no downregulated rxns --> also using own expression data, i had the
same 0 value for downregulated. isn't this a bit weird?


Thank you very much for your help.

Best regards,
Mafalda.

[Aarash]

Hi Mafalda,

The results don't look right. Are you sure, your MILP solver is
properly setup.

Aarash

[mg]

Dear Aarash,

i am checking gurobi's setup. I have an academic license, when i ran
the gurobi test files on the gurobi shell, it worked fine. Then i used
glpkmex inside cobra toolbox2 and when running the test files in cobra
toolbox, also returns the optimal solution.
Could you please paste the correct results from the commands i used
above?

Thank you,
Mafalda.

[mg]

sorry, i meant gurobi mex.
and then using gurobi, the MILP test is passed .

[Aarash Bordbar]

Playing around with the test script in COBRA 2.0, I think there is a bug with the model attached not with the algorithm or the expressionData. When I run iMAT with the testModel in the toolbox I get the same weird results as you. When I run the iMat implementation in COBRA 2.0 using the same commands and expressionData as you with the standard Recon 1 model I typically use for my own research I get very different results.

The results are as follows:

Test Model)     

          Expressed: {1769x1 cell}

        UnExpressed: {498x1 cell}

            unknown: {40x1 cell}

        UpRegulated: {498x1 cell}

      DownRegulated: {0x1 cell}

    UnknownIncluded: {1475x1 cell}

Recon 1 Model) 

         Expressed: {1798x1 cell}

        UnExpressed: {472x1 cell}

            unknown: {37x1 cell}

        UpRegulated: {65x1 cell}

      DownRegulated: {450x1 cell}

    UnknownIncluded: {1328x1 cell}

I checked the MILP solutions and it seems that the test script version is not finding a solution (the MILP solver spits back "stat = 0" which refers to Infeasible MILP). I am not sure why this is happening. I will try to update the test script to use the standard model I am using so you can check if the script is working.

Thanks

Aarash

--
Aarash Bordbar, B.S.
Doctoral Candidate
Systems Biology Research Group
University of California, San Diego
9500 Gilman Dr.
La Jolla, CA 92093-0412
858.822.3181
http://systemsbiology.ucsd.edu/Researchers/Bordbar

[mg]

Hi Aarash,

Thank you for the results.
I have tried to run the createTissuemodel function using own
expression data and Recon1 I downloaded from BiGG.
I had the same king of weird results.
I will try running the commands step by step from the script.

Another thing I was just wondering, Recon1 has on model.c all zeros,
except for one reaction (which i think it's lactate exchange, or
another exchange). Is there a reason why this has to be the case?
is it just that there has to be a '1' on this vector? (i know that
gimme won't work without an objective coeff.)
Because i am using expression data from cells that might not exchange
that much lactate..

Thank you once more,
Mafalda.

> ...
>
> read more »

[Aarash Bordbar]

GIMME requires an objective function to ensure is functional when integrating high-throughput data, hence the 1 in the c vector. If lactate exchange is not a possible objective for the cell of interest, I would change the c vector. This is the main problem with GIMME. For mammalian cells, it is difficult to choose a particular objective of interest.

Aarash