0001 function [rxnScores geneScores hpaScores arrayScores]=scoreModel(model,hpaData,arrayData,tissue,celltype,noGeneScore,multipleGeneScoring,multipleCellScoring,hpaLevelScores)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054 if nargin<3
0055 arrayData=[];
0056 end
0057 if nargin<5
0058 celltype=[];
0059 end
0060 if nargin<6
0061 noGeneScore=-2;
0062 end
0063 if nargin<7
0064 multipleGeneScoring='best';
0065 end
0066 if nargin<8
0067 multipleCellScoring='best';
0068 end
0069 if nargin<9
0070
0071 hpaLevelScores.names={'High' 'Medium' 'Low' 'None' 'Strong' 'Moderate' 'Weak' 'Negative'};
0072 hpaLevelScores.scores=[20 15 10 -8 20 15 10 -8];
0073 end
0074
0075 if isempty(hpaData) && isempty(arrayData)
0076 throw(MException('','Must supply hpaData, arrayData or both'));
0077 end
0078 if ~strcmpi(multipleGeneScoring,'best') && ~strcmpi(multipleGeneScoring,'average')
0079 throw(MException('','Valid options for multipleGeneScoring are "best" or "average"'));
0080 end
0081 if ~strcmpi(multipleCellScoring,'best') && ~strcmpi(multipleCellScoring,'average')
0082 throw(MException('','Valid options for multipleCellScoring are "best" or "average"'));
0083 end
0084
0085
0086
0087 if any(arrayData)
0088 if numel(arrayData.tissues)<2
0089 throw(MException('','arrayData must contain measurements for at least two celltypes/tissues since the score is calculated based on the expression level compared to the overall average'));
0090 end
0091 end
0092
0093
0094
0095 if isempty(arrayData)
0096 arrayData.genes={};
0097 arrayData.tissues={};
0098 arrayData.celltypes={};
0099 arrayData.levels=[];
0100 end
0101 if isempty(hpaData)
0102 hpaData.genes={};
0103 hpaData.tissues={};
0104 hpaData.celltypes={};
0105 hpaData.levels={};
0106 hpaData.types={};
0107 hpaData.reliabilities={};
0108 hpaData.gene2Level=[];
0109 hpaData.gene2Type=[];
0110 hpaData.gene2Reliability=[];
0111 end
0112
0113
0114 if ~ismember(upper(tissue),upper(hpaData.tissues)) && ~ismember(upper(tissue),upper(arrayData.tissues))
0115 throw(MException('','The tissue name does not match'));
0116 end
0117 if any(celltype)
0118
0119 if ~isfield(hpaData,'celltypes') || ~isfield(arrayData,'celltypes')
0120 throw(MException('','Both hpaData and arrayData must contain cell type information if cell type is to be used'));
0121 end
0122 if ~ismember(upper(celltype),upper(hpaData.celltypes)) && ~ismember(upper(celltype),upper(arrayData.celltypes))
0123 throw(MException('','The cell type name does not match'));
0124 end
0125 end
0126
0127
0128
0129 J=~strcmpi(hpaData.tissues,tissue);
0130
0131
0132 if any(celltype)
0133 J=J | ~strcmpi(hpaData.celltypes,celltype);
0134 end
0135
0136 hpaData.tissues(J)=[];
0137 if isfield(hpaData,'celltypes')
0138 hpaData.celltypes(J)=[];
0139 end
0140 if isfield(hpaData,'gene2Level')
0141 hpaData.gene2Level(:,J)=[];
0142 end
0143 if isfield(hpaData,'gene2Type')
0144 hpaData.gene2Type(:,J)=[];
0145 end
0146 if isfield(hpaData,'gene2Reliability')
0147 hpaData.gene2Reliability(:,J)=[];
0148 end
0149
0150
0151
0152 if ~isempty(hpaData.genes)
0153 I=~ismember(hpaData.genes,model.genes) | sum(hpaData.gene2Level,2)==0;
0154 else
0155 I=[];
0156 end
0157 hpaData.genes(I)=[];
0158 if isfield(hpaData,'gene2Level')
0159 hpaData.gene2Level(I,:)=[];
0160 end
0161 if isfield(hpaData,'gene2Type')
0162 hpaData.gene2Type(I,:)=[];
0163 end
0164 if isfield(hpaData,'gene2Reliability')
0165 hpaData.gene2Reliability(I,:)=[];
0166 end
0167
0168 I=strcmpi(arrayData.tissues,tissue);
0169
0170 if any(celltype)
0171 I=I & strcmpi(arrayData.celltypes,celltype);
0172 end
0173
0174
0175
0176 J=~ismember(arrayData.genes,model.genes) | myAll(isnan(arrayData.levels(:,I)),2);
0177 arrayData.genes(J)=[];
0178 arrayData.levels(J,:)=[];
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189 tempArrayLevels=arrayData.levels;
0190 tempArrayLevels(isnan(tempArrayLevels))=0;
0191 average=sum(tempArrayLevels,2)./sum(~isnan(arrayData.levels),2);
0192 if strcmpi(multipleCellScoring,'best')
0193 current=max(tempArrayLevels(:,I),[],2);
0194 else
0195 current=sum(tempArrayLevels(:,I),2)./sum(~isnan(arrayData.levels(:,I)),2);
0196 end
0197 if any(current)
0198 aScores=5*log(current./average);
0199 else
0200 aScores=[];
0201 end
0202 aScores(aScores>0)=min(aScores(aScores>0),10);
0203 aScores(aScores<0)=max(aScores(aScores<0),-5);
0204
0205
0206 [I J]=ismember(upper(hpaData.levels),upper(hpaLevelScores.names));
0207 if ~all(I)
0208 throw(MException('','There are expression level categories that do not match to hpaLevelScores'));
0209 end
0210 [K L M]=find(hpaData.gene2Level);
0211 scores=hpaLevelScores.scores(J);
0212 if strcmpi(multipleCellScoring,'best')
0213 hScores=max(sparse(K,L,scores(M),numel(hpaData.genes),numel(hpaData.tissues)),[],2);
0214 else
0215 hScores=mean(sparse(K,L,scores(M),numel(hpaData.genes),numel(hpaData.tissues)),2);
0216 end
0217
0218
0219 geneScores=inf(numel(model.genes),1)*-1;
0220 hpaScores=geneScores;
0221 arrayScores=geneScores;
0222
0223 [I J]=ismember(model.genes,hpaData.genes);
0224 hpaScores(I)=hScores(J(I));
0225 geneScores(I)=hScores(J(I));
0226 [I J]=ismember(model.genes,arrayData.genes);
0227 arrayScores(I)=aScores(J(I));
0228 geneScores(I & myIsInf(geneScores))=aScores(J(I & myIsInf(geneScores)));
0229
0230
0231 I=ismember(model.genes,hpaData.genes) | ismember(model.genes,arrayData.genes);
0232 model.genes(~I)=[];
0233 model.rxnGeneMat(:,~I)=[];
0234
0235
0236 [hpaExist hpaMap]=ismember(model.genes,hpaData.genes);
0237 [arrayExist arrayMap]=ismember(model.genes,arrayData.genes);
0238
0239
0240 rxnScores=ones(numel(model.rxns),1)*noGeneScore;
0241
0242
0243 for i=1:numel(model.rxns)
0244
0245 I=find(model.rxnGeneMat(i,:));
0246 if any(I)
0247
0248 if any(hpaExist(I))
0249
0250 if strcmpi(multipleGeneScoring,'best')
0251 rxnScores(i)=max(hScores(hpaMap(I(hpaExist(I)))));
0252 else
0253 rxnScores(i)=mean(hScores(hpaMap(I(hpaExist(I)))));
0254 end
0255 else
0256
0257 if any(arrayExist(I))
0258
0259 if strcmpi(multipleGeneScoring,'best')
0260 rxnScores(i)=max(aScores(arrayMap(I(arrayExist(I)))));
0261 else
0262 rxnScores(i)=mean(aScores(arrayMap(I(arrayExist(I)))));
0263 end
0264 end
0265 end
0266 end
0267 end
0268 end
0269
0270
0271
0272 function y=myIsInf(x)
0273 y=isinf(x);
0274 if isempty(y)
0275 y=[];
0276 end
0277 end
0278 function y=myAll(x,dim)
0279 y=all(x,dim);
0280 if isempty(y)
0281 y=[];
0282 end
0283 end