Issues using rMTA function, needed tutorial

[Saiteja-Reddy]

Hi,

I am trying to use the newly added rMTA function in the toolbox. I have a few queries. I am using iMAT from the toolbox to get reconstructed models D (for disease state) and C (for control state) from original model M.

My intended perturbation is from D (Source) to C (Target), in the rMTA analysis.

1. For the model parameter in rMTA function, do I input D or M? As per the documentation, Vref for rMTA is of the source state, so I'm guessing the input is D. 2. for Vref in rMTA parameters, I would be using ACHRSampler to get the reference flux distribution of model D. Is is alright to do this? 3. for getting the rxnFBS parameter, I'm using the diffexprs2rxnFBS function which takes in a model and Vref as an input. Do I give model D and Vref of model D as an input?

It's unclear in the Docs, if I input the original model M or the reconstructed model D for the rMTA and diffexprs2rxnFBS functions.

I understand rMTA came out recently in the toolbox, but it has been really difficult for me to use it without any proper tutorial. Please include the above details in the Docs and/or do add a tutorial for this in COBRA Tutorials.

Thanks.

[lvalcarcel]

Hi @Saiteja-Reddy ,

First of all, thank you for the interest in our tool.

I will try to answer all your questions:

[TSscore,deletedGenes,Vout] = MTA(model, rxnFBS, Vref, alpha, epsilon, varargin)

1. model: is the raw model (M). Reconstructions doesn't make any sense, as once you produce a perturbation, any metabolic model can change its transcriptome and/or epigenome to adapt to that change.
2. rxnFBS: the output from diffexprs2rxnFBS
3. Vref: the mean of n fluxes obtained using ACHRSampler, n usually being 2,000. ACHR is a good way to do it, as the solver works in any OS. You can also try CHCRSampler.
4. alpha: the alpha parameter in the objective function (from 0 to 1)
5. epsilon: the minimum change: unique or by reaction

rxnFBS = diffexprs2rxnFBS(model, diff_exprs, Vref, ‘SeparateTranscript’, ‘’, ‘logFC’, 0, ‘pval’,0.05)

1. model: again, model M. From the model we are only interested in the grRules.
2. diff_exprs: matlab table (as in DOCs https://opencobra.github.io/cobratoolbox/stable/modules/analysis/rMTA/index.html#src.analysis.rMTA.diffexprs2rxnFBS )
3. Vref: same as before.

For both MTA and rMTA we only use the reconstruction D to calculate the reference flux Vref.

I will try to upload a tutorial as soon as possible with all the help and questions that can arise from the use of this tool.

[Saiteja-Reddy]

Hi @lvalcarcel ,

Thanks a lot for the clarification!

Based on your reply, Vref has reference flux of reactions present in reconstruction D. The output mean from ACHRSampler will be of size equal to number of reactions in D.

So, I should be putting Vref as zero for the reactions which are present in model M and not in reconstruction D to get Vref with respect to model M (for same dimensions in diffexprs2rxnFBS), while calculating the rxnFBS for all the reactions in model M. Am I right in assuming this?

[lvalcarcel]

Hi @Saiteja-Reddy ,

Yes, and you need to check the correct order of the reactions. Some of the COBRA reconstruction methods eliminate reactions and invert the sense of other reactions, so it is neccesary to rewrite the Vref vector to correspond to the reactions in model M.

I do recommend to use a previous check: norm(modelM.S*Vref) If is equal or nearly equal to 0, we are ensuring steady-state model.b=0. if not you should check the output of the reconstruction.

Finally, it is highly recoomendable to use the same model for MTA, rMTA and diffexprs2rxnFBS

[lvalcarcel]

Hi @Saiteja-Reddyhttps://github.com/Saiteja-Reddy,

There is a functionality in the develop brach, https://github.com/opencobra/cobratoolbox/blob/b273f8454d18d6b53d6f7a1397b779e9b74d4b44/src/analysis/rMTA/rMTA.m#L82

If you are working with genes, you can input a cell array with the gene names. If you are working with reactions, you can input a cell array with the names of the reactions.

[TSscore,deletedGenes,Vout] = rMTA(model, rxnFBS, Vref, alpha, epsilon, 'rxnKO', 1, 'listKO', {‘rxn1’, ‘rxn2’})

This change is waiting to be updated to the master branch.

From: Sai Teja Reddy Moolamalla [mailto:notifications@github.com] Sent: martes, 7 de mayo de 2019 12:18 To: opencobra/cobratoolbox cobratoolbox@noreply.github.com Cc: Valcarcel, Luis lvalcarcel@tecnun.es; Mention mention@noreply.github.com Subject: Re: [opencobra/cobratoolbox] Issues using rMTA function, needed tutorial (#1457)

Hi @lvalcarcelhttps://github.com/lvalcarcel ,

How do I do the knockout for an input set of reactions? I don't want to do the computation for reactions with no corresponding gene associations.

Can I just skip the calculation of the TS Scores for the unnecessary reactions? At this linehttps://github.com/opencobra/cobratoolbox/blob/master/src/analysis/rMTA/MTA.m#L177. Or do I have to make any other changes to the knockout matrix before?

Thanks.

[Saiteja-Reddy]

Hi @lvalcarcel , I am facing some issues in running the MTA script. For few reactions knockouts I'm getting infeasible solution for,

v_res = MTA_MIQP (CplexModelBest, KOrxn, 'numWorkers', numWorkers, 'timelimit', timelimit, 'printLevel', printLevel);

which is not handled in the script, and the program stops. Do I just don't consider/ignore such knockouts or am I missing something?

[lvalcarcel]

Hi @Saiteja-Reddy ,

If the reaction produces an infeasible solution, I highly recommend to ignore such knock-outs.

The formulation of rMTA shouln't produce infeasible solutions. Maybe these knock-outs are porduceing null fluxes and you have some constrains with fluxes forced to be hgreater than zero, which produces an infeasibility.

I will improve the tool to overcome this kind of problems.

[kath000]

Dear @lvalcarcel ,

I am using rMTA to model heart metabolism (and get some transformation targets). I have ran into two problems (that may or not be related). The first issue is in the diffexprs2rxnFBS function. The model that I am using does not contain equal amounts of genes in the diseased and the healthy state, which is leading to an issue. code below starts from line 74 of diffexpres2rxnFBS()

% As we cannot map the reaction to genes, we can discretize thess values as % (+1, 0 -1), to create geneFBS geneFBS_aux = zeros(size(diff_exprs,1),1); geneFBS_aux(pos_up) = +1; geneFBS_aux(pos_down) = -1;

% Generate geneFBS for all genes in the model geneFBS = zeros(length(model.genes),1); [~,idx] = ismember(diff_exprs.gene,model.genes); geneFBS(idx) = geneFBS_aux; % An error appears here, telling us that idx cannot contain empty/zero values

I have managed to fix this in several ways, but those fixes only let the program run, and I am not sure if any of the fixes are actually correct

In order to fix this code I had to do two things: add an index to geneFBS_aux as (idx), and make idx not contain zeros. I do not understand why geneFBS_aux would also need an index, but the removal of zeros is clear to me: when we take [tilde,idx], (github uses ~ for formatting so I have written out 'tilde') idx becomes a list that is 0 when A is not in B, but otherwise the lowest index of B where A is found. This means (if I understand correctly) that it will not be [1,0,0,1] (which would be fine to index) but [0,15,0,1] instead. We can remove these zeros without affecting the indexing it seems (in [1,0,0,1] removing zeros would lose the order of model.genes, but when using indexing numbers this shouldn't matter right: [0,15,0,1] to [15,1]? So I added: idx = idx(idx tilde=0) before geneFBS(idx) = geneFBS_aux where we add (idx) so that it becomes: geneFBS(idx) = geneFBS_aux(idx). Another fix is to index (idx) and let it go through an if statement before applying a change to geneFBS

idx = idx(idx~=0) geneFBS(idx) = geneFBS_aux(idx)

% or:

for i =1:length(idx) if idx(i) ~= 0 geneFBX(idx(i)) = geneFBS_aux(idx(i)) end end

This may be the whole problem and my fixes may not solve it (these fixes do let the code run, but potentially lose logical coherence and lead to errors down the line).

Either way, the other error is in the running of the rMTA (using gurobi as solver). Without having solveCobraMIQP() in my own MATLAB editor, the code runs up to line 132 of MTA_MIQP: v_res = solution.full(v); %where I am then told that v is exceeding the bounds of solution.full

Afterward, I have saved and loaded solution to examine what the problem is; .full is completely empty and solution.stat is set to 'unfeasible' (I don't have the luxury to ignore some reactions (as you recommended above), because the rMTA program never fully finishes due to the error).

I have then made solveCobraMIQP() into a MATLAB script (so it is not taken from cobra toolbox but ran from my MATLAB folder), and copied the code from the cobra toolbox into it. Running this gives a different issue:

Undefined function or variable 'f'. Error in solveCobraMIQP (line 331) solution.obj = f; Error in MTA_MIQP (line 126) solution = solveCobraMIQP(MIQPproblem, ... Error in rMTA (line 206) v_res = MTA_MIQP (CplexModelBest, KOrxn, 'numWorkers', numWorkers, 'timelimit', timelimit, 'printLevel', printLevel);

I have tried to find out why variable f is not set, but there are several points within the code where it is set to have some value, so I don't understand the error nor how I could fix it

I hope you can shed some light on what the problem is or how to solve it, I do not know enough about the gurobi solver nor about some of the inputs (osense, csense) to think of a way to remake variable f.

Have a nice New Year's Eve and thank you in advance. Kind regards,

Katherine

[lvalcarcel]

Hi @kath000 ,

First issue is correct. I need to change the code, but the elimination of zeros from the idx variable is correct, although it is neccesary to eliminate positions of geneFBS_aux as well, like you have done.

It is more efficient this line of code:

% Generate geneFBS for all genes in the model
geneFBS = zeros(length(model.genes),1);
[~,idx] = ismember(diff_exprs.gene,model.genes);
geneFBS_aux = geneFBS_aux(idx!=0);
idx= idx(idx!=0);
geneFBS(idx) = geneFBS_aux; % An error appears here, telling us that idx cannot contain empty/zero values

I will take a look with the integration with GUROBI and the problem with the f variable

[muhikp]

Dear @lvalcarcel

I am using rMTA with IBM CPLEX 12.10 solver on my COBRA model. When I run this command [TSscore,deletedGenes,Vout] = rMTA(model, rxnFBS, Vrefmean), I encounter the following error:

Total (root+branch&cut) = 0.19 sec. (98.25 ticks) Index exceeds matrix dimensions.

Error in MTA_MIQP (line 132) v_res = solution.full(v);

Error in rMTA (line 206) v_res = MTA_MIQP (CplexModelBest, KOrxn, 'numWorkers', numWorkers, 'timelimit', timelimit, 'printLevel', printLevel);

then the rMTA stops and there are no outputs. I don't know how to resolve this and I was hoping you can help me with this issue. Thank you very much.

Best regards,

Moohebat

[lvalcarcel]

Dear @muhikp ,

This problem is generated because it creates an infeasible problem, and the code for the infeasible problems were not correctly defined in my function. Could you use CPLEX? Add before rMTA the following line of code:

changeCobraSolver('ibm_cplex', 'all');

I have uploaded a pull request to solve this issue. Until then, please activate CPLEX for solving rMTA.

https://github.com/opencobra/cobratoolbox/pull/1727

Best regards,

Luis

[muhikp]

Dear @lvalcarcel,

Thank you for your prompt response. As I have mentioned in my previous comment, I am using IBM CPLEX as my solver and still this error occurs and I don't know how to fix it. I would appreciate your help. Thanks.

Best,

Moohebat

[lvalcarcel]

Dear @muhikp,

I have uploaded a new commit in which there is a new parameter to force Cplex as our solver in rMTA. In addition, there was also an issue with the solver stats. I have solved it.

Please, try with the version I upload here: rMTA.zip

Best,

Luis

[rmtfleming]

Hi Luis, your commit has now been merged with the develop version. Regards, Ronan

On Tue, 9 Mar 2021 at 10:21, Luis Vitores Valcárcel García < notifications@github.com> wrote:

Dear @muhikp https://github.com/muhikp,

I have uploaded a new commit in which there is a new parameter to force Cplex as our solver in rMTA. In addition, there was also an issue with the solver stats. I have solved it.

Please, try with the version I upload here: rMTA.zip https://github.com/opencobra/cobratoolbox/files/6107748/rMTA.zip

Best,

Luis

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/opencobra/cobratoolbox/issues/1457#issuecomment-793673322, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAQMEOV6DO5VRZ5VW2OPEXDTCXZCJANCNFSM4HEA4O7A .

--

Mr. Ronan MT Fleming B.V.M.S. Dip. Math. Ph.D.

Assistant Professor, Division of Systems Biomedicine and Pharmacology, Leiden Academic Centre for Drug Research, Faculty of Science, Leiden University. https://www.universiteitleiden.nl/en/staffmembers/ronan-fleming & H2020 Project Coordinator, Systems Medicine of Mitochondrial Parkinson’s Disease, http://sysmedpd.eu & Senior Lecturer, School of Medicine, National University of Ireland, Galway.

Peer-reviewed publications: https://goo.gl/FZPG23 Mobile: +353 852 109 806 Skype: ronan.fleming

(This message is confidential and may contain privileged information. It is intended for the named recipient only. If you receive it in error please notify me and permanently delete the original message and any copies.)

[muhikp]

Dear @lvalcarcel,

Thank you very much for fixing this issue and providing the code, it worked perfectly. It is much appreciated.

Kind regards,

Moohebat

[lvalcarcel]

Dear @Saiteja-Reddy,

There is now a tutorial for the rMTA tool.

Is it possible to close this thread?

Kind regards, Luis V. Valcárcel