はじめに
COBRA とは
COBRA（Constraint Based Reconstruction and Analysis）とは、ある制約条件をもたせた上でのネットワークの再構築及び解析を行うツールであり、生物学の分野では専ら、微生物や哺乳細胞、微細藻類などの宿主の代謝を数理的にモデリングするツールとして使われています。
代謝ネットワークは、代謝物をノード、酵素（反応）をエッジとしたグラフとして描くことができ、反応の方向性（可逆/不可逆）もあって有向グラフとなります。また、宿主は無尽蔵に栄養を取り込めるわけではなく、ある程度決まった速度で酸素やグルコースを消費しており、実際の計算では酸素消費速度（OUR）やグルコース消費速度（GUR）を設定します。このあたりが「制約条件」と呼ばれる所以ですね。
ライブラリ
openCOBRAというプロジェクトで、各種プログラミング言語に対応したオープンソースライブラリの構築が現在進行形で進められています。具体的には、MATLAB, Python, Julia でオープンソースライブラリが準備されています。
問題意識
オープンソースであるがゆえ（？）、関数として完備されていなかったり、付属ドキュメントの情報量が不足していたりすることも多く、初学者にとってとっつきにくいパッケージです。特に、生物工学の文脈で用いるパッケージであるため、生粋の情報系の方よりは分子生物学・合成生物学の方の利用が期待されます。ちょっとしたプログラミングの難しさでライブラリの使用を断念してしまうのは非常にもったいないです。
本記事で目指すこと
openCOBRA のうち、著者が用いている The COBRA Toolbox (MATLAB) と COBRApy (Python) のパッケージで、ドキュメント説明やチュートリアルが不十分な点を詳説していきたいと思います。本記事は必要に応じて随時更新していく予定です。
The COBRA Toolbox (MATLAB)
findRxnsFromMets
特定の代謝物を含む反応を一覧化する関数です。例えば、大腸菌の代謝モデルにおいて、代謝物としてフマル酸を選択すると、フマル酸が反応物もしくは生成物に含まれる反応を全てリストアップしてくれます（コハク酸デヒドロゲナーゼ：コハク酸→フマル酸、フマラーゼ：フマル酸→L-リンゴ酸、など）。
ただし、この関数は代謝物を1種類しか指定しないと、間違った答えを返すようです（バグ報告例あり）。明示的なエラーを返してくれないのが厄介ですが、明らかに間違った答えを返すので、プログラム上のミスが有ると想定されます。なお、代謝物を2種類以上指定すると、指定した代謝物のうち1種類でも含む反応を全て返してくれますので、正しく動作するようです。
ここで、大元の findRxnsFromMets.m のソースコードの119行目を眺めてみましょう。
totals = sum(model.S(index,:) ~= 0,1); %RF: This seems wrong as it omits production of metabolites.
よく見るとコメントが書かれていますね。"This seems wrong as it omits production of metabolites."（代謝物の生成を無視しているため、このコードは間違っているように見えます。）ということで、この部分の問題が解決されていないようです。
解決策1： 訂正した findRxnsFromMets.m ファイルを作成する
これが真っ当にできればよかったのですが、どのように訂正すべきか分かっておらず、まだ手つかずです。分かり次第、追記します。
解決策2： 無理やり代謝物を2種類以上指定して出力させる
先程も記載したように、代謝物を2種類以上指定すれば、正しく返してくれますので、関与する反応を知りたい本命の代謝物と、関与する反応が既知の代謝物を指定してやることで、一応間違いを回避できそうです。もちろん、出力された反応一覧から、既知の代謝物が関与する反応を後で差し引かなければならないですが……。
COBRApy (Python)
cobra.flux_analysis.variability.flux_variability_analysis
Flux Variability Analysis (FVA；フラックス変動解析）は、ある条件を満たす各反応フラックス解の範囲を提示してくれます。この関数の引数として指定されている fraction_of_optimum ですが、ドキュメントには
fraction_of_optimum (float, optional) – The fraction of optimum which must be maintained. The original objective reaction is constrained to be greater than maximal value times the fraction_of_optimum (default 1.0).
と記載されています。これは、比増殖速度 $μ$ が、fraction_of_optimum で設定した値 $a$（∈ [0, 1]）に最大比増殖速度 $μ_{max}$ を乗じた値以上になる（$μ ≧ aμ_{max}$）フラックス解を求めることになります。極端な例で言えば、fraction_of_optimum = 1 のとき、$μ ≧ μ_{max}$ を満たすフラックス解を求めることになるため、それは $μ = μ_{max}$ を満たすフラックス解を求めることと同義になります。一方で、fraction_of_optimum = 0 のとき、$μ ≧ 0$ を満たすフラックス解を求めることになるため、あらゆるフラックスの組み合わせが解となります。