8  通量取样

8.1 基本用法

开始使用通量采样的最简单的方法是在flux_analysis 子模块的sample 功能。sample至少接受俩个参数：cobra模型以及要生成的样本数量。

from cobra.io import load_model
from cobra.sampling import sample

model = load_model("textbook")
s = sample(model, 100)
s.head()

	ACALD	ACALDt	ACKr	ACONTa	ACONTb	ACt2r	ADK1	AKGDH	AKGt2r	ALCD2x	...	RPI	SUCCt2_2	SUCCt3	SUCDi	SUCOAS	TALA	THD2	TKT1	TKT2	TPI
0	-3.799475	-2.498751	-1.421441	7.540352	7.540352	-1.421441	16.929342	6.483571	-0.338296	-1.300724	...	-1.499339	14.876146	15.591381	566.653991	-6.483571	1.464935	1.833969	1.464935	1.451100	8.457706
1	-1.664563	-1.159660	-0.104808	10.096309	10.096309	-0.104808	21.034965	9.031672	-0.102474	-0.504903	...	-1.458747	2.293237	3.284039	484.075418	-9.031672	1.434034	48.090638	1.434034	1.424095	7.472630
2	-1.803558	-1.667910	-2.144206	8.824075	8.824075	-2.144206	11.796768	2.727459	-0.121884	-0.135648	...	-1.705485	1.930752	2.053874	721.454310	-2.727459	1.681770	75.179992	1.681770	1.672234	7.402820
3	-1.448489	-0.122990	-1.631803	9.828801	9.828801	-1.631803	22.691294	5.339959	-0.201254	-1.325499	...	-1.986551	10.169471	11.442413	684.373359	-5.339959	1.892144	36.876984	1.892144	1.854180	7.635212
4	-0.751012	-0.234146	-1.534250	10.996159	10.996159	-1.534250	2.346065	5.583547	-0.695630	-0.516866	...	-3.448945	1.012854	2.196951	624.064170	-5.583547	3.400375	19.158341	3.400375	3.380843	6.447869

5 rows × 95 columns

默认情况下，sample 使用基于此处 method presented here介绍的方法的optgp方法，因为它适用于较大的模型，并且可以并行运行。默认情况下，采样器使用单个进程。这可以通过使用 processes参数进行更改

print("One process:")
%time s = sample(model, 1000)
print("Two processes:")
%time s = sample(model, 1000, processes=2)

One process:
CPU times: total: 20.8 s
Wall time: 15.8 s
Two processes:
CPU times: total: 1.92 s
Wall time: 34 s

或者，您也可以通过将方法设置为“achr”来使用 Artificial Centering Hit-and-Run 进行采样。 ‘achr’ 不支持并行执行，但收敛性好，几乎是马尔可夫式的。

s = sample(model, 100, method="achr")

一般来说，设置采样器的成本很高，因为初始搜索方向是通过解决许多线性规划问题生成的。因此，我们建议一次性生成尽可能多的样品。但是，这可能需要对采样过程进行更精细的控制，如下一节所述。

8.2 高级使用方法

8.2.1 采样目标

通过直接使用采样器类，可以在较低级别上控制采样过程。

from cobra.sampling import OptGPSampler, ACHRSampler

这两个采样器类都具有标准化的接口，并采用一些额外的参数。例如，“变薄”因素。“减薄”是指每 n 次迭代仅记录样本。较高的稀化系数意味着较少的相关样本，但计算时间也较长。默认情况下，采样器使用 100 的稀疏因子，这会创建大致不相关的样本。如果您想要更少的样品但更好的混合，请随时增加此参数。如果要研究自己的模型的收敛性，则可能需要将其设置为 1 以获取所有迭代。

achr = ACHRSampler(model, thinning=10)

OptGPSampler有一个额外的 processes参数，用于指定用于创建并行采样链的进程数。这应该按照 CPU 内核的顺序排列，以实现最大效率。如前所述，由于生成初始搜索方向，类初始化可能需要几分钟时间。另一方面，采样速度很快。

optgp = OptGPSampler(model, processes=4)

8.2.2 采样和验证

两个采样器都有一个示例函数，该函数从初始化的对象生成样本，其作用类似于上面描述的 sample函数，只是这次它只接受一个参数，即样本数。对于OptGPSampler，样本数应为进程数的倍数，否则将自动增加到最接近的倍数。

s1 = achr.sample(100)

s2 = optgp.sample(100)

You cansampleamplerepeatedly and both samplers are optimized to generate large amount of samples without falling into "numerical traps". All sampler objects have avalidate` function in order to check if a set of points are feasible and give detailed information about feasibility violations in a form of a short code denoting feasibility. Here the short code is a combination of any of the following letters:

• “v” - valid point
• “l” - lower bound violation
• “u” - upper bound violation
• “e” - equality violation (meaning the point is not a steady state)

For instance for a random flux distribution (should not be feasible):

import numpy as np

bad = np.random.uniform(-1000, 1000, size=len(model.reactions))
achr.validate(np.atleast_2d(bad))

array(['le'], dtype='<U3')

对于我们生成的样本：

achr.validate(s1)

array(['v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',
       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v'], dtype='<U3')

尽管大多数模型在数值上足够稳定，采样器应该只生成有效样本，但我们仍然敦促检查这一点。validate 非常快，即使对于大型模型和许多样本也能快速工作。如果发现无效样本，则不必重新运行整个样本，但可以从示例 DataFrame 中排除它们。

s1_valid = s1[achr.validate(s1) == "v"]
len(s1_valid)

100

8.2.3 批量抽样

采样器对象用于生成数十亿个样本，但是在使用基因组规模模型时，使用“样本”功能可能会很快填满您的 RAM。在这里，采样器对象的sample方法可能会派上用场。 batch 有两个参数，即每批中的样本数和批数。举个小例子，这是有道理的。

让我们假设我们想要量化样本中将增长的比例。为此，我们可能需要生成 10 批，每批 50 个样本，并测量单个 100 个样本中显示大于 0.1 的百分比。最后，我们要计算这些单个百分比的平均值和标准差。

counts = [np.mean(s.Biomass_Ecoli_core > 0.1) for s in optgp.batch(100, 10)]
print("Usually {:.2f}% +- {:.2f}% grow...".format(
    np.mean(counts) * 100.0, np.std(counts) * 100.0))

Usually 3.00% +- 0.77% grow...

8.3 添加约束

通量采样将遵循模型中定义的其他约束条件。例如，我们可以添加一个约束，以与上一节类似的方式强制增长。

co = model.problem.Constraint(model.reactions.Biomass_Ecoli_core.flux_expression, lb=0.1)
model.add_cons_vars([co])

请注意，这仅用于演示目的。通常，您可以直接设置反应的下限，而不是创建新约束。

s = sample(model, 10)
print(s.Biomass_Ecoli_core)

0    0.105278
1    0.126172
2    0.108318
3    0.105436
4    0.100796
5    0.154830
6    0.155173
7    0.177615
8    0.193415
9    0.148553
Name: Biomass_Ecoli_core, dtype: float64

正如我们所看到的，我们的新约束被接受了。