Unverified Commit 8ca41feb authored by Laurent Heirendt's avatar Laurent Heirendt Committed by GitHub
Browse files

Merge pull request #151 from LCSB-BioCore/mva-knockout-2

Mva knockout 2
parents feea0fce e56c8a01
......@@ -7,6 +7,7 @@ id :: String
name :: Union{String, Nothing}
notes :: Dict{String, Vector{String}}
annotation :: Dict{String, Union{Vector{String}, String}}
associated_reactions :: Set{String}
````
"""
mutable struct Gene
......@@ -14,7 +15,13 @@ mutable struct Gene
name::Maybe{String}
notes::Notes
annotations::Annotations # everything is a String[]
associated_reactions::Set{String}
Gene(id::String = ""; name = nothing, notes = Notes(), annotations = Annotations()) =
new(id, name, notes, annotations)
Gene(
id::String = "";
name = nothing,
notes = Notes(),
annotations = Annotations(),
associated_reactions = Set{String}(),
) = new(id, name, notes, annotations, associated_reactions)
end
......@@ -65,6 +65,7 @@ Return the number of reactions contained in `model`.
"""
n_reactions(model::StandardModel)::Int = length(model.reactions)
"""
metabolites(model::StandardModel)
......@@ -253,6 +254,15 @@ function gene_annotations(model::StandardModel, id::String)::Maybe{Annotations}
model.genes[id].annotations
end
"""
gene_associated_reactions(model::StandardModel, id::String)::Set{String}
Return the reactions associated with gene `id` in `model`.
"""
function gene_associated_reactions(model::StandardModel, id::String)::Set{String}
model.genes[id].associated_reactions
end
"""
reaction_notes(model::StandardModel, id::String)::Notes
......
......@@ -105,3 +105,46 @@ function change_objective(
@objective(opt_model, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
end
end
"""
knockout(gene_ids::Array{String,1})
Callback function to set bounds of all reactions to zero which are affected by knocking out their respective genes
"""
function knockout(gene_ids::Vector{String})
# Dev note: the three nested for loops are inefficiency. However:
# - gene_ids (user input) will be probably only very few items
# - model.genes[gene_id].reactions are just a few reactions (most genes don't code for a lot of reactions)
# - reaction.grr also should only hold few items (reactions aren't coded by many different combinations of genes)
# Let's avoid premature optimization for now and see if anyone ever has problems with this
return (model, opt_model) -> begin
all_reactions = reactions(model)
for gene_id in gene_ids
for reaction_id in gene_associated_reactions(model, gene_id)
if all(
any(occursin.(gene_ids, gene_array)) > 0 for
gene_array in reaction_gene_association(model, reaction_id)
)
set_bound(
first(indexin([reaction_id], all_reactions)),
opt_model,
ub = 0,
lb = 0,
)
# Also set bounds for model object
set_bound(model, reaction_id, ub = 0, lb = 0)
end
end
end
end
end
"""
knockout(gene_id::String)
Callback function to set bounds of a reaction to zero which are affected by knocking out their respective genes
"""
function knockout(gene_id::String)
return knockout([gene_id])
end
function _map_reaction_to_genes!(model::StandardModel, rxn::Reaction)
if !isnothing(rxn.grr)
for gene_array in rxn.grr
for gene in gene_array
push!(model.genes[gene].associated_reactions, rxn.id)
end
end
end
end
"""
add!(model::StandardModel, rxns::Union{Vector{Reaction}, Reaction})
......@@ -11,6 +21,7 @@ end
function add!(model::StandardModel, rxn::Reaction)
model.reactions[rxn.id] = rxn
_map_reaction_to_genes!(model, rxn)
end
"""
......@@ -115,7 +126,14 @@ function rm!(::Type{Reaction}, model::StandardModel, ids::Vector{String})
end
function rm!(::Type{Reaction}, model::StandardModel, id::String)
delete!(model.reactions, id)
rxn = pop!(model.reactions, id)
if !isnothing(rxn.grr)
for gene_array in rxn.grr
for gene_id in gene_array
delete!(model.genes[gene_id].associated_reactions, id)
end
end
end
end
"""
......@@ -146,12 +164,43 @@ Remove all genes with `ids` from `model`.
rm!(Gene, model, ["g1", "g2"])
rm!(Gene, model, "g1")
"""
function rm!(::Type{Gene}, model::StandardModel, ids::Vector{String})
function rm!(
::Type{Gene},
model::StandardModel,
ids::Vector{String};
knockout::Bool = false,
)
for id in ids
rm!(Gene, model, id)
rm!(Gene, model, id, knockout = knockout)
end
end
function rm!(::Type{Gene}, model::StandardModel, id::String)
function rm!(::Type{Gene}, model::StandardModel, id::String; knockout::Bool = false)
if knockout
# Dev note: the three nested for loops are inefficiency. However:
# - gene_ids (user input) will be probably only very few items
# - model.genes[gene_id].reactions are just a few reactions (most genes don't code for a lot of reactions)
# - reaction.grr also should only hold few items (reactions aren't coded by many different combinations of genes)
# Let's avoid premature optimization for now and see if anyone ever has problems with this
for reaction_id in gene_associated_reactions(model, id)
reaction = get(model.reactions, reaction_id, nothing)
if isnothing(reaction)
return nothing
end
# AND inside the gene_array, so destroy as soon as one is missing
reaction.grr = filter(x -> !any(occursin.(id, x)), reaction.grr)
# OR outside, so all have to be deleted for the reaction to be deleted
if length(reaction.grr) == 0
rm!(Reaction, model, reaction.id)
end
end
end
delete!(model.genes, id)
end
function set_bound(model::StandardModel, reaction_id::String; ub, lb)
reaction = model.reactions[reaction_id]
reaction.lb = lb
reaction.ub = ub
end
@testset "single_knockout" begin
optimizer = Tulip.Optimizer
m = StandardModel()
add!(m, Metabolite("A"))
add!(m, Metabolite("B"))
add!(m, Gene("g1"))
add!(m, Gene("g2"))
add!(m, Reaction("v1", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1"]]))
add!(
m,
Reaction("v2", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1", "g2"]]),
)
add!(
m,
Reaction("v3", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1"], ["g2"]]),
)
add!(
m,
Reaction(
"v4",
metabolites = Dict("A" => -1.0, "B" => 1.0),
grr = [["g1", "g2"], ["g2"]],
),
)
opt_model = make_optimization_model(m, optimizer)
knockout("g1")(m, opt_model)
# Knockout should remove v1
@test normalized_rhs(opt_model[:lbs][1]) == 0
@test normalized_rhs(opt_model[:ubs][1]) == 0
# Knockout should remove [g1, g2] (AND) and thus remove reaction
@test normalized_rhs(opt_model[:lbs][2]) == 0
@test normalized_rhs(opt_model[:ubs][2]) == 0
# Knockout should remove [g1], but keep reaction (OR)
@test normalized_rhs(opt_model[:lbs][3]) == 1000
@test normalized_rhs(opt_model[:ubs][3]) == 1000
# Knockout should remove [g1, g2] (AND), but keep reaction (OR)
@test normalized_rhs(opt_model[:lbs][4]) == 1000
@test normalized_rhs(opt_model[:ubs][4]) == 1000
end
@testset "multiple_knockouts" begin
optimizer = Tulip.Optimizer
m = StandardModel()
add!(m, Metabolite("A"))
add!(m, Metabolite("B"))
add!(m, Gene("g1"))
add!(m, Gene("g2"))
add!(m, Gene("g3"))
add!(
m,
Reaction("v1", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1"], ["g3"]]),
)
add!(
m,
Reaction(
"v2",
metabolites = Dict("A" => -1.0, "B" => 1.0),
grr = [["g1", "g2"], ["g3"]],
),
)
add!(
m,
Reaction(
"v3",
metabolites = Dict("A" => -1.0, "B" => 1.0),
grr = [["g1"], ["g2"], ["g3"]],
),
)
opt_model = make_optimization_model(m, optimizer)
knockout(["g1", "g3"])(m, opt_model)
# Reaction 1 should be knocked out, because both
# gene1 and gene 3 are knocked out
@test normalized_rhs(opt_model[:lbs][1]) == 0
@test normalized_rhs(opt_model[:ubs][1]) == 0
# Reaction 2 should be knocked out, because both
# [g1, g2] is an AND relationship
@test normalized_rhs(opt_model[:lbs][1]) == 0
@test normalized_rhs(opt_model[:ubs][1]) == 0
# Reaction 3 should stay, because gene2 is still
# available (the arrays have an OR relationship)
@test normalized_rhs(opt_model[:lbs][3]) == 1000
@test normalized_rhs(opt_model[:ubs][3]) == 1000
end
@testset "single gene - single reaction" begin
m = StandardModel()
add!(m, Gene("g1"))
add!(m, Reaction("v1", metabolites = Dict{String,Int}(), grr = [["g1"]]))
@test m.genes["g1"].associated_reactions == Set(("v1",))
rm!(Reaction, m, "v1")
@test m.genes["g1"].associated_reactions == Set{String}()
end
@testset "single gene - multiple reactions" begin
m = StandardModel()
add!(m, Gene("g1"))
add!(m, Gene("g2"))
add!(m, Reaction("v1", metabolites = Dict{String,Int}(), grr = [["g1"]]))
add!(m, Reaction("v2", metabolites = Dict{String,Int}(), grr = [["g1"]]))
@test m.genes["g1"].associated_reactions == Set(("v1", "v2"))
rm!(Reaction, m, "v1")
@test m.genes["g1"].associated_reactions == Set(("v2",))
end
@testset "multiple genes - single reactions" begin
m = StandardModel()
add!(m, Gene("g1"))
add!(m, Gene("g2"))
add!(m, Reaction("v1", metabolites = Dict{String,Int}(), grr = [["g1"], ["g2"]]))
@test m.genes["g1"].associated_reactions == Set(("v1",))
@test m.genes["g2"].associated_reactions == Set(("v1",))
rm!(Reaction, m, "v1")
@test m.genes["g1"].associated_reactions == Set{String}()
@test m.genes["g1"].associated_reactions == Set{String}()
end
@testset "multiple genes - multiple reactions" begin
m = StandardModel()
add!(m, Gene("g1"))
add!(m, Gene("g2"))
add!(m, Reaction("v1", metabolites = Dict{String,Int}(), grr = [["g1"], ["g2"]]))
add!(m, Reaction("v2", metabolites = Dict{String,Int}(), grr = [["g1"], ["g2"]]))
@test m.genes["g1"].associated_reactions == Set(("v1", "v2"))
@test m.genes["g2"].associated_reactions == Set(("v1", "v2"))
rm!(Reaction, m, "v1")
@test m.genes["g1"].associated_reactions == Set(("v2",))
@test m.genes["g2"].associated_reactions == Set(("v2",))
end
"""
The gene-reaction rules (grr) are written as
[[gene1 and gene2...] or [gene3 and...] ...]
so for a reaction to be available it is sufficient that one group
is available, but inside a group all of the genes need to be available
"""
@testset "knockout_single_reaction" begin
m = StandardModel()
add!(m, Metabolite("A"))
add!(m, Metabolite("B"))
add!(m, Gene("g1"))
add!(m, Gene("g2"))
# Knockout should remove v1
add!(m, Reaction("v1", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1"]]))
# Knockout should remove [g1, g2] (AND) and thus remove reaction v2
add!(
m,
Reaction("v2", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1", "g2"]]),
)
# Knockout should remove [g1], but keep reaction v3 (OR)
add!(
m,
Reaction("v3", metabolites = Dict("A" => -1.0, "B" => 1.0), grr = [["g1"], ["g2"]]),
)
# Knockout should remove [g1, g2] (AND), but keep reaction v4 (OR)
add!(
m,
Reaction(
"v4",
metabolites = Dict("A" => -1.0, "B" => 1.0),
grr = [["g1", "g2"], ["g2"]],
),
)
rm!(Gene, m, "g1", knockout = true)
@test length(m.reactions) == 2
@test !haskey(m.reactions, "v1")
@test !haskey(m.reactions, "v2")
rxn = m.reactions["v3"]
@test length(rxn.grr) == 1
@test rxn.grr[1][1] == "g2"
rxn = m.reactions["v4"]
@test length(rxn.grr) == 1
@test rxn.grr[1][1] == "g2"
end
@testset "knockout_multiple_reactions" begin
m = StandardModel()
add!(m, Metabolite("A"))
add!(m, Metabolite("B"))
add!(m, Gene("g1"))
add!(m, Gene("g2"))
add!(m, Gene("g3"))
add!(
m,
Reaction(
"v1",
metabolites = Dict("A" => -1.0, "B" => 1.0),
grr = [["g1"], ["g2"], ["g3"]],
),
)
rm!(Gene, m, ["g1", "g3"], knockout = true)
rxn = m.reactions["v1"]
@test length(rxn.grr) == 1
@test rxn.grr[1][1] == "g2"
end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment