From 5ab9330925c9f4637e6fcb4c74c6c88e24b35192 Mon Sep 17 00:00:00 2001 From: Eduard Kerkhoven Date: Fri, 30 Oct 2026 21:28:16 +0100 Subject: [PATCH] loadProtData refuses a replicate count that does not match the file; calculateFfactor averages over the measured samples --- src/io/loadProtData.m | 13 ++++++++++++- src/limit_proteins/calculateFfactor.m | 5 +++-- test/VerificationMatrix.tsv | 2 +- test/unit_tests/protDataTests.m | 24 ++++++++++++++++++++++++ 4 files changed, 40 insertions(+), 4 deletions(-) diff --git a/src/io/loadProtData.m b/src/io/loadProtData.m index e4e3be27..53f10013 100644 --- a/src/io/loadProtData.m +++ b/src/io/loadProtData.m @@ -11,7 +11,9 @@ % replPerCond : double % vector with number of replicates for each condition in the dataset. % Example: [3, 2] if first condition has triplicates and second -% condition has duplicates. +% condition has duplicates. The sum must equal the number of replicate +% columns of the file (all columns after the UniProt ID), otherwise an +% error is raised. % % Name-Value Arguments % -------------------- @@ -109,6 +111,15 @@ if ~isstruct(protDataFile) checkFileExistence(protDataFile); fID = fopen(protDataFile); + nHeaderCols = numel(strsplit(fgetl(fID),'\t','CollapseDelimiters',false)); + if nHeaderCols ~= sum(replPerCond)+1 + fclose(fID); + error('loadProtData:replicateMismatch', ... + ['The header of %s has %d columns but replPerCond sums to %d, ' ... + 'so %d columns are expected (the UniProt ID and one per replicate).'], ... + protDataFile, nHeaderCols, sum(replPerCond), sum(replPerCond)+1) + end + frewind(fID); protDataRaw = textscan(fID,format,'Delimiter','\t','HeaderLines',1,'TreatAsEmpty',{'NA','na','NaN','#VALUE!'}); uniprotIDs = protDataRaw{1}; abundances = cell2mat(protDataRaw(2:end)); diff --git a/src/limit_proteins/calculateFfactor.m b/src/limit_proteins/calculateFfactor.m index 7cb42764..16796059 100644 --- a/src/limit_proteins/calculateFfactor.m +++ b/src/limit_proteins/calculateFfactor.m @@ -13,7 +13,8 @@ % -------------------- % protData : struct % structure with proteome data, from loadProtData (by default it instead -% attempts to load data/paxDB.tsv). +% attempts to load data/paxDB.tsv). With several samples, the abundance +% of a protein is averaged over the samples in which it is measured. % enzymes : cell % list of enzymes (default model.ec.enzymes). % modelAdapter : ModelAdapter @@ -99,7 +100,7 @@ protData.abundances = protData.level .* protData.MW; end -avgAbundances = mean(protData.abundances,2); +avgAbundances = mean(protData.abundances,2,'omitnan'); totalProt = sum(avgAbundances,'omitnan'); % Get enzymes in model diff --git a/test/VerificationMatrix.tsv b/test/VerificationMatrix.tsv index 0ec808fd..bf792183 100644 --- a/test/VerificationMatrix.tsv +++ b/test/VerificationMatrix.tsv @@ -33,7 +33,7 @@ gather_kcats/kcatsFromHomologs.m tc0095, tc0096 gather_kcats/extractEnzymeSubstratePairs.m tc0096, tc0015 gather_kcats/convertCustomKcats.m tc0091, tc0092, tc0094 change_model/upgradeEcModel.m tc0093, tc0094, tc0099 -io/loadProtData.m tc0101, tc0102 +io/loadProtData.m tc0101, tc0102, tc0260 change_model/addNewRxnsToEC.m tc0211, tc0232, tc0237 change_model/mapRxnsToConv.m tc0213 flux_analysis/ecFVA.m tc0214 diff --git a/test/unit_tests/protDataTests.m b/test/unit_tests/protDataTests.m index 70d1a6bb..2254c766 100644 --- a/test/unit_tests/protDataTests.m +++ b/test/unit_tests/protDataTests.m @@ -31,6 +31,30 @@ function testLoadProtDataCutLowestIsUsed_tc0102(testCase) verifyEqual(testCase, min(p20.abundances), 9) end +function testLoadProtDataReplicateCountMustMatchColumns_tc0260(testCase) + % replPerCond has to add up to the replicate columns of the file; columns + % beyond the sum are not read silently. + adapter = getTestAdapter(); + file = writeProteomics({'X1','X2'}, [1 2 3; 4 5 6]); + verifyError(testCase, @() loadProtData(2, 'protDataFile', file, 'modelAdapter', adapter), ... + 'loadProtData:replicateMismatch') + verifyError(testCase, @() loadProtData(4, 'protDataFile', file, 'modelAdapter', adapter), ... + 'loadProtData:replicateMismatch') + protData = loadProtData(3, 'protDataFile', file, 'modelAdapter', adapter, ... + 'addStdevs', 0, 'cutLowest', 0); + verifyEqual(testCase, protData.uniprotIDs, {'X1';'X2'}) +end + +function testCalculateFfactorAveragesOverMeasuredSamples_tc0261(testCase) + % A protein measured in one of two samples counts with that measurement + % in both sums of the f factor instead of dropping out of them. + adapter = getTestAdapter(); + protData.uniprotIDs = {'A';'B';'C'}; + protData.abundances = [10 NaN; 20 20; 30 30]; + f = calculateFfactor([], 'protData', protData, 'enzymes', {'A';'B'}, 'modelAdapter', adapter); + verifyEqual(testCase, f, 0.5, 'AbsTol', 1e-12) +end + function testApplyComplexDataStoichiometryByProteinId_tc0103(testCase) % The complex lists its proteins in another order than model.ec.enzymes; % each protein receives its own stoichiometry.