6. Simulating Deletions

[1]:
import pandas
from time import time

from cobra.io import load_model
from cobra.flux_analysis import (
    single_gene_deletion, single_reaction_deletion, double_gene_deletion,
    double_reaction_deletion)

cobra_model = load_model("textbook")
ecoli_model = load_model("iJO1366")

6.1. Knocking out single genes and reactions

A commonly asked question when analyzing metabolic models is what will happen if a certain reaction was not allowed to have any flux at all. This can tested using cobrapy by

[2]:
print('complete model: ', cobra_model.optimize())
with cobra_model:
    cobra_model.reactions.PFK.knock_out()
    print('pfk knocked out: ', cobra_model.optimize())
complete model:  <Solution 0.874 at 0x7fd813bf3390>
pfk knocked out:  <Solution 0.704 at 0x7fd813bf31d0>

For evaluating genetic manipulation strategies, it is more interesting to examine what happens if given genes are knocked out as doing so can affect no reactions in case of redundancy, or more reactions if gene when is participating in more than one reaction.

[3]:
print('complete model: ', cobra_model.optimize())
with cobra_model:
    cobra_model.genes.b1723.knock_out()
    print('pfkA knocked out: ', cobra_model.optimize())
    cobra_model.genes.b3916.knock_out()
    print('pfkB knocked out: ', cobra_model.optimize())
complete model:  <Solution 0.874 at 0x7fd813bf3c18>
pfkA knocked out:  <Solution 0.874 at 0x7fd813bf35c0>
pfkB knocked out:  <Solution 0.704 at 0x7fd813bf3588>

6.2. Single Deletions

Perform all single gene deletions on a model

[4]:
deletion_results = single_gene_deletion(cobra_model)

These can also be done for only a subset of genes

[5]:
single_gene_deletion(cobra_model, cobra_model.genes[:20])
[5]:
ids growth status
0 {b1849} 0.873922 optimal
1 {b1478} 0.873922 optimal
2 {b1276} 0.873922 optimal
3 {b0351} 0.873922 optimal
4 {b3733} 0.374230 optimal
5 {b1241} 0.873922 optimal
6 {b0116} 0.782351 optimal
7 {b0727} 0.858307 optimal
8 {b0356} 0.873922 optimal
9 {b2587} 0.873922 optimal
10 {b0118} 0.873922 optimal
11 {b0726} 0.858307 optimal
12 {b3736} 0.374230 optimal
13 {b3735} 0.374230 optimal
14 {b2296} 0.873922 optimal
15 {b3115} 0.873922 optimal
16 {b3732} 0.374230 optimal
17 {b0474} 0.873922 optimal
18 {s0001} 0.211141 optimal
19 {b3734} 0.374230 optimal

This can also be done for reactions

[6]:
single_reaction_deletion(cobra_model, cobra_model.reactions[:20])
[6]:
ids growth status
0 {AKGt2r} 8.739215e-01 optimal
1 {D_LACt2} 8.739215e-01 optimal
2 {Biomass_Ecoli_core} 0.000000e+00 optimal
3 {ETOHt2r} 8.739215e-01 optimal
4 {ADK1} 8.739215e-01 optimal
5 {ACONTb} 3.279963e-17 optimal
6 {ACt2r} 8.739215e-01 optimal
7 {ACONTa} 3.344590e-15 optimal
8 {ACKr} 8.739215e-01 optimal
9 {ALCD2x} 8.739215e-01 optimal
10 {ATPS4r} 3.742299e-01 optimal
11 {ACALD} 8.739215e-01 optimal
12 {ENO} 1.357454e-16 optimal
13 {EX_ac_e} 8.739215e-01 optimal
14 {ACALDt} 8.739215e-01 optimal
15 {CS} 4.757918e-15 optimal
16 {AKGDH} 8.583074e-01 optimal
17 {CO2t} 4.616696e-01 optimal
18 {CYTBD} 2.116629e-01 optimal
19 {ATPM} 9.166475e-01 optimal

6.3. Double Deletions

Double deletions run in a similar way.

[7]:
double_gene_deletion(
    cobra_model, cobra_model.genes[-5:]).round(4)
[7]:
ids growth status
0 {b2465, b3919} 0.7040 optimal
1 {b2935, b2464} 0.8739 optimal
2 {b0008, b2465} 0.8739 optimal
3 {b2465} 0.8739 optimal
4 {b0008} 0.8739 optimal
5 {b0008, b2464} 0.8648 optimal
6 {b2465, b2935} -0.0000 optimal
7 {b0008, b3919} 0.7040 optimal
8 {b2465, b2464} 0.8739 optimal
9 {b0008, b2935} 0.8739 optimal
10 {b3919} 0.7040 optimal
11 {b2464} 0.8739 optimal
12 {b2935} 0.8739 optimal
13 {b2935, b3919} 0.7040 optimal
14 {b3919, b2464} 0.7040 optimal

By default, the double deletion function will automatically use multiprocessing, splitting the task over up to 4 cores if they are available. The number of cores can be manually specified as well. Setting use of a single core will disable use of the multiprocessing library, which often aids debugging.

[8]:
start = time()  # start timer()
double_gene_deletion(
    ecoli_model, ecoli_model.genes[:25], processes=2)
t1 = time() - start
print("Double gene deletions for 200 genes completed in "
      "%.2f sec with 2 cores" % t1)

start = time()  # start timer()
double_gene_deletion(
    ecoli_model, ecoli_model.genes[:25], processes=1)
t2 = time() - start
print("Double gene deletions for 200 genes completed in "
      "%.2f sec with 1 core" % t2)

print("Speedup of %.2fx" % (t2 / t1))
Double gene deletions for 200 genes completed in 1.29 sec with 2 cores
Double gene deletions for 200 genes completed in 1.74 sec with 1 core
Speedup of 1.35x

Double deletions can also be run for reactions.

[9]:
double_reaction_deletion(
    cobra_model, cobra_model.reactions[2:7]).round(4)
[9]:
ids growth status
0 {ACONTa} 0.0000 optimal
1 {ACt2r, ACONTa} 0.0000 optimal
2 {ACKr} 0.8739 optimal
3 {ACt2r, ACKr} 0.8739 optimal
4 {ADK1, ACONTb} 0.0000 optimal
5 {ACONTa, ACONTb} 0.0000 optimal
6 {ACKr, ACONTa} 0.0000 optimal
7 {ACKr, ACONTb} 0.0000 optimal
8 {ACt2r, ADK1} 0.8739 optimal
9 {ACt2r, ACONTb} 0.0000 optimal
10 {ADK1, ACKr} 0.8739 optimal
11 {ACONTb} 0.0000 optimal
12 {ADK1, ACONTa} 0.0000 optimal
13 {ACt2r} 0.8739 optimal
14 {ADK1} 0.8739 optimal

6.4. Accessing individual deletion results

Note that the indices for deletions are python set objects. This is the appropriate type since the order of deletions does not matter. Deleting reaction 1 and reaction 2 will have the same effect as deleting reaction 2 and reaction 1.

To make it easier to access results all DataFrames returned by COBRAPpy deletion functions have a knockout indexer that makes that a bit simpler. Each entry in the indexer is treated as a single deletion entry. So you need to pass sets for double deletions.

[11]:
single = single_reaction_deletion(cobra_model)
double = double_reaction_deletion(cobra_model)

print(single.knockout["ATPM"])
print(double.knockout[{"ATPM", "TKT1"}])
       ids    growth   status
12  {ATPM}  0.916647  optimal
              ids   growth   status
859  {ATPM, TKT1}  0.90584  optimal

This can be used to get several deletions at once and will also work for Reaction or Gene objects (depending on what you deleted) directly.

[12]:
atpm = cobra_model.reactions.ATPM
tkt1 = cobra_model.reactions.TKT1
pfk = cobra_model.reactions.PFK

print(single.knockout[atpm, tkt1, pfk])
print(double.knockout[{atpm, tkt1}, {atpm, pfk}, {atpm}])

       ids    growth   status
12  {ATPM}  0.916647  optimal
44  {TKT1}  0.864759  optimal
71   {PFK}  0.704037  optimal
               ids    growth   status
425    {ATPM, PFK}  0.704037  optimal
859   {ATPM, TKT1}  0.905840  optimal
4125        {ATPM}  0.916647  optimal
[ ]: