Bio Systems Biology Gene Essentiality

Version Compatibility

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.

Gene Essentiality Analysis

"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)

The governing principle: essentiality is a model-and-medium prediction, not a gene property

An in-silico knockout answers "can THIS network still make biomass on THIS medium after this perturbation?" Everything else follows from taking that literally:

• A gene knockout is NOT a reaction knockout. Genes map to reactions through the gene-protein-reaction (GPR) rule: AND = protein complex (all subunits needed), OR = isozymes (any one suffices). A gene in an OR with a viable partner does nothing when deleted - isozyme masking is the dominant source of false non-essential calls and the reason synthetic lethals exist. Always delete GENES and let the GPR decide which reactions close; deleting reactions directly is a different, usually wrong, analysis.
• The essential/non-essential CUTOFF is a policy choice, not a library default. COBRApy returns raw knockout growth rates; there is no 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.
• Essentiality is medium-dependent. Rich vs minimal media give different essential sets (a gene for a biosynthetic pathway is essential on minimal medium but dispensable when the product is supplied). To compare to an experiment, set the SAME medium the experiment used (LB vs M9), or the comparison is meaningless.
• Predictions have a ceiling (~85-93% accuracy on E. coli metabolic genes). False essentials come from a missing bypass/isozyme in the model; false non-essentials come from biology FBA cannot see (regulation, toxicity, essential non-metabolic or structural roles). Report MCC, not accuracy - essential genes are the minority class, so accuracy is inflated by the true-negative pile.

Decision: which knockout method, which order

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.

Single-Gene Deletion Screen

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')

Classify with a threshold sweep (report low-confidence calls)

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

MOMA / ROOM: the immediate, non-re-optimized mutant

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])