main.m~

close all force; clear variables; clc

%% Startup the COBRA Toolbox
initCobraToolbox(false) % Don't update the toolbox
changeCobraSolver ('gurobi', 'all', 1); % For large models
clc

%% Load data and model
addpath('functions')  
% Load omic data
abundance = readtable('data/sup_material_8.xlsx');
expresion = readtable('data/DESeq2_normalised_counts.xls');
% Get list of proteinECNumbers
proteinECNumbers_FileName = 'ProteinECNumbers.mat';
if isfile(proteinECNumbers_FileName)
    load(proteinECNumbers_FileName, 'ProteinECNumbers');
else
    ProteinECNumbers = getProteinCodes(abundance);
    save(proteinECNumbers_FileName, 'ProteinECNumbers');
end

Osorio_model=readCbModel('models/Astrocyte_Osorio2019.xml');
% Load Recon3D model
load('models/Recon3D_301.mat')
% Recon3D stats
N_filled = find(~cellfun(@isempty,Recon3D.rxnECNumbers));
percent_of_filled = length(N_filled)/length(Recon3D.rxnECNumbers) * 100;
fprintf('%.02f%% of the reactions in Recon3D are associated with EC numbers\n', percent_of_filled)
N_filled = find(~cellfun(@isempty,Recon3D.rules));
percent_of_filled = length(N_filled)/length(Recon3D.rules) * 100;
fprintf('%.02f%% of the reactions in Recon3D are associated with gene rules\n', percent_of_filled)
% Standarize Recon3D rxnECNumber rules
Recon3D = harmonizeReconECnumbers(Recon3D);
Recon3D = makePrules(Recon3D);

%% Map omic data to Recon3D reactions
disp('Mapping omic data to Recon3D reactions ...')

abundanceTable_FileName = 'out/abundanceTable.mat';
expressionTable_FileName = 'out/expressionTable.mat';
if isfile(expressionTable_FileName) && isfile(abundanceTable_FileName)
    load(abundanceTable_FileName, 'abundanceTable');
    load(expressionTable_FileName, 'expressionTable');
else
    profile on
    abundanceTable = mapAbundance(Recon3D, abundance, ProteinECNumbers);
    expresion = EnsemblToEntrez(expresion, 'data/mart_export.txt');
    expressionTable = mapExpression(Recon3D, expresion);
    profile viewer

    save(abundanceTable_FileName, 'abundanceTable');
    save(expressionTable_FileName, 'expressionTable');
end


%% Transcriptomic and proteomic data integration
disp('Integrating Transcriptomic and Proteomic data ...')
omicIntegratedData = omicIntegrationPCA(abundanceTable, expressionTable);
omicIntegratedData = log(abs(min(omicIntegratedData)) + omicIntegratedData);

%% Recon3D Model reduction to astrocite specific GEM
disp('Creating Astrocite Specific Model ...')
options.solver = 'iMAT';
options.expressionRxns = omicIntegratedData;
options.threshold_lb = 500;
options.threshold_ub = -500;

astrociteModel = createTissueSpecificModel(Recon3D, options);
astrociteModel.description = 'MendozaAstrociteModel';
save('out/dirtyAstrociteModel.mat', 'astrociteModel');

%% First FBA to check flux
FBAsolution = optimizeCbModel(astrociteModel,'max');
fprintf("The first FBA flux after reduction is %f\n", FBAsolution.f)

%% Cleaning the model
disp('Searching and removing Leaking metabolites ...')
LeakRxns = searchLeakMetabolites(astrociteModel);
if ~isempty(LeakRxns)
    disp('Removing leaking metabolites');
    modelNoLeak = removeRxns(astrociteModel, LeakRxns);
end
save('out/NoLeakAstrociteModel.mat', 'modelNoLeak');

% FBA after removing the leaking metabolites
FBAsolution = optimizeCbModel(modelNoLeak,'max');
fprintf("The FBA flux after removing the leaking metabolites is %f\n", FBAsolution.f)

%% Set DMEM medium
astrociteModelDMEM = setBoundriesDMEMedium(modelNoLeak);
FBAsolution = optimizeCbModel(astrociteModelDMEM ,'max');

% How many exchange reactions are active
modelExchanges = strmatch('EX_',astrociteModelDMEM.rxns);
% modelExchanges = findExchangeReactions(astrociteModelDMEM);
exRxns = find(ismember(astrociteModelDMEM.rxns,astrociteModelDMEM.rxns(modelExchanges)) & FBAsolution.v ~= 0);
astrociteModelDMEM.rxns(exRxns)

allExchanges = findExchangeReactions(astrociteModelDMEM, 1);
modelExchangesLogic = zeros(1, length(astrociteModelDMEM.rxns));
modelExchangesLogic(modelExchanges) = 1;
% printRxnFormula(astrociteModelDMEM, 'rxnAbbrList', astrociteModelDMEM.rxns(allExchanges & ~modelExchangesLogic));
% Close other exchange reactions
astrociteModelDMEM.lb(ismember(astrociteModelDMEM.rxns,astrociteModelDMEM.rxns(allExchanges == 1 & ~modelExchangesLogic)))=0;
FBAsolution = optimizeCbModel(astrociteModelDMEM ,'max');

only_recon3d = compareModelsRxns(Recon3D, astrociteModelDMEM);
only_astrcyte = compareModelsRxns(astrociteModelDMEM, Recon3D);
fprintf("reactions only recon: %d only astrocyte %d shared %d \n", ...
    length(only_recon3d), length(only_astrcyte), length(Recon3D.rxns) - length(only_recon3d))

profile on
only_osorio = compareModelsRxns(Osorio_model, astrociteModelDMEM);
only_astrcyte = compareModelsRxns(astrociteModelDMEM, Osorio_model);
fprintf("reactions only recon: %d only astrocyte %d shared %d %d\n", ...
    length(only_osorio), length(only_astrcyte), length(Osorio_model.rxns) - length(only_osorio), length(astrociteModelDMEM.rxns) - length(only_astrcyte))
profile viewer
%% Contexto
omicIntegratedData(omicIntegratedData < 0) = 0;

%% Graficas
figure;
subplot(2,3,1);
histogram(table2array(abundanceTable(:,:)))
ylim([0 14000])
title('Abundances in Rxns')
subplot(2,3,2);
histogram(table2array(expressionTable(:,:)))
ylim([0 14000])
title('Expressions in Rxns')
subplot(2,3,3);
histogram(omicIntegratedData)
ylim([0 14000])
title('Integrated data in Rxns')

subplot(2,3,4);
histogram(real(log(table2array(abundanceTable(:,:)))))
ylim([0 8000])
title('Log Abundances in Rxns')
subplot(2,3,5);
histogram(real(log(table2array(expressionTable(:,:)))))
ylim([0 8000])
title('Log Expressions in Rxns')
subplot(2,3,6);
histogram(real(log(omicIntegratedData)))
ylim([0 8000])
title('Log Integrated data in Rxns')