Biopython 中的 Genepop 支持

The Genepop 模块允许使用 Python 接口访问 Genepop 功能。这意味着 Genepop 的绝大多数方法（Hardy-Weinberg 平衡的精确检验、群体分化、基因型不平衡、F 统计量、零等位基因频率、微卫星的基于等位基因大小的统计量等等）现在都可以从 Python 中访问。由于代码只是一个包装器，因此需要安装 Genepop。

Genepop 文件的解析器也可用（并在教程中进行了说明）。Genepop 格式的文件可以在没有 Genepop 的情况下进行处理。

提供了两个接口：一个通用、更复杂、更高效的接口 (GenePopController) 和一个简化、更易于使用、不完整、效率不高的版本 (EasyController)。EasyController 可能无法处理非常大的文件，因为它有接口。另一方面，它提供了计算一些非常简单的统计信息（如等位基因计数）的实用函数，这些函数在通用接口中不可直接使用。

更复杂的接口假设更熟练的 Python 开发人员（例如，通过使用迭代器），目前还没有文档。但即使对于经验丰富的 Python 开发人员，只要在 EasyController 中公开所需的功能，并且其性能被认为是可以接受的，EasyController 也可以很方便。

有关用于计算的方法的详细信息，请查看 Genepop 文档，该文档提供了指向所有派生计算论文的指针。

EasyController 教程

安装

为了使用控制器，必须在系统中安装 Genepop，可以从 这里 下载。

在我们开始之前，让我们测试一下安装（为此，您需要一个 Genepop 格式的文件）

from Bio.PopGen.GenePop.EasyController import EasyController

ctrl = EasyController(your_file_here)
print(ctrl.get_basic_info())

用您的文件名和路径替换 your_file_here。如果您收到 IOError: Genepop not found，则 Biopython 找不到您的 Genepop 可执行文件。如果 Genepop 不在 PATH 上，您可以将其添加到构造函数行中，例如

ctrl = EasyController(your_file_here, path_to_genepop_here)

如果一切正常，现在我们可以继续使用 Genepop。对于以下示例，我们将使用 Genepop 文件 big.gen，该文件与单元测试一起提供。我们还将假设存在一个用选定的相关文件初始化的 ctrl 对象。

我们从获取一些基本信息开始

pop_names, loci_names = ctrl.get_basic_info()

返回文件中可用的群体名称和基因座名称列表。

警告：大多数现有的 Genepop 文件提供了关于群体名称的错误数据。在许多情况下，可能不信任这些信息。评估群体信息，大多数情况下是通过文件中群体的相对位置来完成的，而不是名称。因此，文件中的第一个群体是索引 0，第二个是索引 1，依此类推…

统计数据

杂合性

让我们获取特定群体和特定等位基因的杂合性信息

(exp_homo, obs_homo, exp_hetero, obs_hetero) = ctrl.get_heterozygosity_info(0, "Locus2")

将获得群体 0 和基因座 2 的预期和观察到的纯合性和杂合性（对于文件 big.gen，如果您使用的是其他文件，请相应地调整群体位置和基因座名称）。

现有等位基因

可以获得特定群体中特定基因座的所有等位基因列表

allele_list = ctrl.get_alleles(0, "Locus2")

allele_list 将为 [3, 20]（即，等位基因 3 和 20 在该群体上）。

等位基因的数量只需使用 len(allele_list) 获取即可。

还可以获取所有群体中特定基因座的所有等位基因列表

all_allele_list = ctrl.get_alleles_all_pops("Locus2")

all_allele_list 将为 [3, 20]。

等位基因和基因型频率

可以获取特定群体中等位基因的频率

allele_data = ctrl.get_allele_frequency(0, "Locus2")

allele_data 将为 (62, {3: 0.88700000000000001, 20: 0.113})。也就是说，有 62 个基因。88.7% 为 3，11.3% 为 20。

我们可以获取基因型（二倍体数据）的类似信息。还将报告预期频率

genotype_list = ctrl.get_genotype_frequency(0, "Locus2")

genotype_list 将为：[(3, 3, 24, 24.3443), (20, 3, 7, 6.3114999999999997), (20, 20, 0, 0.34429999999999999)]

让我们解释第一个元素：有 24 个个体具有 (3, 3) 的基因型，而具有该基因型的个体预期数量为 24.2443。

F 统计量

让我们从一般的多基因座 F 统计量开始

Fis, Fst, Fit = ctrl.get_multilocus_f_stats()

这将获取多基因座 Fis、Fst 和 Fit。

让我们获取这些数据（以及更多）每个基因座

Fis, Fst, Fit, Qintra, Qinter = ctrl.get_f_stats("Locus2")

这将获取单基因座 Fis、Fst 和 Fit、Qintra 和 Qinter。

在 Fis 部分下面有专门的 Fst 和 Fis 部分（其中引入了成对和群体特定变体）。在 Fis 部分中，解释了 Qintra 和 Qinter。

Fst

让我们获取特定基因座的成对 Fst

pair_fst = ctrl.get_avg_fst_pair_locus("Locus4")

将返回一个映射，其中键是 population1、population2（群体 ID）组成的对。population2 始终小于 population1。例如：群体 0 和群体 3 之间基因座 4 的成对 Fst 由 pair_fst[(3,0)] 给出。

您还可以获取多基因座成对 Fst 估计

multilocus_fst = ctrl.get_avg_fst_pair()

这将返回与上面相同的结构，但使用多基因座成对 Fst。

Fis

现在我们将获取特定基因座/群体的 Fis 以及其他一些统计数据

allele_dict, summary_fis = ctrl.get_fis(0, "Locus2")

让我们详细了解 get_fis 的输出

summary_fis = (62, -0.1111, -0.11269999999999999)

allele_dict = {3: (55, 0.8871, -0.1111), 20: (7, 0.1129, -0.1111)}

summary_fis 持有一个三元组，其中包含：等位基因总数、Cockerham 和 Weir Fis、Robertson 和 Hill Fis。

allele_dict 为每个等位基因（每个等位基因都是键）保存：等位基因的重复次数、频率以及 Cockerham 和 Weir Fis。

因此，从上面的结果中可以读出以下内容：有 62 个基因，有 2 个不同的等位基因（55 个为 3 型，7 个为 20 型）。3 型的频率为 0.89，20 型的频率为 0.11。所有 CW Fis 都是 -0.111，RH Fis 是 -0.112。

现在让我们获取多基因座 Fis

pop_list = ctrl.get_avg_fis()

pop_list 将返回每个群体的一个元素。每个元素都是一个四元组，包含

1. 群体名称（同样，群体名称不可信）
2. 1 - QIntra：个体之间的基因多样性
3. 1 - QInter：群体内部个体之间的基因多样性
4. Fis

迁移

我们可以获取移民人数的估计值

samp_size, priv_allele_freq, mig10, mig25, mig50, migcorr = ctrl.estimate_nm()

samp_size 是平均样本量，priv_allele_freq 是私有等位基因的平均频率，mig10 是 Ne=10 时移民人数，mig25 是 Ne=25 时移民人数，mig50 是 Ne=50 时移民人数，migcorr 是在校正预期规模后的移民人数。

测试

测试通常计算量很大，因为它们通常基于马尔可夫链算法。在某些情况下，可以使用完全枚举方法，但这些方法只能应用于等位基因数量非常少的基因座。这意味着大多数测试需要相当长的时间才能完成。

有关下面马尔可夫链参数（去记忆、分批和迭代）的更多详细信息，请参阅 Genepop 手册。还要咨询手册以了解完全枚举何时适用。

Hardy-Weinberg 平衡

让我们从测试每个群体中每个基因座的 Hardy-Weinberg 平衡开始

loci_map = ctrl.test_hw_pop(1, "excess")

第二个参数可以是 probability、excess 或 deficiency。probability 是标准的 Haldane HW 检验。当您对杂合子缺乏感兴趣时，使用 deficiency；当您对过量感兴趣时，使用 excess。

输出是一个映射，其中键是基因座名称。内容是一个元组，包含 P 值、标准误差、Fis（Weir 和 Cockerham）、Fis（Robertson 和 Hill）以及步骤。

pop_test, loc_test, all_test = ctrl.test_hw_global("deficiency")

当您对杂合子缺乏感兴趣时，使用 deficiency；当您对过量感兴趣时，使用 excess。probability 在这里不适用，就像在 test_hw_pop 中一样。

输出是一个三元组

1. pop_test 是一个列表，每个群体有一个元素，包括 P 值、标准误差和开关。
2. loc_test 是同一个列表，但每个基因座有一个元素，包括基因座名称、P 值、标准误差和开关。
3. all_test 是总体结果，包含一个 P 值、标准误差和开关三元组。

连锁不平衡

我们可以使用对数似然比统计量（G 检验）测试两个基因座是否处于连锁不平衡状态。

chi2, df, pval = ctrl.test_ld_all_pair(
    "Locus1", "Locus2", dememorization=1000, batches=10, iterations=100
)

返回 G 统计量的卡方值、自由度和 P 值。

距离隔离 (IBD)

距离隔离 (IBD) 分析需要一种特殊的 Genepop 文件形式

1. 每个群体只有一个个体
2. 个体的名称必须是其坐标

示例

...
Pop
0 15, 0201  0303 0102 0302 1011
Pop
0 30, 0202  0301 0102 0303 1111
Pop
0 45, 0102  0401 0202 0102 1010
Pop
0 60, 0103  0202 0101 0202 1011
Pop
0 75, 0203  0204 0101 0102 1010
POP
15 15, 0102 0202 0201 0405 0807
...

请注意，我们正在使用的示例文件不能用于这种情况。

IBD 分析只有一个调用

estimate, distance, (a, b), (bb, bblow, bbhigh) = ctrl.calc_ibd(
    self, is_diplo=True, stat="a", scale="Log", min_dist=0.00001
)

• is_diplo 指定数据是二倍体 (True) 还是单倍体 (False)。
• stat 既可以是 'a' 也可以是 'e'（有关详细信息，请参阅 Genepop 手册）。
• scale 既可以是 'Log' 也可以是 'Linear'。'Log' 用于二维坐标，'Linear' 用于一维坐标。

仅使用大于 min_dist 的成对比较来计算回归系数。

该方法返回

• estimate，一个三角矩阵，包含根据所选统计量计算的样本之间的遗传距离。
• distance，一个三角矩阵，包含样本之间的距离（对数或线性）。
• a 和 b *是回归的拟合参数。*bblow 和 bbhigh 是 *b* 参数的引导置信区间（*bb* 应该非常接近 *b*）。

对三角矩阵的解释应该这样进行：在 Python 中，矩阵是用数字的列表列表实现的，就像这样

[[0.1], [0.2, 0.3], [0.4, 0.5, 0.6]]

上面的数据结构对应于以下三角矩阵

      1    2    3
2   0.1
3   0.2  0.3
4   0.4  0.5  0.6