Skip to content

Commit

Permalink
style: all RAVEN functions
Browse files Browse the repository at this point in the history 
-applied MATLAB smart indent (EDITOR -> EDIT -> Indent -> Smart)
-rewrapped the code comments to the new indents (EDITOR -> EDIT -> Comment -> Wrap Comments)
-fixed several cases when function did not have the ending line
  • Loading branch information
Simonas Marcišauskas committed Apr 16, 2018
1 parent 7456bd5 commit 5690fe4
Show file tree
Hide file tree
Showing 143 changed files with 5,150 additions and 5,120 deletions.
66 changes: 33 additions & 33 deletions INIT/getINITModel.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,7 +134,7 @@

%Create the task structure if not supplied
if any(taskFile) && isempty(taskStructure)
taskStructure=parseTaskList(taskFile);
taskStructure=parseTaskList(taskFile);
end

if printReport==true
Expand All @@ -156,7 +156,7 @@
fprintf('-Using metabolic tasks\n');
end
fprintf('\n');

printScores(refModel,'Reference model statistics',hpaData,arrayData,tissue,celltype);
end

Expand All @@ -175,9 +175,9 @@
%reactions
if ~isempty(taskStructure)
[taskReport, essentialRxnMat]=checkTasks(cModel,[],printReport,true,true,taskStructure);

essentialRxnsForTasks=cModel.rxns(any(essentialRxnMat,2));

%Remove tasks that cannot be performed
taskStructure(taskReport.ok==false)=[];
if printReport==true
Expand All @@ -192,11 +192,11 @@

%Run the INIT algorithm. The exchange reactions that are used in the final
%reactions will be open, which doesn't fit with the last step. Therefore
%delete reactions from the original model instead of taking the output.
%The default implementation does not constrain reversible reactions to only
%carry flux in one direction.
%Runs without the constraints on reversibility and with all output allowed.
%This is to reduce the complexity of the problem.
%delete reactions from the original model instead of taking the output. The
%default implementation does not constrain reversible reactions to only
%carry flux in one direction. Runs without the constraints on reversibility
%and with all output allowed. This is to reduce the complexity of the
%problem.
[~, deletedRxnsInINIT, metProduction]=runINIT(simplifyModel(cModel),rxnScores,metabolomicsData,essentialRxnsForTasks,0,true,false,params);
initModel=removeReactions(cModel,deletedRxnsInINIT,true,true);
if printReport==true
Expand All @@ -217,7 +217,7 @@
%Remove exchange reactions and reactions already included in the INIT
%model
refModelNoExc=removeReactions(refModel,union(initModel.rxns,getExchangeRxns(refModel)),true,true);

%At this stage the model is fully connected and most of the genes with
%good scores should have been included. The final gap-filling should
%take the scores of the genes into account, so that "rather bad"
Expand All @@ -235,7 +235,7 @@
printScores(outModel,'Functional model statistics',hpaData,arrayData,tissue,celltype);
printScores(removeReactions(outModel,intersect(outModel.rxns,initModel.rxns),true,true),'Reactions added to perform the tasks',hpaData,arrayData,tissue,celltype);
end

addedRxnsForTasks=refModelNoExc.rxns(any(addedRxnMat,2));
else
outModel=initModel;
Expand All @@ -258,19 +258,19 @@
for i=1:numel(model.rxns)
ids=find(model.rxnGeneMat(i,:));
if numel(ids)>1
scores=geneScores(I(ids));
%Only keep the positive ones if possible
model.rxnGeneMat(i,ids(~(scores>0 | scores==max(scores))))=0;
scores=geneScores(I(ids));
%Only keep the positive ones if possible
model.rxnGeneMat(i,ids(~(scores>0 | scores==max(scores))))=0;
end
%Rewrite the grRules to be only OR
if isfield(model,'grRules')
J=find(model.rxnGeneMat(i,:));
for j=1:numel(J)
model.grRules{i}=[model.grRules{i} '(' model.genes{J(j)} ')'];
if j<numel(J)
model.grRules{i}=[model.grRules{i} ' or '];
end
end
J=find(model.rxnGeneMat(i,:));
for j=1:numel(J)
model.grRules{i}=[model.grRules{i} '(' model.genes{J(j)} ')'];
if j<numel(J)
model.grRules{i}=[model.grRules{i} ' or '];
end
end
end
end

Expand All @@ -291,10 +291,10 @@
model.geneComps(I)=[];
end

%At this stage the model will contain some exchange reactions but probably not all
%(and maybe zero). This can be inconvenient, so all exchange reactions from the
%reference model are added, except for those which involve metabolites that
%are not in the model.
%At this stage the model will contain some exchange reactions but probably
%not all (and maybe zero). This can be inconvenient, so all exchange
%reactions from the reference model are added, except for those which
%involve metabolites that are not in the model.

%First delete and included exchange reactions in order to prevent the order
%from changing
Expand Down Expand Up @@ -343,12 +343,12 @@

%This is for printing a summary of a model
function [rxnS, geneS]=printScores(model,name,hpaData,arrayData,tissue,celltype)
[a, b]=scoreModel(model,hpaData,arrayData,tissue,celltype);
rxnS=mean(a);
geneS=mean(b(~isinf(b)));
fprintf([name ':\n']);
fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']);
fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']);
fprintf(['\tMean gene score: ' num2str(geneS) '\n']);
fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']);
[a, b]=scoreModel(model,hpaData,arrayData,tissue,celltype);
rxnS=mean(a);
geneS=mean(b(~isinf(b)));
fprintf([name ':\n']);
fprintf(['\t' num2str(numel(model.rxns)) ' reactions, ' num2str(numel(model.genes)) ' genes\n']);
fprintf(['\tMean reaction score: ' num2str(rxnS) '\n']);
fprintf(['\tMean gene score: ' num2str(geneS) '\n']);
fprintf(['\tReactions with positive scores: ' num2str(100*sum(a>0)/numel(a)) '%%\n\n']);
end
72 changes: 34 additions & 38 deletions INIT/runINIT.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@
end

if numel(presentMets)~=numel(unique(presentMets))
EM='Duplicate metabolite names in presentMets';
EM='Duplicate metabolite names in presentMets';
dispEM(EM);
end

Expand All @@ -123,13 +123,13 @@
dispEM(EM,false);
end

%The irreversible reactions that are essential must have a flux and are therefore not
%optimized for using MILP, which reduces the problem size. However, reversible
%reactions must have a flux in one direction, so they have to stay in
%the problem. The essentiality constraint on reversible reactions is
%implemented in the same manner as for reversible reactions when
%noRevLoops==true, but with the additional constraint that C ub=-1. This
%forces one of the directions to be active.
%The irreversible reactions that are essential must have a flux and are
%therefore not optimized for using MILP, which reduces the problem size.
%However, reversible reactions must have a flux in one direction, so they
%have to stay in the problem. The essentiality constraint on reversible
%reactions is implemented in the same manner as for reversible reactions
%when noRevLoops==true, but with the additional constraint that C ub=-1.
%This forces one of the directions to be active.
revRxns=find(model.rev~=0);
essentialReversible=find(ismember(model.rxns(revRxns),essentialRxns));
essentialRxns=intersect(essentialRxns,model.rxns(model.rev==0));
Expand Down Expand Up @@ -164,7 +164,7 @@
metsToAdd.mets=strcat({'FAKEFORPM'},num2str(pmIndexes));
metsToAdd.metNames=metsToAdd.mets;
metsToAdd.compartments=irrevModel.comps{1};

%There is no constraints on the metabolites yet, since maybe not all of
%them could be produced
irrevModel=addMets(irrevModel,metsToAdd);
Expand All @@ -174,16 +174,16 @@
for i=1:numel(pmIndexes)
%Get the matching mets
I=ismember(irrevModel.metNames,presentMets(pmIndexes(i)));

%Find the reactions where any of them are used.
[~, K, L]=find(irrevModel.S(I,:));

%This ugly loop is to avoid problems if a metabolite occurs several
%times in one reaction
KK=unique(K);
LL=zeros(numel(KK),1);
for j=1:numel(KK)
LL(j)=sum(L(K==KK(j)));
LL(j)=sum(L(K==KK(j)));
end
irrevModel.S(numel(irrevModel.mets)-numel(pmIndexes)+i,KK)=LL;
end
Expand All @@ -196,7 +196,8 @@
nonEssentialIndex=setdiff(1:nRxns,essentialIndex);
S=irrevModel.S;

%Add so that each non-essential reaction produces one unit of a fake metabolite
%Add so that each non-essential reaction produces one unit of a fake
%metabolite
temp=sparse(1:nRxns,1:nRxns,1);
temp(essentialIndex,:)=[];
S=[S;temp];
Expand All @@ -220,31 +221,25 @@
%Add constraints so that reversible reactions can only be used in one
%direction. This is done by adding the fake metabolites A, B, C for each
%reversible reaction in the following manner
% forward: A + .. => ...
% backwards: B + ... => ...
% int1: C => 1000 A
% int2: C => 1000 B
% A ub=999.9
% B ub=999.9
% C lb=-1
% int1 and int2 are binary
% forward: A + .. => ... backwards: B + ... => ... int1: C => 1000 A int2:
% C => 1000 B A ub=999.9 B ub=999.9 C lb=-1 int1 and int2 are binary
if any(forwardIndexes)
nRevBounds=numel(forwardIndexes);

%Add the A metabolites for the forward reactions and the B
%metabolites for the reverse reactions
%Add the A metabolites for the forward reactions and the B metabolites
%for the reverse reactions
I=speye(numel(irrevModel.rxns))*-1;
temp=[I(forwardIndexes,:);I(backwardIndexes,:)];

%Padding
temp=[temp sparse(size(temp,1),size(S,2)-numel(irrevModel.rxns))];

%Add the int1 & int2 reactions that produce A and B
temp=[temp speye(nRevBounds*2)*1000];

%And add that they also consume C
temp=[temp;[sparse(nRevBounds,size(S,2)) speye(nRevBounds)*-1 speye(nRevBounds)*-1]];

%Add the new reactions and metabolites
S=[S sparse(size(S,1),nRevBounds*2)];
S=[S;temp];
Expand All @@ -253,23 +248,25 @@
end

%Add so that the essential reactions must have a small flux and that the
%binary ones (and net-production reactions) may have zero flux. The
%integer reactions for reversible reactions have [0 1]
%binary ones (and net-production reactions) may have zero flux. The integer
%reactions for reversible reactions have [0 1]
prob.blx=[irrevModel.lb;zeros(nNonEssential+nNetProd+nRevBounds*2,1)];
prob.blx(essentialIndex)=max(0.1,prob.blx(essentialIndex));

%Add so that the binary ones and net-production reactions can have at the most flux 1.0
%Add so that the binary ones and net-production reactions can have at the
%most flux 1.0
prob.bux=[irrevModel.ub;ones(nNonEssential+nNetProd+nRevBounds*2,1)];

%Add that the fake metabolites must be produced in a small amount and that
%the A and B metabolites for reversible reactions can be [0 999.9] and C
%metabolites [-1 0]
prob.blc=[irrevModel.b(:,1);ones(nNonEssential,1);zeros(nRevBounds*2,1);ones(nRevBounds,1)*-1];

%Add that normal metabolites can be freely excreted if allowExcretion==true,
%and that the fake ones can be excreted 1000 units at most. C metabolites
%for essential reversible reactions should have an upper bound of -1.
%If noRevLoops is false, then add this constraint for all the reactions instead.
%Add that normal metabolites can be freely excreted if
%allowExcretion==true, and that the fake ones can be excreted 1000 units at
%most. C metabolites for essential reversible reactions should have an
%upper bound of -1. If noRevLoops is false, then add this constraint for
%all the reactions instead.
if noRevLoops==true
revUB=zeros(nRevBounds,1);
revUB(essentialReversible)=-1;
Expand All @@ -284,9 +281,8 @@
prob.buc=[metUB;ones(nNonEssential,1)*1000;ones(nRevBounds*2,1)*999.9;revUB];

%Add objective coefficients for the binary reactions. The negative is used
%since we're minimizing. The negative is taken for the prodWeight
%as well, in order to be consistent with the syntax that positive scores
%are good
%since we're minimizing. The negative is taken for the prodWeight as well,
%in order to be consistent with the syntax that positive scores are good
prob.c=[zeros(nRxns,1);rxnScores;ones(nNetProd,1)*prodWeight*-1;zeros(nRevBounds*2,1)];
prob.a=S;

Expand Down
57 changes: 29 additions & 28 deletions core/FSEOF.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,42 +55,43 @@
for i=1:iterations
n=i*targetMax/iterations;
model=setParam(model,'lb',targetRxn,n);

sol=solveLP(model,1);

fseof.results(:,i)=sol.x;

%Loop through all fluxes and identify the ones that increased upon the enforced objective flux

%Loop through all fluxes and identify the ones that increased upon the
%enforced objective flux
for j=1:length(fseof.results)
if fseof.results(j,1) > 0 %Check the positive fluxes

if i == 1 %The initial round
rxnDirection(j,1)=1;
fseof.target(j,1)=1;
if fseof.results(j,1) > 0 %Check the positive fluxes
if i == 1 %The initial round
rxnDirection(j,1)=1;
fseof.target(j,1)=1;
else

if (fseof.results(j,i) > fseof.results(j,i-1)) & fseof.target(j,1)
fseof.target(j,1)=1;
else
fseof.target(j,1)=0;
end
end

elseif fseof.results(j,1) < 0 %Check the negative fluxes

if i == 1 %The initial round
rxnDirection(j,1)=-1;
fseof.target(j,1)=1;
fseof.target(j,1)=1;
else
fseof.target(j,1)=0;
end
end
elseif fseof.results(j,1) < 0 %Check the negative fluxes
if i == 1 %The initial round
rxnDirection(j,1)=-1;
fseof.target(j,1)=1;
else
if (fseof.results(j,i) < fseof.results(j,i-1)) & fseof.target(j,1)
fseof.target(j,1)=1;
else
fseof.target(j,1)=0;
end
end

fseof.target(j,1)=1;
else
fseof.target(j,1)=0;
end
end
end

end
end

Expand Down Expand Up @@ -128,4 +129,4 @@
targets.logical=logical(fseof.target);
targets.slope=abs(fseof.results(:,iterations)-fseof.results(:,1))/abs(targetMax-targetMax/iterations); %Slope calculation
end
end
end
4 changes: 2 additions & 2 deletions core/addExchangeRxns.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,8 +75,8 @@
model.rxnGeneMat=[model.rxnGeneMat;sparse(numel(J),numel(model.genes))];
end
if isfield(model,'rxnComps')
model.rxnComps=[model.rxnComps;ones(numel(J),1)];
fprintf('NOTE: The exchange reactions are assigned to the first compartment\n');
model.rxnComps=[model.rxnComps;ones(numel(J),1)];
fprintf('NOTE: The exchange reactions are assigned to the first compartment\n');
end
if isfield(model,'rxnNotes')
model.rxnNotes=[model.rxnNotes;filler];
Expand Down
Loading

0 comments on commit 5690fe4

Please sign in to comment.