Push to master for pre-release v0.6.1

ComplementaryData/modelCuration/inactivationRxns.tsv

@@ -235,3 +235,5 @@ TAG_HSad	reaction is artificial and involves dead-end artificial metabolite with
 TAG_HSad_NE	reaction is artificial and involves dead-end artificial metabolite with no apparent purpose, and should therefore be DELETED
 TAG_HSad_E	reaction is artificial and involves dead-end artificial metabolite with no apparent purpose, and should therefore be DELETED
 CHOLESTle	reaction is artificial and involves dead-end artificial metabolite with no apparent purpose, and should therefore be DELETED
+r1330	reaction appears to be an attempt to fix ATP leakage, but not necessary with the updated handling of mitochondrial proton transport and should be DELETED
+r1331	reaction appears to be an attempt to fix ATP leakage, but not necessary with the updated handling of mitochondrial proton transport and should be DELETED

ComplementaryScripts/Functions/getCompNetwork.m

@@ -0,0 +1,35 @@
+function [compNetwork, I]=getCompNetwork(model,comp,includePartial)
+% getCompNetwork
+%   Gets the metabolic network for a specified compartment
+%
+% Input:
+%   model           a model structure
+%   comp            string with the compartment id
+%   includePartial  if true, include reactions with metabolites partially
+%                   present in the specified compartment (opt, default false)
+%
+% Output:
+%   compNetwork     a model structure for the specified compartment
+%   I               index of reactions in the specified compartment
+%
+% Usage: [compNetwork, I]=getCompNetwork(model,comp,includePartial)
+%
+% Hao Wang, 2018-11-07
+%
+
+if ischar(comp)
+    comp={comp};
+else
+    error('Incorrect compartment id!');
+end
+if nargin<3
+    includePartial=false;
+end
+
+% get the reaction list for the specified compartment
+I=getRxnsInComp(model,comp,includePartial);
+
+rxnToRemove=setdiff(transpose(1:length(model.rxns)), I);
+compNetwork=removeReactions(model,rxnToRemove,true,true,true);
+
+end

ComplementaryScripts/Functions/getSubNetworksInComp.m

@@ -0,0 +1,107 @@
+function [subNetworks]=getSubNetworksInComp(model,comp,metsToRemove,includePartial)
+% getSubNetworksInComp
+%	Check network connectivity in specified compartment and output
+%   the sub-networks in JSON fommat for map generation purpose
+%
+% Input:
+%   model           a model structure
+%   comp            string with the compartment id
+%   metsToRemove    cell array of mets to be excluded from the network
+%   includePartial  if true, include reactions with metabolites partially
+%                   present in the specified compartment (opt, default false)
+%
+% Output:
+%   subNetworks     a structure for the sub-networks
+%      id           the id for identifying the sub-networks
+%      number       number of reactions in each sub-network
+%      rxns         cell array of reaction identifiers in each sub-network
+%      subSystems   cell array of subSystems for reactions in each sub-network
+%
+%
+% NOTE: The output filename "SubNetworks.json" and in defined JSON format 
+%
+% Usage: [subNetworks]=getSubNetworksInComp(model,comp,metsToRemove,includePartial)
+%
+% Hao Wang, 2018-11-12
+%
+
+% handle input arguments
+if nargin < 2
+    error('Missing input!');
+end
+if nargin < 3
+    metsToRemove={''};
+end
+if nargin < 4
+    includePartial=false;
+end
+
+% get the subnetworks in specified compartment
+compNetwork = getCompNetwork(model,comp,includePartial);
+
+% refine the network by excluding any metabolites to be removed and
+% identifiy any ractions that have only currency metabolites
+excludeMets = intersect(metsToRemove, compNetwork.mets);
+if isempty(excludeMets)
+    frintf('The provided metabolites are not found in the model\n');
+    reducedNetwork = compNetwork;
+    currencyRxns = '';
+else
+    reducedNetwork = removeMets(compNetwork,excludeMets,0,1);
+    currencyRxns = setdiff(compNetwork.rxns, reducedNetwork.rxns);
+    [~, currencyRxnInd] = ismember(currencyRxns, model.rxns);
+end
+
+% get the sub-graphs
+subGraphs = getAllSubGraphs(reducedNetwork);
+
+% generate output
+graphNum = size(subGraphs,2);
+subNetworks.id = transpose(1:graphNum);
+
+% write in JSON format
+fid = fopen('Subnetworks.json','w');
+fprintf(fid, '{\n');
+
+for i=1:graphNum
+    index = find(subGraphs(:,i));
+    [~, rxnList]=find(reducedNetwork.S(index,:));
+    I = unique(rxnList);
+    subNetworks.number(i,1) = length(I);
+    subNetworks.rxns{i,1} = reducedNetwork.rxns(I);
+    subNetworks.subSystems{i,1} = reducedNetwork.subSystems(I);
+
+    fprintf(fid,['\t"' num2str(i) '":[\n']);
+    writeSubNetworks(fid, reducedNetwork.rxns(I), reducedNetwork.subSystems(I));
+    if i==graphNum && isempty(currencyRxns)
+        fprintf(fid, '\t]\n');
+    else
+        fprintf(fid, '\t],\n');
+    end
+end
+
+if ~isempty(currencyRxns)
+    fprintf(fid,['\t"currencyRxns":[\n']);
+    writeSubNetworks(fid, currencyRxns, model.subSystems(currencyRxnInd));
+    fprintf(fid, '\t]\n');
+end
+fprintf(fid, '}\n');
+fclose(fid);
+
+end
+
+
+function writeSubNetworks(file, rxnList, subSystemList)
+if ~isequal(numel(rxnList), numel(subSystemList))
+    error('Input reaction and subSystems lists do not match!');
+else
+    for i=1:numel(rxnList)
+        if i==numel(rxnList)
+            fprintf(file,['\t\t["' rxnList{i} '", "' subSystemList{i} '"]\n']);
+        else
+            fprintf(file,['\t\t["' rxnList{i} '", "' subSystemList{i} '"],\n']);
+        end
+    end
+end
+
+end

ComplementaryScripts/Functions/loadHumanGEM.m

@@ -0,0 +1,37 @@
+function model = loadHumanGEM
+% loadHumanGEM
+%   Load humanGEM and prepare for a simulation-ready model by contraining 
+%   problematic reactions archived in incactivationRxns.tsv
+%
+% Usage: model = loadHumanGEM
+%
+% Hao Wang, 2018-11-14
+%
+
+% get model path
+[ST, I]=dbstack('-completenames');
+modelPath=fileparts(fileparts(fileparts(ST(I).file)));
+
+% load model
+matFile=fullfile(modelPath,'ModelFiles','mat','humanGEM.mat');
+load(matFile);
+
+% get reactions need to be constrained
+rxnsToConstrain = '';
+inactivateRxnsFile=fullfile(modelPath,'ComplementaryData','modelCuration','inactivationRxns.tsv');
+if exist(inactivateRxnsFile, 'file') == 2
+    fid = fopen(inactivateRxnsFile,'r');
+    input = textscan(fid,'%s %s','Delimiter','\t','Headerlines',1);
+    fclose(fid);
+    rxnsToConstrain = input{1};
+end
+
+% conduct constraining, if any
+if ~isempty(rxnsToConstrain)
+    model = setParam(ihuman, 'eq', rxnsToConstrain, 0);
+else
+    model = ihuman;
+end
+
+end
+