prepINITModel not ignoring custom reactions in model pruning

[wshao1]

Hi! I am trying to apply ftINIT to only a subset of the Human1 model, by adding the reactions I want to ignore to the customRxnsToIgnore argument of prepINITModel. However, after running ftINIT with and without my customRxnsToIgnore, I see no difference in reaction inclusion. What might be causing this? Am I misinterpreting how to use customRxnsToIgnore?

[mihai-sysbio]

Hi @wshao1! Just judging by the documentation, I too find it unclear. I'm wondering if the explanation could be more simple, e.g. the reaction identifier is not in the right format. Could you share a bit of the sample code?

[wshao1]

Yes! As a bit of background, I am trying to create a GEM for the HCT116 cancer cell line, while protecting the reactions associated with essential genes. My code below is based on the ftINIT guide.

First I prepare Human1 for ftINIT by running a modified version of prepHumanModelForftINIT (prepHumanModelForftINIT.txt), where I add an additional argument essentialRxns to append reactions to customRxnsToIgnore.

prepDataHCT116 = prepHumanModelForftINIT(ihuman, false, human_tasks_path, human1_rxns_path, essentialRxns)
save('prepDataHCT116.mat', 'prepDataHCT116')

Then I proceed to run ftINIT using the resulting prepData structure.

%% Load DepMap transcriptomics data
hct116_tpms = readtable('../raw_data/hct116_tpms/hct116_tpms.csv')

%% Convert data into input structure
hct116_struct.genes = hct116_tpms{:, 1};    % gene names
hct116_struct.tissues = {'HCT116'};    % sample name
hct116_struct.levels = hct116_tpms{:, 2};    % gene TPM values
hct116_struct.threshold = 0.02;

%% Run ftINIT
hct116_model = ftINIT(prepDataHCT116, hct116_struct.tissues{1}, [], [], hct116_struct, {}, getHumanGEMINITSteps('1+1'), true, true);

I don't believe the issue stems from incorrect reaction identifiers, because I do see more reactions marked with a 1 in the prepData.toIgnoreCustomRxns field. The issue seems to be how ftINIT processes the toIgnoreCustomRxns field to protect those reactions from removal.

[mihai-sysbio]

Having had a deeper look, I'm actually confused that as params '1+1' is passed for the steps, which is the same as the default option but then that might not fully correspond with the steps that are required to ignore the custom reactions

RAVEN/INIT/ftINIT.m

Lines 168 to 183 in da86f62

	rxnsToIgnoreLastStep = [1;1;1;1;1;1;1;1];
	%We assume that all essential rxns are irrev - this is taken care of in
	%prepINITModel. We then use an initial flux "from last run" of 0.1 for all
	%reactions. This is used for knowing what flux should be forced through an
	%essential rxn.
	fluxes = ones(length(prepData.minModel.rxns), 1).*0.1;
	for initStep = 1:length(INITSteps)
	disp(['ftINIT: Running step ' num2str(initStep)])
	stp = INITSteps{initStep};
	if any ((rxnsToIgnoreLastStep - stp.RxnsToIgnoreMask) < 0)
	dispEM('RxnsToIgnoreMask may not cover rxns not covered in previous steps, but the other way around is fine.');
	end
	rxnsToIgnoreLastStep = stp.RxnsToIgnoreMask;

@johan-gson could you have a look at this?

[johan-gson]

Hi @wshao1 , customRxnsToIgnore does not do this, these are the reactions that are not included in the first step in the algorithm, so something different. I have not thought of this case, there is no direct support for it, it would eventually be a good thing to add as it would be a very simple fix. I'm very short on time right now, so I cannot engage in any development right now. If you want, you could modify the code to take care of this. A suggestion is to test the following: In prepINITModel, add two parameters:
function prepData = prepINITModel(origRefModel, taskStruct, spontRxnNames, convertGenes, customRxnsToIgnore, extComp, skipScaling, extraEssentialRxns, extraEssentialDirs)

you send in these extra reaction ids when you call the prepINITModel function. You need to decide on a direction in which these reactions are essential (you cannot easily do both). If irrev, then it is easy, always forward (1), otherwise it can also be backwards (-1)
ex.:
extraEssentialRxns = {'MARXXXX';'MARYYYY'}; % a vertical cell array, not horizontal, the reaction ids of the extra essential reactions
extraEssentialDirs = [1;-1]; %vertical vector

Then, on line 155 or somewhere around there, add these lines:
%Flip extra essential reactions with reversed essential direction
minModel1 = reverseRxns(minModel1, extraEssentialRxns(extraEssentialDirs < 0));
%Make the extra essential reactions irrev
sel = ismember(minModel1.rxns, extraEssentialDirs);
minModel1.rev(sel) = 0;
minModel1.lb(sel) = 0;

i.e., insert afterthese lines:
minModel1 = reverseRxns(minModel1, minModel1.rxns(selInd(essentialRevDir)));
%constructEquations(minModel2,minModel2.rxns(selInd(essentialRevDir))) %looks good

%Then make them irreversible
minModel1.rev(selInd) = 0;
minModel1.lb(selInd) = 0;

I think this should work, definitely worth a try! Let me know how it went if you try it! Note that the code does not mess up the model that you get out (making irreversible etc.), we are just changing the specialized model used in the MILP, so it is just internal for the algorithm!

Best of luck!