Skip to content

Commit

Permalink
Merge pull request #2354 from rosieluo2021/metaPD
Browse files Browse the repository at this point in the history 
update overlap maps to symmetric maps
  • Loading branch information
@rmtfleming
rmtfleming authored Nov 7, 2024
2 parents ebbe615 + b600a6b commit 823b6e9
Show file tree
Hide file tree
Showing 2 changed files with 158 additions and 198 deletions.
150 changes: 83 additions & 67 deletions src/dataIntegration/XomicsToModel/analysis/compareXomicsModels.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
% This function is used to compare generated models from xomicsToModel
% pipeline
%
%Author(s): - Xi Luo
%Author(s): - Xi Luo, update at 2024/10
% - Hanneke Leegwater (2022)

%% Check input params
Expand All @@ -36,89 +36,105 @@
error('printFlag should be a number or a bool')
end

%%
models=multiModels;
%% Initialize
models = multiModels;

if isstruct(models) && length(fieldnames(models))>1
fields=fieldnames(models);
statMets=cell(length(fields));
statRxns=cell(length(fields));
statGenes=cell(length(fields));
for i=1:length(fields)
for j=1:length(fields)
%mets
overlap_mets.(fields{i}).(fields{j})=models.(fields{i}).mets(ismember(models.(fields{i}).mets,models.(fields{j}).mets));
statMets{i,j}=sum(ismember(models.(fields{i}).mets,models.(fields{j}).mets));
%rxns
overlap_rxns.(fields{i}).(fields{j})=models.(fields{i}).rxns(ismember(models.(fields{i}).rxns,models.(fields{j}).rxns));
statRxns{i,j}=sum(ismember(models.(fields{i}).rxns,models.(fields{j}).rxns));
%genes
overlap_genes.(fields{i}).(fields{j})=models.(fields{i}).genes(ismember(models.(fields{i}).genes,models.(fields{j}).genes));
statGenes{i,j}=sum(ismember(models.(fields{i}).genes,models.(fields{j}).genes));
fields = fieldnames(models);
numModels = length(fields);

% Initialize matrices for statistics and results
statMets = cell(numModels);
statRxns = cell(numModels);
statGenes = cell(numModels);
overlap_mets = struct();
overlap_rxns = struct();
overlap_genes = struct();

% Loop over all model pairs to compute overlaps
for i = 1:numModels
for j = 1:numModels
% Mets overlap
overlap_mets.(fields{i}).(fields{j}) = models.(fields{i}).mets(ismember(models.(fields{i}).mets, models.(fields{j}).mets));
statMets{i,j} = length(overlap_mets.(fields{i}).(fields{j}));

% Rxns overlap
overlap_rxns.(fields{i}).(fields{j}) = models.(fields{i}).rxns(ismember(models.(fields{i}).rxns, models.(fields{j}).rxns));
statRxns{i,j} = length(overlap_rxns.(fields{i}).(fields{j}));

% Genes overlap
overlap_genes.(fields{i}).(fields{j}) = models.(fields{i}).genes(ismember(models.(fields{i}).genes, models.(fields{j}).genes));
statGenes{i,j} = length(overlap_genes.(fields{i}).(fields{j}));
end
end
colname=cell2table(fields);
statMets=cell2table(statMets,"VariableNames",fields);
statMets=[colname statMets];
statRxns=cell2table(statRxns,"VariableNames",fields);
statRxns=[colname statRxns];
statGenes=cell2table(statGenes,"VariableNames",fields);
statGenes=[colname statGenes];
%all overlapped
overlap_mets.alloverlap=models.(fields{1}).mets(ismember(models.(fields{1}).mets,models.(fields{2}).mets));
overlap_rxns.alloverlap=models.(fields{1}).rxns(ismember(models.(fields{1}).rxns,models.(fields{2}).rxns));
overlap_genes.alloverlap=models.(fields{1}).genes(ismember(models.(fields{1}).genes,models.(fields{2}).genes));
for i=1:length(fields)-2
overlap_mets.alloverlap=overlap_mets.alloverlap(ismember(overlap_mets.alloverlap,models.(fields{i+2}).mets));
overlap_rxns.alloverlap=overlap_rxns.alloverlap(ismember(overlap_rxns.alloverlap,models.(fields{i+2}).rxns));
overlap_genes.alloverlap=overlap_genes.alloverlap(ismember(overlap_genes.alloverlap,models.(fields{i+2}).genes));
end
overlapResults.mets=overlap_mets;
overlapResults.rxns=overlap_rxns;
overlapResults.genes=overlap_genes;
% Convert statistics to tables
statMets = cell2table(statMets, 'VariableNames', fields, 'RowNames', fields);
statRxns = cell2table(statRxns, 'VariableNames', fields, 'RowNames', fields);
statGenes = cell2table(statGenes, 'VariableNames', fields, 'RowNames', fields);

% Store results in output structure
overlapResults.mets = overlap_mets;
overlapResults.rxns = overlap_rxns;
overlapResults.genes = overlap_genes;
% Calculate all overlap across all models
overlapResults.mets.alloverlap = intersectMultipleModels(models, fields, 'mets');
overlapResults.rxns.alloverlap = intersectMultipleModels(models, fields, 'rxns');
overlapResults.genes.alloverlap = intersectMultipleModels(models, fields, 'genes');

statistic.overlapnumber_mets=statMets;
statistic.overlapnumber_rxns=statRxns;
statistic.overlapnumber_genes=statGenes;
%met proportion
metsdata=statistic.overlapnumber_mets{:,2:end};
[max_a,index]=max(metsdata);
xa=repmat(max_a',[1 length(statistic.overlapnumber_mets{:,1})]);
pro=round(metsdata./xa*100,2);
pro=num2cell(pro);
pro=cell2table(pro,"VariableNames",fields);
pro=[colname pro];
statistic.overlaporportion_mets=pro;
%rxns propotion
metsdata=statistic.overlapnumber_rxns{:,2:end};
[max_a,index]=max(metsdata);
xa=repmat(max_a',[1 length(statistic.overlapnumber_rxns{:,1})]);
pro=round(metsdata./xa*100,2);
pro=num2cell(pro);
pro=cell2table(pro,"VariableNames",fields);
pro=[colname pro];
statistic.overlaporportion_rxns=pro;
%gene propotion
metsdata=statistic.overlapnumber_genes{:,2:end};
[max_a,index]=max(metsdata);
xa=repmat(max_a',[1 length(statistic.overlapnumber_genes{:,1})]);
pro=round(metsdata./xa*100,2);
pro=num2cell(pro);
pro=cell2table(pro,"VariableNames",fields);
pro=[colname pro];
statistic.overlaporportion_genes=pro;

% Calculate symmetric proportion matrix (A ∩ B) / (A ∪ B)
statistic.overlapproportion_mets = calculateProportionMatrix(models, fields, 'mets', statMets);
statistic.overlapproportion_rxns = calculateProportionMatrix(models, fields, 'rxns', statRxns);
statistic.overlapproportion_genes = calculateProportionMatrix(models, fields, 'genes', statGenes);

else
disp('please check the input variable')
end


%% Print tables with output if printFlag = 1
% Print tables with output if printFlag = 1
if printFlag ==1
disp('Number of overlapping mets between models is:')
statistic.overlapnumber_mets

disp('Number of overlapping rxns between models is:')
statistic.overlapnumber_rxns

disp('Number of overlapping genes between models is:')
statistic.overlapnumber_genes
end

end

% Function to calculate the intersection across multiple models
function allOverlap = intersectMultipleModels(models, fields, featureType)
allOverlap = models.(fields{1}).(featureType);
for i = 2:length(fields)
allOverlap = intersect(allOverlap, models.(fields{i}).(featureType));
end
end

% Function to calculate proportion matrix (A ∩ B) / (A ∪ B)
function proportionTable = calculateProportionMatrix(models, fields, featureType, statTable)
numModels = length(fields);
proportionMatrix = cell(numModels);

for i = 1:numModels
for j = 1:numModels
% Calculate the union (A ∪ B) for proportion
unionCount = length(union(models.(fields{i}).(featureType), models.(fields{j}).(featureType)));
if unionCount > 0
proportionMatrix{i,j} = round(100 * statTable{i,j} / unionCount, 2); % (A ∩ B) / (A ∪ B)
else
proportionMatrix{i,j} = NaN; % Avoid division by zero
end
end
end

% Convert to table for easy reading
proportionMatrix = cell2table(proportionMatrix, 'VariableNames', fields, 'RowNames', fields);
proportionTable = proportionMatrix;
end
Loading

0 comments on commit 823b6e9

Please sign in to comment.