2細菌が同じ環境にいるときの代謝を推定する
前回の記事の続き。2種類の細菌が同じ環境にいるときに、Cooperative/Competitiveな関係になるのか、どのように代謝物を交換するのか、といった細菌間相互作用をCOBRApyを使ったFBA (Flux Balance Analysis)で推定してみる。
ここでは以下の論文のFig. 4dで示されているように、Bacteroides caccae ATCC 43185とLactobacillus rhamnosus GGをDMEM 6429という培地で共培養した場合の細菌間相互作用を推定してみる。
本記事の検証に使ったコードはgithubにアップした。
２生物種のGEMsを作る方法
この論文ではBacteroides caccae ATCC 43185とLactobacillus rhamnosus GGのGEMs (Genome-scale metabolic models) を連結させた代謝モデルを作り、両細菌の生育率の合計が最大になる場合の代謝流束をFBAで求めている。メソッドを読むと筆者らが2015年の論文に用いた細菌間に共通する区画*1を設ける手法を利用していた。この手法のオリジナルは以下の2010年の論文のようだ。このオリジナル論文のFig. 2が仕様を理解する上でわかりやすかった。
引用元：Kitgord and Segre. PLoS Comput Biol. 2010 Nov 18;6(11):e1001002. doi: 10.1371/journal.pcbi.1001002.
Fig. ABは各細菌のGEMsである。区画CYT1とCYT2は各細菌の細胞内領域、EX1とEX2は各細菌の細胞外領域、X・Y・Zは代謝物質を表す。この手法ではFig. Cのように２つのGEMsに追加して、細菌間で共通の区画ENVを設定する。
これはStoichiometric matrixで表すとFig. EFGのように表せる。各細菌のStoichiometric matrix (Fig. EF) に若干手を加えるだけで２種共培養のStoichiometric matrix (Fig. G) を作ることができる。
MATLAB版のCOBRAには2細菌のモデルを作るための機能が用意されているが*2、COBRApyでは見当たらない。すでにMMINTEというライブラリが同様の機能をPythonで実装しているようなので、それを使ってみる。
MMINTEのgithubページにおいてあったスクリプトを参考に以下のコードを走らせるとエラーが出た。
import mminte species = ["Bacteroides_caccae_ATCC_43185", "Lactobacillus_rhamnosus_GG_ATCC_53103"] model_filenames = ["AGORA/CurrentVersion/AGORA_1_03/AGORA_1_03_sbml/{s}.xml".format(s=s) for s in species] pairs = mminte.get_all_pairs(model_filenames) pair_model_filenames = mminte.create_interaction_models(pairs, ".") pair_model_filenames
RuntimeError: The heuristic for discovering an external compartment relies on names and boundary reactions. Yet, there are neither compartments with recognized names nor boundary reactions in the model.
エラーメッセージを見ても原因がよくわからなかったので、MMINTEのコードをチェックしてみた。上記の関数mminte.create_interaction_modelsの中では関数create_community_modelが主要な働きをしている様子。
import six community_compartment = "u" # lumen # MetaboliteのSuffixをチェックする # e.g. ****[e] はExtracellular metaboliteを表すbigg形式のSuffix modelseed_suffix_re = re.compile(r'_([ce])$') bigg_suffix_re = re.compile(r'\[([cepu])\]$') def _id_type(object_id): if re.search(bigg_suffix_re, object_id) is not None: return 'bigg' if re.search(modelseed_suffix_re, object_id) is not None: return 'modelseed' return 'unknown' # Metaboliteの場所 (Compartment) を変更する関数 def _change_compartment(id_str, compartment, id_type): if id_type == 'bigg': return re.sub(bigg_suffix_re, '[{0}]'.format(compartment), id_str) elif id_type == 'modelseed': return re.sub(modelseed_suffix_re, '_{0}'.format(compartment), id_str) return '{0}_{1}'.format(id_str, compartment) def create_community_model(model_files, model_names): # 空のGEMsを作成する community = cobra.Model('Community') # Exchange reaction IDのリスト (後に重複しないようにリストを作成する) exchange_reaction_ids = set() # 各生物のGEMsを追加する for model_file, model_name in zip(model_files, model_names): # GEMの読み込み model = cobra.io.read_sbml_model(model_file) # linear_coefficients = cobra.util.solver.linear_reaction_coefficients(model) objective_id = '{0}_{1}'.format(model_name, six.next(six.iterkeys(linear_coefficients)).id) objective_value = six.next(six.itervalues(linear_coefficients)) # 各MetaboliteのCompartmentの情報が存在するかMetabolite IDのSuffixをもとに確認 for metabolite in model.metabolites: metabolite.notes['type'] = _id_type(metabolite.id) unknown = model.metabolites.query(lambda x: 'unknown' in x['type'], 'notes') if len(unknown) > 0: raise Exception('Unknown compartment suffixes found in metabolites for {0}'.format(model_file)) # 各細菌のGEMsのExchange reactionのリストを取得する exchange_reactions = model.exchanges # 細菌間で共通の区画`ENV`を設定する for reaction in exchange_reactions: if reaction.id not in exchange_reaction_ids: exchange_reaction_ids.add(reaction.id) metabolite = six.next(six.iterkeys(reaction.metabolites)).copy() # MetaboliteのCompartmentを"e"からlumen ("u")へ変更 metabolite.compartment = community_compartment metabolite.id = _change_compartment(metabolite.id, community_compartment, metabolite.notes['type']) community.add_metabolites(metabolite) rxn = community.add_boundary(metabolite) # This is slow on cobra06 rxn.id = reaction.id # Keep same ID as species model for reaction in model.reactions: if reaction in exchange_reactions: # 区画`ENV`と`EX`の間の物質交換を定義する e_metabolite = six.next(six.iterkeys(reaction.metabolites)) u_metabolite = community.metabolites.get_by_id( _change_compartment(e_metabolite.id, community_compartment, e_metabolite.notes['type'])) reaction.add_metabolites({u_metabolite: 1.}) # Reaction ID の名前を変える (各細菌の反応を区別するため) reaction.id = '{0}_TR_{1}'.format(model_name, reaction.id[3:]) # Strip "EX_" prefix reaction.bounds = (-1000., 1000.) else: # Reaction ID の名前を変える (各細菌の反応を区別するため) reaction.id = '{0}_{1}'.format(model_name, reaction.id) # Metabolite ID の名前を変える (各細菌の物質を区別するため) for metabolite in model.metabolites: if metabolite.compartment != community_compartment: metabolite.id = '{0}_{1}'.format(model_name, metabolite.id) # あたらしいGEMに追加する community += model # 各細菌のBiomass functionを目的関数に追加 objective_reaction = community.reactions.get_by_id(objective_id) cobra.util.solver.set_objective(community, {objective_reaction: objective_value}, additive=True) # GEMに名前をつける community.id = 'x'.join(model_names) return community
cobra.Model.add_metabolitesで上記のエラーメッセージが出ていた。調べると以下のメソッドが関わっていそうなことがわかった。cobra.Model.compartments == "u"が細胞外を表す区画だと認識されていないようだ。おそらくMMINTEが開発された頃からCOBRAの仕様が変わったものと思われる。
import logging LOGGER = logging.getLogger(__name__) # cobra.Model.compartmentsから細胞外のCompartmentを見つける def find_external_compartment(model): if model.boundary: counts = pd.Series(tuple(r.compartments)[0] for r in model.boundary) most = counts.value_counts() most = most.index[most == most.max()].to_series() else: most = None like_external = cobra.medium.annotations.compartment_shortlist["e"] + ["e"] matches = pd.Series( [co in like_external for co in model.compartments], index=model.compartments ) if matches.sum() == 1: compartment = matches.index[matches][0] LOGGER.info( "Compartment `%s` sounds like an external compartment. " "Using this one without counting boundary reactions" % compartment ) return compartment elif most is not None and matches.sum() > 1 and matches[most].sum() == 1: compartment = most[matches[most]][0] LOGGER.warning( "There are several compartments that look like an " "external compartment but `%s` has the most boundary " "reactions, so using that as the external " "compartment." % compartment ) return compartment elif matches.sum() > 1: raise RuntimeError( "There are several compartments (%s) that look " "like external compartments but we can't tell " "which one to use. Consider renaming your " "compartments please." ) if most is not None: return most[0] LOGGER.warning( "Could not identify an external compartment by name and" " choosing one with the most boundary reactions. That " "might be complete nonsense or change suddenly. " "Consider renaming your compartments using " "`Model.compartments` to fix this." ) # No info in the model, so give up raise RuntimeError( "The heuristic for discovering an external compartment " "relies on names and boundary reactions. Yet, there " "are neither compartments with recognized names nor " "boundary reactions in the model." )
cobra.Model.copartmentsが以下のどれかだと、外部環境を表す区画として認識してくれる様子。
cobra.medium.annotations.compartment_shortlist["e"] + ["e"]
['extracellular', 'extraorganism', 'out', 'extracellular space', 'extra organism', 'extra cellular', 'extra-organism', 'external', 'external medium', 'e']
community_compartment = "out"としたら問題なく走るようになった。
結果
species = ["Bacteroides_caccae_ATCC_43185", "Lactobacillus_rhamnosus_GG_ATCC_53103"] model_filenames = ["AGORA/CurrentVersion/AGORA_1_03/AGORA_1_03_sbml/{s}.xml".format(s=s) for s in species] # 共培養GEMをつくる community = create_community_model(model_filenames, species) # Mediumの設定：Arabinogalactoseを除外した場合 medium = community.medium for k in medium.keys(): medium[k] = 0 for ex,val in zip(dmem["Exchange reaction ID"], dmem["Uptake rate"]): if community.reactions.has_id(ex): medium[ex] = val medium["EX_ala_L(e)"] = 0 medium["EX_arabinogal(e)"] = 0 community.medium = medium op = community.optimize() print("DMEM with Vitamin K, Hemin: {g}".format(g=op.objective_value)) print(" {s}: {g}".format(s=community.notes["species"][0]["name"], g=op.fluxes[community.notes["species"][0]["objective"]])) print(" {s}: {g}".format(s=community.notes["species"][1]["name"], g=op.fluxes[community.notes["species"][1]["objective"]])) # Mediumの設定：Arabinogalactoseありの場合 medium = community.medium for k in medium.keys(): medium[k] = 0 dmem = pd.read_table("Magnusdottir-2017-Nat_Biotechnol/Medium-DMEM_6429.tsv") for ex,val in zip(dmem["Exchange reaction ID"], dmem["Uptake rate"]): if community.reactions.has_id(ex): medium[ex] = val medium["EX_ala_L(e)"] = 0 community.medium = medium op = community.optimize() print("DMEM with Vitamin K, Hemin, Arabinogalactan: {g}".format(g=op.objective_value)) print(" {s}: {g}".format(s=community.notes["species"][0]["name"], g=op.fluxes[community.notes["species"][0]["objective"]])) print(" {s}: {g}".format(s=community.notes["species"][1]["name"], g=op.fluxes[community.notes["species"][1]["objective"]]))
各細菌の生育率は以下のように推定された。Lactobacillus rhamnosus GG単独のGEMにFBAを適用すると、Arabinogalactanが無いとDMEM培地では生育しないと推定されたが*3、共培養するとArabinogalactanなしでも生育すると推定された。
- Arabinogalactose あり
- Bacteroides_caccae_ATCC_43185: 0.4386544423746743
- Lactobacillus_rhamnosus_GG_ATCC_53103: 0.22080916149500934
- Arabinogalactose なし
- Bacteroides_caccae_ATCC_43185: 0.24543405369998675
- Lactobacillus_rhamnosus_GG_ATCC_53103: 0.11140717070316093
*1:論文ではlumenと呼んでいる
*2:https://opencobra.github.io/cobratoolbox/latest/tutorials/tutorialMicrobeMicrobeInteractions.html
*3:論文と同じ＆前回記事参照
