function [translatedModel] = translateDraftReconstruction(model)
% This function translates reaction and metabolite IDs in KBase draft 
% reconstructions to VMH nomenclature so the test suite can be performed on
% them. It does not make changes to the reaction content otherwise.
%
% USAGE: [translatedModel] = translateDraftReconstruction(model)
%
% INPUT
% model                     COBRA model structure generated by KBase 
%                           reconstruction pipeline
% OUTPUT
% translatedModel:          COBRA model structure with metabolites and
%                           reactions translated to VMH identifiers (but 
%                           not ran through refinement pipeline)
%
% AUTHOR
%       - Almut Heinken 09/2020

tol=0.0000001;

% Load metabolite and reaction translation tables
database=loadVMHDatabase;

translateRxns = readtable('ReactionTranslationTable.txt', 'Delimiter', '\t');
translateRxns = [translateRxns.Properties.VariableDescriptions;table2cell(translateRxns)];
translateRxns(find(~strncmp(translateRxns(:,1),'EX_',3)),:)=[];
translateMets = readtable('MetaboliteTranslationTable.txt', 'Delimiter', '\t');
translateMets = [translateMets.Properties.VariableDescriptions;table2cell(translateMets)];

model.mets=strrep(model.mets,'[c0]','[c]');
model.mets=strrep(model.mets,'[e0]','[e]');
model.mets=strrep(model.mets,'_c0','[c]');
model.mets=strrep(model.mets,'_e0','[e]');

for k=1:length(translateMets)
    model.mets = strrep(model.mets,translateMets{k,1}, translateMets{k,2});
end

% special case: two biomass reactions in model (bug in KBase?)
if length(intersect(model.rxns,{'bio1','bio2'}))==2
    model=removeRxns(model,'bio2');
end

% add some reactions that need to be in every reconstruction
essentialRxns={'DM_atp_c_','sink_PGPm1[c]'};

for i=1:length(essentialRxns)
    model = addReaction(model, essentialRxns{i}, 'reactionFormula', database.reactions{find(ismember(database.reactions(:, 1), essentialRxns{i})), 3});
end

% translate while keeping track of original reaction
rxns = strrep(model.rxns, '_c0', '');
rxns = strrep(rxns, '_e0', '');
%     rxns = strrep(rxns, 'R_', '');

for k=1:length(translateRxns)
    rxns = strrep(rxns,translateRxns{k,1}, translateRxns{k,2});
end

% rebuild model to get reactions in order-will get "reactions not
% sorted" error otherwise
[rxnsSorted,oldInd]=sort(rxns);

newModel=struct;
newModel.rxns={};
newModel.rxnNames={};
newModel.genes={};
newModel.grRules={};
newModel.rules={};
newModel.subSystems={};
newModel.lb=zeros(1);
newModel.ub=zeros(1);
newModel.S=[];
newModel.mets={};
newModel.metFormulas={};
newModel.metCharges=[];

for k=1:length(rxnsSorted)
    form=printRxnFormula(model,model.rxns{oldInd(k),1});
    name=model.rxnNames{oldInd(k)};
    lb=model.lb(oldInd(k));
    ub=model.ub(oldInd(k));
    geneRule=model.grRules{oldInd(k)};
    newModel=addReaction(newModel,rxnsSorted{k},'reactionName',name,'reactionFormula',...
        form{1},'lowerBound',lb,'upperBound',ub,'geneRule',geneRule,'subSystem','');
end
model=newModel;

% add metabolite formulas and charges
for k=1:length(model.mets)
    met=strrep(model.mets{k},'[c]','');
    met=strrep(met,'[e]','');
    if ~isempty(find(strcmp(database.metabolites(:,1),met)))
        model.metFormulas{k,1}=database.metabolites{find(strcmp(database.metabolites(:,1),met)),4};
        model.metCharges(k,1)=str2double(database.metabolites{find(strcmp(database.metabolites(:,1),met)),5});
    else
        model.metFormulas{k,1}='X';
        model.metCharges(k,1)=0;
    end
end
model=convertOldStyleModel(model);
model=changeObjective(model,model.rxns{find(strncmp(model.rxns,'bio',3))});

% Delete sink for petidoglycan if not needed
modelTest = removeRxns(model, 'sink_PGPm1[c]');
FBA = optimizeCbModel(modelTest, 'max');
if FBA.f>tol
    model= modelTest;
end

translatedModel=model;

end

