equilibrator_tools.jl 4.16 KB
Newer Older
St. Elmo's avatar
St. Elmo committed
1
2
3
4
5
"""
buildrxnstring(rxn)

Get rxn in string format for Equilibrator.
"""
6
function build_rxn_string(rxn::Reaction, compoundtype="kegg")
St. Elmo's avatar
St. Elmo committed
7
8
    pos_s = []
    neg_s = []
9
10
11
12
13
14
15
16
17
18

    if compoundtype == "kegg"
        cid = "kegg.compound"
    else
        cid = "bigg.metabolite"
    end

    for (met, coeff) in rxn.metabolites
        metid = get(met.annotation, cid, [""])[1]
        metid == "" && continue 
St. Elmo's avatar
St. Elmo committed
19
        if coeff > 0.0
20
21
22
23
24
            if compoundtype == "kegg" 
                push!(pos_s, "$(coeff) KEGG:$metid")
            else
                push!(pos_s, "$(coeff) bigg.metabolite:$metid")
            end
St. Elmo's avatar
St. Elmo committed
25
        else
26
27
28
29
30
            if compoundtype == "kegg" 
                push!(neg_s, "$(abs(coeff)) KEGG:$metid")
            else
                push!(neg_s, "$(abs(coeff)) bigg.metabolite:$metid")
            end 
St. Elmo's avatar
St. Elmo committed
31
32
33
34
35
        end
    end
    return join(neg_s, " + ")*" = "*join(pos_s, " + ") # keep order for ease of use later
end

36

St. Elmo's avatar
St. Elmo committed
37
38
39
"""
gibbs_arr = mapGibbs(rxns; dgtype="zero", ph=7.0, ionic_str="100 mM")

40
Return an dict of rxn.id => ΔG of the specified dgtype.
St. Elmo's avatar
St. Elmo committed
41
"""
42
43
44
45
46
47
48
49
function map_gibbs_rxns(rxns::Array{Reaction, 1}; dgtype="zero", ph=7.0, ionic_str="100 mM", usekegg=true) 
    
    if usekegg
        rxns_strings = [build_rxn_string(rxn, "kegg") for rxn in rxns]
    else
        rxns_strings = [build_rxn_string(rxn, "bigg") for rxn in rxns]
    end
    
St. Elmo's avatar
St. Elmo committed
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
    if dgtype == "phys"
        bals, gs, errs = py"pygetdgprimephys"(rxns_strings, ph, ionic_str)
    elseif dgtype == "prime" 
        bals, gs, errs = py"pygetdgprime"(rxns_strings, ph, ionic_str)
    else # "zero"
        bals, gs, errs = py"pygetdg0"(rxns_strings, ph, ionic_str)
    end

    gibbs = Dict{String, Measurement{Float64}}()
    for (rxn, g, err) in zip(rxns, gs, errs)
        if err < 1000.0 # ignore crazy errors
            gibbs[rxn.id] = g ± err 
        else
            gibbs[rxn.id] = 0.0 ± 0.0    
        end
    end
    return gibbs
end

"""
70
map_gibbs_external(fluxres, gibbs)
St. Elmo's avatar
St. Elmo committed
71
72

Calculate the Gibbs free energy change taking only the external fluxes into account.
73
NB: you need to account for the biomass function separately.
St. Elmo's avatar
St. Elmo committed
74
75

Fluxres can be both a ReactionFluxes object or a Dict with rxnid -> flux.
St. Elmo's avatar
St. Elmo committed
76
"""
77
function map_gibbs_external(fluxres::Dict{String, Float64}, gibbs)
St. Elmo's avatar
St. Elmo committed
78
    total_ΔG = 0.0 ± 0.0
St. Elmo's avatar
St. Elmo committed
79
    missing_flux = 0.0
St. Elmo's avatar
St. Elmo committed
80
81
    for (rxnid, v) in fluxres
        if startswith(rxnid, "EX_")
St. Elmo's avatar
St. Elmo committed
82
            if gibbs[rxnid]  0.0
83
                missing_flux += abs(v)
St. Elmo's avatar
St. Elmo committed
84
85
            end    
            total_ΔG -= v * gibbs[rxnid] # negative here because "combustion" is actually Gibbs value not formation  
St. Elmo's avatar
St. Elmo committed
86
87
        end
    end
St. Elmo's avatar
St. Elmo committed
88
    return total_ΔG, missing_flux/sum(abs, values(fluxres)) # units J/gDW/h
89
90
91
end

"""
St. Elmo's avatar
St. Elmo committed
92
map_gibbs_internal(fluxres, gibbs, biomassid="BIOMASS")
93
94

Calculate the Gibbs free energy change taking only the internal fluxes into account.
St. Elmo's avatar
St. Elmo committed
95
NB: you need to account for the biomass function separately. NB: 
96

St. Elmo's avatar
St. Elmo committed
97
Fluxres are a Dict with rxnid -> flux.
98
"""
St. Elmo's avatar
St. Elmo committed
99
function map_gibbs_internal(fluxres::Dict{String, Float64}, gibbs, biomassid="BIOMASS")
100
    total_ΔG = 0.0 ± 0.0
St. Elmo's avatar
St. Elmo committed
101
    missing_flux = 0.0
St. Elmo's avatar
St. Elmo committed
102
    found_flux = 0.0
103
    for (rxnid, v) in fluxres
St. Elmo's avatar
St. Elmo committed
104
        if !startswith(rxnid, "EX_") && !contains(rxnid, biomassid) # ignore exchange reactions and biomass function
St. Elmo's avatar
St. Elmo committed
105
            if gibbs[rxnid]  0.0
106
                missing_flux += abs(v)
St. Elmo's avatar
St. Elmo committed
107
108
            else
                found_flux += abs(v)
St. Elmo's avatar
St. Elmo committed
109
110
            end 
            total_ΔG += v * gibbs[rxnid] # add because this is not formation but rather just adding equations (the flux direction sign compensates)
111
112
        end
    end
St. Elmo's avatar
St. Elmo committed
113
    return total_ΔG, missing_flux/(missing_flux+found_flux) # units J/gDW/h
114
115
116
end

"""
St. Elmo's avatar
St. Elmo committed
117
save_gibbs(path, gibbs)
118

St. Elmo's avatar
St. Elmo committed
119
Save gibbs dict. Saved as String => [mag, err]
120
"""
St. Elmo's avatar
St. Elmo committed
121
function save_gibbs(path, gibbs)
122
123
124
125
    decomp = Dict{String, Array{Float64,1}}()
    for (k, v) in gibbs
        decomp[k] = [Measurements.value(v), Measurements.uncertainty(v)]
    end
St. Elmo's avatar
St. Elmo committed
126
127
128
    open(path, "w") do io
        JSON.print(io, decomp)
    end
129
130
131
132
133
134
135
end

"""
load_Gibbs(path)

Load Gibbs dict. Loads String => [mag, err]
"""
St. Elmo's avatar
St. Elmo committed
136
137
function load_gibbs(path)
    decomp = JSON.parsefile(path)
138
139
140
141
142
    gibbs = Dict{String, Measurement{Float64}}()
    for (k, vs) in decomp
        gibbs[k] = vs[1] ± vs[2]
    end
    return gibbs
St. Elmo's avatar
St. Elmo committed
143
end