printModelStats

 printModelStats
   prints some statistics about a model to the screen

   model               a model structure
   printModelIssues    true if information about unconnected
                       reactions/metabolites and elemental balancing
                       should be printed (optional, default false)
   printDetails        true if detailed information should be printed
                       about model issues. Only used if printModelIssues
                       is true (optional, default true)

 Usage: printModelStats(model,printModelIssues, printDetails)

0001 function printModelStats(model, printModelIssues, printDetails)
0002 % printModelStats
0003 %   prints some statistics about a model to the screen
0004 %
0005 %   model               a model structure
0006 %   printModelIssues    true if information about unconnected
0007 %                       reactions/metabolites and elemental balancing
0008 %                       should be printed (optional, default false)
0009 %   printDetails        true if detailed information should be printed
0010 %                       about model issues. Only used if printModelIssues
0011 %                       is true (optional, default true)
0012 %
0013 % Usage: printModelStats(model,printModelIssues, printDetails)
0014 
0015 if nargin<2
0016     printModelIssues=false;
0017 end
0018 if nargin<3
0019     printDetails=true;
0020 end
0021 
0022 fprintf(['Network statistics for ' model.id ': ' model.name '\n']);
0023 
0024 %Get which reactions are present in each compartment
0025 rxnComps=sparse(numel(model.rxns),numel(model.comps));
0026 
0027 %For each compartment, find the metabolites that are present in that
0028 %compartment and then the reactions they are involved in
0029 for i=1:numel(model.comps)
0030     [~, I]=find(model.S(model.metComps==i,:));
0031     rxnComps(I,i)=1;
0032 end
0033 
0034 if isfield(model,'eccodes')
0035     fprintf(['EC-numbers\t\t\t' num2str(numel(unique(model.eccodes))) '\n']);
0036 end
0037 
0038 %Print information about genes
0039 if isfield(model,'genes')
0040     fprintf(['Genes*\t\t\t\t' num2str(numel(model.genes)) '\n']);
0041     %Find the genes in each compartment
0042     for i=1:numel(model.comps)
0043         [~, I]=find(model.rxnGeneMat(rxnComps(:,i)==1,:));
0044         fprintf(['\t' model.compNames{i} '\t' num2str(numel(unique(I))) '\n']);
0045     end
0046 end
0047 
0048 %Print information about reactions
0049 fprintf(['\nReactions*\t\t\t' num2str(numel(model.rxns)) '\n']);
0050 for i=1:numel(model.comps)
0051     fprintf(['\t' model.compNames{i} '\t' num2str(sum(rxnComps(:,i))) '\n']);
0052 end
0053 
0054 %Removes the effect of compartments and removes duplicate reactions
0055 temp=model;
0056 temp.comps(:)={'s'}; %Set all compartments to be the same
0057 equ=constructEquations(sortModel(temp,true,true),temp.rxns,false);
0058 
0059 fprintf(['Unique reactions**\t' num2str(numel(unique(equ))) '\n']);
0060 
0061 %Print information about metabolites
0062 fprintf(['\nMetabolites\t\t\t' num2str(numel(model.mets)) '\n']);
0063 for i=1:numel(model.comps)
0064     fprintf(['\t' model.compNames{i} '\t' num2str(sum(model.metComps==i)) '\n']);
0065 end
0066 
0067 fprintf(['Unique metabolites\t' num2str(numel(unique(model.metNames))) '\n']);
0068 
0069 fprintf('\n* Genes and reactions are counted for each compartment if any of the corresponding metabolites are in that compartment. The sum may therefore not add up to the total number.\n');
0070 fprintf('** Unique reactions are defined as being biochemically unique (no compartmentalization)\n');
0071 
0072 %Also print some potential problems if there are any
0073 if printModelIssues==true
0074     fprintf(['\nShort model quality summary for ' model.id ': ' model.name '\n']);
0075     
0076     %Check that all the metabolites are being used
0077     involvedMat=model.S;
0078     involvedMat(involvedMat~=0)=1;
0079     usedMets=sum(involvedMat,2);
0080     notPresent=find(usedMets==0);
0081     if ~isempty(notPresent)
0082         errorText=['Non-used metabolites\t' num2str(numel(notPresent)) '\n'];
0083         if printDetails==true
0084             for i=1:numel(notPresent)
0085                 errorText=[errorText '\t(' model.mets{notPresent(i)} ') ' model.metNames{notPresent(i)} '\n'];
0086             end
0087             errorText=[errorText '\n'];
0088         end
0089         fprintf(errorText);
0090     end
0091     
0092     %Check if there are empty reactions
0093     usedRxns=sum(involvedMat,1);
0094     notUsed=find(usedRxns==0);
0095     if ~isempty(notUsed)
0096         errorText=['Empty reactions\t' num2str(numel(notUsed)) '\n'];
0097         if printDetails==true
0098             for i=1:numel(notUsed)
0099                 errorText=[errorText '\t' model.rxns{notUsed(i)} '\n'];
0100             end
0101             errorText=[errorText '\n'];
0102         end
0103         fprintf(errorText);
0104     end
0105     
0106     %Check if there are dead-end reactions/metabolites
0107     [~, deletedReactions, deletedMetabolites]=simplifyModel(model,true,false,false,true);
0108     
0109     if ~isempty(deletedReactions)
0110         errorText=['Dead-end reactions\t' num2str(numel(deletedReactions)) '\n'];
0111         if printDetails==true
0112             for i=1:numel(deletedReactions)
0113                 errorText=[errorText '\t' deletedReactions{i} '\n'];
0114             end
0115             errorText=[errorText '\n'];
0116         end
0117         fprintf(errorText);
0118     end
0119     
0120     %Ignore non-used metabolites
0121     deletedMetabolites=setdiff(deletedMetabolites,model.mets(notPresent));
0122     %Must map to indexes in order to print names
0123     deletedMetabolites=find(ismember(model.mets,deletedMetabolites));
0124     if ~isempty(deletedMetabolites)
0125         errorText=['Dead-end metabolites\t' num2str(numel(deletedMetabolites)) '\n'];
0126         if printDetails==true
0127             for i=1:numel(deletedMetabolites)
0128                 errorText=[errorText '\t(' model.mets{deletedMetabolites(i)} ') ' model.metNames{deletedMetabolites(i)} '\n'];
0129             end
0130             errorText=[errorText '\n'];
0131         end
0132         fprintf(errorText);
0133     end
0134     
0135     balanceStructure=getElementalBalance(model);
0136     
0137     notParsed=find(balanceStructure.balanceStatus<0);
0138     notBalanced=find(balanceStructure.balanceStatus==0);
0139     
0140     if ~isempty(notParsed)
0141         errorText=['Reactions which could not be elementally balanced\t' num2str(numel(notParsed)) '\n'];
0142         if printDetails==true
0143             for i=1:numel(notParsed)
0144                 errorText=[errorText '\t' model.rxns{notParsed(i)} '\n'];
0145             end
0146             errorText=[errorText '\n'];
0147         end
0148         fprintf(errorText);
0149     end
0150     if ~isempty(notBalanced)
0151         errorText=['Reactions which are elementally unbalanced\t' num2str(numel(notBalanced)) '\n'];
0152         if printDetails==true
0153             names=strcat(balanceStructure.elements.names,{', '});
0154             for i=1:numel(notBalanced)
0155                 badOnes=sprintf('%s', names{abs(balanceStructure.leftComp(notBalanced(i),:)-balanceStructure.rightComp(notBalanced(i),:))>10^-7});
0156                 errorText=[errorText '\t' model.rxns{notBalanced(i)} '\t' badOnes(1:end-2) '\n'];
0157             end
0158             errorText=[errorText '\n'];
0159         end
0160         fprintf(errorText);
0161     end
0162 end
0163 end