options

methods = {'carveme','gapseq','kbase','coral_proj'};

outDir = fullfile('data/gap-filling/iterative');

% medium that has been used for gap filling
load(mediumFile)

% gap-filling database which contains permeability information
load(dbFile)

% output directory for figures
figOutDir = 'figures/exchanged_metabolites';

% specify the number of comparisons
n_comp = 8;
exchanged_met = cell(n_comp, 1);

for i=1:numel(habitat)
    for l = 1:numel(methods)
        % load model workspace
        modelWorkspace = fullfile(outDir, habitat{i},[methods{l}, '.mat']);
        load(modelWorkspace, 'GF')
        n = numel(GF);

        %% extract exported metabolites from sink reactions that have been added
        % during iterative gap filling
        exported_per_model = cell(n, 1);
        imported_per_model = cell(n, 1);

        for j=1:n
    
            % find exported metabolites sink reactions
            exported = regexp(GF{j}.rxns(contains(GF{j}.rxns, 'sink_')),...
                'MNXM\d+\[e]', 'match');
            exported = [exported{:}]';
            exported_per_model{j} = setdiff(exported, medium); clear exported
    
            % find imported metabolites from exchange reactions
            imported = regexp(GF{j}.rxns(contains(GF{j}.rxns, 'EX_')),...
                'MNXM\d+\[e]', 'match');
            imported = [imported{:}]';
            imported_per_model{j} = setdiff(imported, medium); clear imported
    
        end

        % get a list of exchange metabolites respectively
        exchanged_IDs = intersect(vertcat(imported_per_model{:}),vertcat(exported_per_model{:}));
        
        if i == 1
            exchanged_met{l,1} = exchanged_IDs;
        else
            exchanged_met{4+l,1} = exchanged_IDs;
        end
    end
end

% calculate the Jaccard index between each pair of exchanged
% metabolites for every approach

% loop over models
for j=1:size(exchanged_met, 2)
    tmp = double(size(exchanged_met, 1));
    % loop over methods
    for k=1:size(exchanged_met, 1)
        % loop over methods
        for l=1:size(exchanged_met, 1)
            tmp(k, l) = numel(intersect(vertcat(exchanged_met{k, j}),...
            vertcat(exchanged_met{l, j}))) / ...
            numel(union(vertcat(exchanged_met{k, j}),...
            vertcat(exchanged_met{l, j})));
        end
    end
end

row_names = {'c_carveme', 'c_gapseq', 'c_kbase','c_consensus','s_carveme', 's_gapseq', 's_kbase','s_consensus'};
writetable(array2table(tmp,...
        'VariableNames', row_names, 'RowNames', row_names),...
        fullfile(figOutDir, 'dist_methods.txt'),...
        'WriteVariableNames', true, 'WriteRowNames', true, 'Delimiter', '\t');
