applyComplexData

 applyComplexData
   Apply stochiometry for complex in an ecModel

 Input:
   model           an ecModel in GECKO 3 format (with ecModel.ec structure)
   complexInfo     structure as generated by getComplexData. If nothing
                   is provided, an attempt will be made to read
                   data/ComplexPortal.json from the obj.params.path folder
                   specified in the modelAdapter.
   modelAdapter    a loaded model adapter (Optional, will otherwise use the
                   default model adapter).
   verbose         logical if a summary should be shown in the Command
                   Window (Optional, default true)

 Output:
   model           ecModel where model.ec.rxnEnzMat is populated with
                   subunit stochiometries
   foundComplex    complexes that fully matched between the model and the
                   complex data
   proposedComplex complexes where the model contained >75% but <100% of
                   the proteins indicated by Complex Portal, or where the
                   model contained more proteins than indicated for that
                   complex in Complex Portal.

 Usage:
   [model, foundComplex, proposedComplex] = applyComplexData(ecModel, complexInfo, modelAdapter);

0001 function [model, foundComplex, proposedComplex] = applyComplexData(model, complexInfo, modelAdapter, verbose)
0002 % applyComplexData
0003 %   Apply stochiometry for complex in an ecModel
0004 %
0005 % Input:
0006 %   model           an ecModel in GECKO 3 format (with ecModel.ec structure)
0007 %   complexInfo     structure as generated by getComplexData. If nothing
0008 %                   is provided, an attempt will be made to read
0009 %                   data/ComplexPortal.json from the obj.params.path folder
0010 %                   specified in the modelAdapter.
0011 %   modelAdapter    a loaded model adapter (Optional, will otherwise use the
0012 %                   default model adapter).
0013 %   verbose         logical if a summary should be shown in the Command
0014 %                   Window (Optional, default true)
0015 %
0016 % Output:
0017 %   model           ecModel where model.ec.rxnEnzMat is populated with
0018 %                   subunit stochiometries
0019 %   foundComplex    complexes that fully matched between the model and the
0020 %                   complex data
0021 %   proposedComplex complexes where the model contained >75% but <100% of
0022 %                   the proteins indicated by Complex Portal, or where the
0023 %                   model contained more proteins than indicated for that
0024 %                   complex in Complex Portal.
0025 %
0026 % Usage:
0027 %   [model, foundComplex, proposedComplex] = applyComplexData(ecModel, complexInfo, modelAdapter);
0028 
0029 if nargin < 4 || isempty(verbose)
0030     verbose = true;
0031 end
0032 
0033 if nargin < 3 || isempty(modelAdapter)
0034     modelAdapter = ModelAdapterManager.getDefault();
0035     if isempty(modelAdapter)
0036         error('Either send in a modelAdapter or set the default model adapter in the ModelAdapterManager.')
0037     end
0038 end
0039 params = modelAdapter.params;
0040 
0041 if nargin<2 || isempty(complexInfo)
0042     complexInfo = fullfile(params.path,'data','ComplexPortal.json');
0043     if ~isfile(complexInfo)
0044         complexInfo = getComplexData([], modelAdapter);
0045     end
0046 end
0047 
0048 if ischar(complexInfo) || isstring(complexInfo)
0049     jsonStr = fileread(complexInfo);
0050     complexData = jsondecode(jsonStr);
0051 else
0052     complexData = complexInfo;
0053 end
0054 
0055 %Remove prefixes on rxn names for gecko light
0056 if ~model.ec.geckoLight
0057     rxnNames = model.ec.rxns;
0058 else
0059     rxnNames = extractAfter(model.ec.rxns,4);
0060 end
0061 
0062 foundComplex = cell(numel(rxnNames),7);
0063 foundComplexCount = 0;
0064 proposedComplex = cell(numel(rxnNames)*2,8); % Might propose more
0065 proposedComplexCount = 0;
0066 
0067 complexProts = repmat({''}, numel({complexData.complexID})*10, 1);
0068 complexMatrix = zeros(numel({complexData.complexID}),numel(complexProts));
0069 lastProt = 0;
0070 for i=1:size(complexMatrix,1)
0071     % Locate proteins (and add if new)
0072     protIDs     = complexData(i).protID;
0073     [isOld, protMatIdx] = ismember(protIDs,complexProts);
0074     newProts    = protIDs(~isOld);
0075     complexProts(lastProt+1:lastProt+numel(newProts))=newProts;
0076     protMatIdx(~isOld)  = lastProt+1 : lastProt+numel(newProts);
0077     lastProt = lastProt + numel(newProts);
0078     % Add subunit numbers
0079     if all(complexData(i).stochiometry == 0)
0080         complexData(i).stochiometry(:) = 1;
0081     end
0082     complexMatrix(i,protMatIdx) = complexData(i).stochiometry;
0083 end
0084 
0085 % Remove excess space
0086 complexMatrix(:,lastProt+1:end) = [];
0087 complexProts(lastProt+1:end) = [];
0088 progressbar('Assign complexes to reactions')
0089 for i = 1:numel(rxnNames)
0090     % Get the proteins from the model
0091     modelProts = model.ec.enzymes(find(model.ec.rxnEnzMat(i,:)));
0092     [protsInMat, protsIdx] = ismember(modelProts,complexProts);
0093     if ~any(protsInMat)
0094         continue
0095     end
0096     % Find complexes that includes the model proteins
0097     [potComplex,~] = find(complexMatrix(:,protsIdx(protsInMat)));
0098     potComplex = unique(potComplex);
0099     % Some stats on number of subunits, percentage match etc.
0100     matchComplexUnits = sum(logical(complexMatrix(potComplex,protsIdx(protsInMat))),2);
0101     totalComplexUnits = sum(logical(complexMatrix(potComplex,:)),2);
0102     modComplexUnits = numel(modelProts);
0103     percMatch  = matchComplexUnits/modComplexUnits;
0104     totalMatch = totalComplexUnits/modComplexUnits;
0105     if any(percMatch == 1 & totalMatch == 1)
0106         % Exact match
0107         complexMatch = find(percMatch == 1 & totalMatch == 1);
0108         if numel(complexMatch)>1 % If more than one match, take first match
0109             complexMatch = complexMatch(1);
0110         end
0111         complexMatch = potComplex(complexMatch);
0112         foundComplexCount = foundComplexCount + 1;
0113         foundComplex{foundComplexCount,1} = model.ec.rxns{i};
0114         foundComplex{foundComplexCount,2} = complexData(complexMatch).complexID;
0115         foundComplex{foundComplexCount,3} = complexData(complexMatch).name;
0116         foundComplex{foundComplexCount,4} = complexData(complexMatch).geneName;
0117         foundComplex{foundComplexCount,5} = modelProts;
0118         foundComplex{foundComplexCount,6} = complexData(complexMatch).protID;
0119         foundComplex{foundComplexCount,7} = complexData(complexMatch).stochiometry;
0120 
0121         % Apply in rxnEnzMat
0122         [~, modelProtsIdx] = ismember(modelProts, model.ec.enzymes);
0123         model.ec.rxnEnzMat(i,modelProtsIdx) = complexData(complexMatch).stochiometry;
0124     
0125     elseif modComplexUnits > 1 % Only suggest potential complex if grRule > 1 protein
0126         moreUnitsInData = find(percMatch == 1 & totalMatch > 1);
0127         if any(moreUnitsInData)
0128             % All ecModel subunits match, but Complex Portal has more subunits
0129             % Take the complex with the minimum number of extra subunits
0130             [~,propComplex] = min(totalMatch(moreUnitsInData));
0131             propComplex     = moreUnitsInData(propComplex);
0132             selectComplex   = potComplex(propComplex);
0133             proposedComplexCount = proposedComplexCount + 1;
0134             proposedComplex{proposedComplexCount,1} = model.ec.rxns{i};
0135             proposedComplex{proposedComplexCount,2} = complexData(selectComplex).complexID;
0136             proposedComplex{proposedComplexCount,3} = complexData(selectComplex).name;
0137             proposedComplex{proposedComplexCount,4} = complexData(selectComplex).geneName;
0138             proposedComplex{proposedComplexCount,5} = modelProts;
0139             proposedComplex{proposedComplexCount,6} = complexData(selectComplex).protID;
0140             proposedComplex{proposedComplexCount,7} = complexData(selectComplex).stochiometry;
0141             proposedComplex{proposedComplexCount,8} = totalMatch(propComplex)*100;
0142         end
0143         otherUnits = find(percMatch >= 0.75 & percMatch < 1 & totalMatch <= 1);
0144         if any(otherUnits)
0145             % At least 75% of ecModels subunits match, Complex Portal maybe have less or different subunits
0146             % Take the complex with highest % match.
0147             [~,propComplex] = max(percMatch(otherUnits));
0148             propComplex     = otherUnits(propComplex);
0149             selectComplex   = potComplex(propComplex);
0150             proposedComplexCount = proposedComplexCount + 1;
0151             proposedComplex{proposedComplexCount,1} = model.ec.rxns{i};
0152             proposedComplex{proposedComplexCount,2} = complexData(selectComplex).complexID;
0153             proposedComplex{proposedComplexCount,3} = complexData(selectComplex).name;
0154             proposedComplex{proposedComplexCount,4} = complexData(selectComplex).geneName;
0155             proposedComplex{proposedComplexCount,5} = modelProts;
0156             proposedComplex{proposedComplexCount,6} = complexData(selectComplex).protID;
0157             proposedComplex{proposedComplexCount,7} = complexData(selectComplex).stochiometry;
0158             proposedComplex{proposedComplexCount,8} = percMatch(propComplex)*100;
0159         end
0160     end
0161     progressbar(i/numel(rxnNames))    
0162 end
0163 progressbar(1) % Make sure it closes
0164 
0165 % Remove empty space
0166 foundComplex(foundComplexCount+1:end,:) = [];
0167 proposedComplex(proposedComplexCount+1:end,:) = [];
0168 
0169 rowHeadings = {'rxn', 'complexID','name','genes','protID_model','protID_complex','stochiometry'};
0170 
0171 foundComplex = cell2table(foundComplex, 'VariableNames', rowHeadings);
0172 
0173 proposedComplex = cell2table(proposedComplex, 'VariableNames', [rowHeadings 'match']);
0174 if verbose
0175     disp(['A total of ' int2str(numel(foundComplex(:,1))) ' complex have full match, and ' int2str(numel(proposedComplex(:,1))) ' proposed.'])
0176 end
0177 end