name: bio-systems-biology-gene-essentiality description: Performs in-silico single and double gene deletions, condition-dependent essentiality, and synthetic-lethality screens on genome-scale metabolic models with COBRApy, evaluating gene-protein-reaction rules and comparing FBA re-optimization against MOMA/ROOM minimal-adjustment. Use when predicting essential genes, finding synthetic-lethal pairs for drug targets, choosing a growth cutoff, deciding FBA vs MOMA vs ROOM for a knockout, making essentiality medium-specific to match an experiment, or validating predictions against Keio/Tn-seq/CRISPR screens with MCC. origin: openai4s category: bioskills/systems-biology metadata: tool_type: python primary_tool: cobrapy third_party: name: GPTomics/bioSkills repository: https://github.com/GPTomics/bioSkills commit: d91ed3d563019e649dc854c56ccd62551359488a license: MIT
Reference examples tested with: COBRApy 0.29+, Python 3.10+
Before using code patterns, verify installed versions match. If versions differ:
- Python: pip show <package> then help(module.function) to check signatures
If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
Note: single_gene_deletion / double_gene_deletion spawn worker processes, so call them inside if __name__ == '__main__': (or pass processes=1) or a spawn platform will recurse. delete_model_genes is deprecated since 0.25; use knock_out_model_genes / remove_genes.
"Which genes are essential for growth in my organism?" -> Delete each gene (via its GPR rule), re-solve the model, and call a gene essential when its knockout drops predicted growth below a chosen cutoff - always relative to a specific model and medium.
- Python: cobra.flux_analysis.single_gene_deletion(), double_gene_deletion(), moma(), room() (COBRApy)
An in-silico knockout answers "can THIS network still make biomass on THIS medium after this perturbation?" Everything else follows from taking that literally:
growth_cutoff argument. The cutoff (commonly <1-10% of wild type) moves the essential set in the "sick but alive" tail. It should be reported and swept (1/2/5/10%); genes whose call flips are low-confidence hypotheses.| Goal | Method | Assumption / when |
|---|---|---|
| Essential genes, evolved/adapted strain, or only hard lethality | FBA single_gene_deletion |
mutant re-optimizes to max growth; cheapest; lethality calls agree with MOMA anyway |
| Immediate/fresh transposon or CRISPR mutant (one growth cycle) | MOMA moma |
mutant stays closest in flux space to wild type (QP); fits fresh-mutant fluxes better |
| Fresh mutant where response is a few regulatory on/off switches | ROOM room |
minimizes the NUMBER of significantly changed fluxes (MILP); recovers short bypasses |
| Synthetic-lethal PAIRS | double_gene_deletion on viable singles |
both single KOs viable, double lethal; O(n^2), restrict the gene list |
| Higher-order lethal sets (triples/quads) | Fast-SL (flux-support pruning) | brute-force O(n^3+) infeasible; Fast-SL prunes by flux support |
| Condition/medium contrast | per-medium single_gene_deletion in with model: |
essentiality re-computed under each defined medium |
MOMA/ROOM classify lethality similarly to FBA; they differ mainly on the quantitative growth of sick-but-alive mutants. Match the method to the timescale of the actual experiment.
Goal: Rank every gene by the growth defect of its knockout and flag essential and growth-reducing genes.
Approach: Take wild-type growth once, then single_gene_deletion clamps each gene's reactions to zero through the GPR, re-optimizes, and returns a DataFrame with columns ids (a set holding the deleted gene id(s)), growth, and status. The caller applies the cutoff.
import cobra
from cobra.flux_analysis import single_gene_deletion
model = cobra.io.load_model('textbook')
wt_growth = model.slim_optimize()
results = single_gene_deletion(model) # DataFrame: ids (set of gene ids), growth, status
results['gene'] = results['ids'].apply(lambda s: list(s)[0]) # ids elements are gene-id STRINGS
results['relative'] = results['growth'] / wt_growth
ESSENTIAL_CUTOFF = 0.01 # KO grows < 1% of WT -> essential (policy, not a default)
essential = results[results['relative'] < ESSENTIAL_CUTOFF]
print(f'Essential genes: {len(essential)} / {len(model.genes)} on this medium')
def classify_essentiality(results, wt_growth, cutoffs=(0.01, 0.02, 0.05, 0.10)):
'''Classify genes and report how many calls flip across cutoffs (the sick-tail sensitivity).'''
rel = results['growth'] / wt_growth
calls = {c: set(results.loc[rel < c, 'gene']) for c in cutoffs}
core = set.intersection(*calls.values()) # essential at every cutoff -> high confidence
boundary = set.union(*calls.values()) - core # call depends on the cutoff -> low confidence
return core, boundary
from cobra.flux_analysis import moma, room
# FBA assumes the mutant re-optimizes; a fresh knockout has not re-wired its regulation yet.
# MOMA keeps mutant flux closest (Euclidean) to wild type; ROOM minimizes the count of changed
# fluxes. Both need a wild-type reference solution and a QP/MILP-capable solver.
from cobra.util.solver import linear_reaction_coefficients
biomass = list(linear_reaction_coefficients(model))[0] # the objective (biomass) reaction
wt = model.optimize()
with model:
model.genes.get_by_id('b2276').knock_out() # context-aware; reverts on block exit
moma_sol = moma(model, solution=wt, linear=True) # linear=True = fast LP approximation (lMOMA)
# moma_sol.objective_value is the MINIMIZED ADJUSTMENT, not growth; read the biomass flux.
print('MOMA mutant growth:', moma_sol.fluxes[biomass.id])
Goal: Find gene PAIRS that are viable singly but lethal together - redundant pathways and isozymes, and candidate combination drug targets.
Approach: Restrict to genes whose single knockout is viable (a synthetic lethal requires both singles viable), run pairwise double_gene_deletion, and keep pairs whose double-KO growth falls below the cutoff. Cost is O(n^2), so subset the gene list. Score interactions against the multiplicative neutral expectation (independent effects on an exponential growth rate).
from cobra.flux_analysis import double_gene_deletion
viable = list(results.loc[results['relative'] > ESSENTIAL_CUTOFF, 'gene'])[:60] # cap the O(n^2) sweep
dbl = double_gene_deletion(model, gene_list1=viable, gene_list2=viable)
dbl['n'] = dbl['ids'].apply(len)
sl_pairs = dbl[(dbl['n'] == 2) & (dbl['growth'] / wt_growth < ESSENTIAL_CUTOFF)]
print(f'Synthetic-lethal pairs: {len(sl_pairs)}')
# Genome-scale and higher-order (triple/quad) sets: use Fast-SL flux-support pruning (Pratapa 2015),
# not brute force.
Goal: Compare essential-gene sets across defined media to separate core-essential genes from condition-specific ones.
Approach: Apply each medium inside a with model: block (so it reverts), run the deletion screen, and take intersections/differences of the essential sets. Define media with real functions, not lambdas (a lambda cannot contain an assignment).
def aerobic(m):
m.reactions.EX_o2_e.lower_bound = -20
def anaerobic(m):
m.reactions.EX_o2_e.lower_bound = 0
def essential_set(model, setup):
with model:
setup(model)
wt = model.slim_optimize()
res = single_gene_deletion(model)
return set(res.loc[res['growth'] / wt < ESSENTIAL_CUTOFF, 'ids'].apply(lambda s: list(s)[0]))
sets = {name: essential_set(model, fn) for name, fn in [('aerobic', aerobic), ('anaerobic', anaerobic)]}
core = set.intersection(*sets.values())
condition_specific = {k: v - core for k, v in sets.items()}
# Compare predicted essentials to an experimental set (Keio single-KO, Tn-seq, or CRISPR fitness),
# on the SAME medium. Use MCC, not accuracy: essential genes are a minority class, so accuracy is
# inflated by the large true-negative pile.
from sklearn.metrics import matthews_corrcoef
def score(predicted_essential, experimental_essential, all_genes):
y_pred = [g in predicted_essential for g in all_genes]
y_true = [g in experimental_essential for g in all_genes]
return matthews_corrcoef(y_true, y_pred)
| Symptom | Cause | Fix |
|---|---|---|
| Central gene predicted non-essential | it has an isozyme (OR in the GPR) that stays open | expected; that gene is a synthetic-lethal candidate, not truly dispensable |
| Deleting a reaction gives different results than deleting its gene | reaction KO ignores GPR; a gene may map to several reactions or share them | delete GENES (single_gene_deletion / gene.knock_out()), not reactions |
TypeError: 'set' object ... .id |
ids column holds sets of gene-id STRINGS, not gene objects |
list(s)[0] gives the id string directly; no .id |
| Essential set disagrees with the paper | medium mismatch (LB vs M9) or a different cutoff | set the experiment's medium; report and sweep the cutoff |
| Script recurses / spawns endlessly | deletion functions parallelize; no __main__ guard |
wrap in if __name__ == '__main__': or pass processes=1 |
AttributeError: delete_model_genes |
deprecated since cobra 0.25 | use knock_out_model_genes / remove_genes, or gene.knock_out() |
| High accuracy but poor agreement on real essentials | accuracy inflated by true negatives (minority class) | report MCC and sensitivity, not accuracy |