Python-生物信息学秘籍第三版-二-

布客飞龙
122 阅读51分钟

Python 生物信息学秘籍第三版（二）

原文：annas-archive.org/md5/9694cf42f7d741c69225ff1cf52b0efe

译者：飞龙

协议：CC BY-NC-SA 4.0

第五章：高级 NGS 数据处理

如果您使用下一代测序（NGS）数据，您会知道质量分析和处理是获取结果中的两个主要时间消耗。在本章的第一部分，我们将通过使用包含亲属信息的数据集来深入探讨 NGS 分析 - 在我们的情况下，是一个母亲、一个父亲和大约 20 个后代的数据集。这是进行质量分析的常见技术，因为家系信息将允许我们推断出我们的过滤规则可能产生的错误数量。我们还将利用同一数据集来查找基于现有注释的基因组特征。

本章的最后一个配方将深入探讨使用 NGS 数据的另一个高级主题：宏基因组学。我们将使用 QIIME2，一个用于宏基因组学的 Python 包，来分析数据。

如果您使用 Docker，请使用 tiagoantao/bioinformatics_base 镜像。有关 QIIME2 内容的特殊设置过程将在相关配方中讨论。

本章包括以下配方：

  • 准备用于分析的数据集

  • 使用门利因错误信息进行质量控制

  • 使用标准统计方法探索数据

  • 从测序注释中找到基因组特征

  • 使用 QIIME2 进行宏基因组学

准备用于分析的数据集

我们的起点将是一个 VCF 文件（或等效文件），其中包含由基因分析工具（在我们的情况下是基因组分析工具包（GATK））进行的调用，包括注释信息。因为我们将过滤 NGS 数据，我们需要可靠的决策标准来调用一个位点。那么，我们如何获取这些信息？一般来说，我们做不到，但如果我们确实需要这样做，有三种基本方法：

  • 使用更强大的测序技术进行比较 - 例如，使用 Sanger 测序验证 NGS 数据集。这种方法成本高昂，只能用于少数基因位点。

  • 对于测序密切相关的个体，例如两个父母及其后代。在这种情况下，我们使用门利因遗传规则来决定某个调用是否可接受。这是人类基因组计划和安非蚊 1000 基因组计划均采用的策略。

  • 最后，我们可以使用模拟。这种设置不仅相当复杂，而且可靠性存疑。这更多是一个理论选项。

在本章中，我们将使用第二个选项，基于安非蚊 1000 基因组计划。该项目提供了基于蚊子杂交的信息。一个杂交会包括父母（母亲和父亲）以及最多 20 个后代。

在这个配方中，我们将准备我们的数据，以便在后续配方中使用。

准备就绪

我们将以 HDF5 格式下载文件以加快处理速度。请注意，这些文件相当大；您需要良好的网络连接和充足的磁盘空间：

wget -c ftp://ngs.sanger.ac.uk/production/ag1000g/phase1/AR3/variation/main/hdf5/ag1000g.phase1.ar3.pass.3L.h5
wget -c ftp://ngs.sanger.ac.uk/production/ag1000g/phase1/AR3/variation/main/hdf5/ag1000g.phase1.ar3.pass.2L.h5

这些文件有四个交叉，每个交叉大约有 20 个后代。我们将使用染色体臂 3L 和 2L。在这一阶段，我们也计算孟德尔误差（这是下一个食谱的主题，因此我们将在那时详细讨论）。

相关的笔记本文件是Chapter04/Preparation.py。目录中还有一个名为samples.tsv的本地样本元数据文件。

如何操作……

下载数据后，按照以下步骤操作：

  1. 首先，从一些导入开始：

    import pickle
import gzip
import random
import numpy as np
import h5py
import pandas as pd

  2. 让我们获取样本元数据：

    samples = pd.read_csv('samples.tsv', sep='\t')
print(len(samples))
print(samples['cross'].unique())
print(samples[samples['cross'] == 'cross-29-2'][['id', 'function']])
print(len(samples[samples['cross'] == 'cross-29-2']))
print(samples[samples['function'] == 'parent'])

我们还打印一些关于我们将要使用的交叉和所有父母的基本信息。

  1. 我们准备根据其 HDF5 文件处理染色体臂 3L：

    h5_3L = h5py.File('ag1000g.crosses.phase1.ar3sites.3L.h5', 'r')
samples_hdf5 = list(map(lambda sample: sample.decode('utf-8'), h5_3L['/3L/samples']))
calldata_genotype = h5_3L['/3L/calldata/genotype']
MQ0 = h5_3L['/3L/variants/MQ0']
MQ = h5_3L['/3L/variants/MQ']
QD = h5_3L['/3L/variants/QD']
Coverage = h5_3L['/3L/variants/Coverage']
CoverageMQ0 = h5_3L['/3L/variants/CoverageMQ0']
HaplotypeScore = h5_3L['/3L/variants/HaplotypeScore']
QUAL = h5_3L['/3L/variants/QUAL']
FS = h5_3L['/3L/variants/FS']
DP = h5_3L['/3L/variants/DP']
HRun = h5_3L['/3L/variants/HRun']
ReadPosRankSum = h5_3L['/3L/variants/ReadPosRankSum']
my_features = {
    'MQ': MQ,
    'QD': QD,
    'Coverage': Coverage,
    'HaplotypeScore': HaplotypeScore,
    'QUAL': QUAL,
    'FS': FS,
    'DP': DP,
    'HRun': HRun,
    'ReadPosRankSum': ReadPosRankSum
}
num_features = len(my_features)
num_alleles = h5_3L['/3L/variants/num_alleles']
is_snp = h5_3L['/3L/variants/is_snp']
POS = h5_3L['/3L/variants/POS']

  2. 计算孟德尔误差的代码如下：

    #compute mendelian errors (biallelic)
def compute_mendelian_errors(mother, father, offspring):
    num_errors = 0
    num_ofs_problems = 0
    if len(mother.union(father)) == 1:
        # Mother and father are homogenous and the            same for ofs in offspring:
            if len(ofs) == 2:
                # Offspring is het
                num_errors += 1
                num_ofs_problems += 1
            elif len(ofs.intersection(mother)) == 0:
                # Offspring is homo, but opposite from parents
                num_errors += 2
                num_ofs_problems += 1
    elif len(mother) == 1 and len(father) == 1:
        # Mother and father are homo and different
        for ofs in offspring:
            if len(ofs) == 1:
                # Homo, should be het
                num_errors += 1
                num_ofs_problems += 1
    elif len(mother) == 2 and len(father) == 2:
        # Both are het, individual offspring can be anything
        pass
    else:
        # One is het, the other is homo
        homo = mother if len(mother) == 1 else father
        for ofs in offspring:
            if len(ofs) == 1 and ofs.intersection(homo):
                # homo, but not including the allele from parent that is homo
                num_errors += 1
                num_ofs_problems += 1
    return num_errors, num_ofs_problems

我们将在下一个食谱中讨论这个问题，使用孟德尔误差信息进行质量控制。

  1. 现在，我们定义一个支持生成器和函数，用于选择可接受的位置并累积基本数据：

    def acceptable_position_to_genotype():
    for i, genotype in enumerate(calldata_genotype):
        if is_snp[i] and num_alleles[i] == 2:
            if len(np.where(genotype == -1)[0]) > 1:
                # Missing data
                continue
            yield i
def acumulate(fun):
    acumulator = {}
    for res in fun():
        if res is not None:
            acumulator[res[0]] = res[1]
    return acumulator

  2. 现在，我们需要找到在 HDF5 文件中交叉的索引（母亲、父亲和 20 个后代）：

    def get_family_indexes(samples_hdf5, cross_pd):
    offspring = []
    for i, individual in cross_pd.T.iteritems():
        index = samples_hdf5.index(individual.id)
        if individual.function == 'parent':
            if individual.sex == 'M':
                father = index
            else:
                mother = index
        else:
            offspring.append(index)
    return {'mother': mother, 'father': father, 'offspring': offspring}
cross_pd = samples[samples['cross'] == 'cross-29-2']
family_indexes = get_family_indexes(samples_hdf5, cross_pd)

  3. 最后，我们将实际计算孟德尔误差并将其保存到磁盘：

    mother_index = family_indexes['mother']
father_index = family_indexes['father']
offspring_indexes = family_indexes['offspring']
all_errors = {}
def get_mendelian_errors():
    for i in acceptable_position_to_genotype():
        genotype = calldata_genotype[i]
        mother = set(genotype[mother_index])
        father = set(genotype[father_index])
        offspring = [set(genotype[ofs_index]) for ofs_index in offspring_indexes]
        my_mendelian_errors = compute_mendelian_errors(mother, father, offspring)
        yield POS[i], my_mendelian_errors
mendelian_errors = acumulate(get_mendelian_errors)
pickle.dump(mendelian_errors, gzip.open('mendelian_errors.pickle.gz', 'wb'))

  4. 我们现在将生成一个高效的带注释和孟德尔误差信息的 NumPy 数组：

    ordered_positions = sorted(mendelian_errors.keys())
ordered_features = sorted(my_features.keys())
num_features = len(ordered_features)
feature_fit = np.empty((len(ordered_positions), len(my_features) + 2), dtype=float)
for column, feature in enumerate(ordered_features):  # 'Strange' order
    print(feature)
    current_hdf_row = 0
    for row, genomic_position in enumerate(ordered_positions):
        while POS[current_hdf_row] < genomic_position:
            current_hdf_row +=1
        feature_fit[row, column] = my_features[feature][current_hdf_row]
for row, genomic_position in enumerate(ordered_positions):
    feature_fit[row, num_features] = genomic_position
    feature_fit[row, num_features + 1] = 1 if mendelian_errors[genomic_position][0] > 0 else 0
np.save(gzip.open('feature_fit.npy.gz', 'wb'), feature_fit, allow_pickle=False, fix_imports=False)
pickle.dump(ordered_features, open('ordered_features', 'wb'))

在这段代码中埋藏着整个章节中最重要的决定之一：我们如何权衡孟德尔误差？在我们的案例中，如果有任何误差，我们只存储 1，如果没有误差，我们存储 0。另一种选择是计数错误的数量——因为我们有最多 20 个后代，这将需要一些复杂的统计分析，而我们这里不进行这种分析。

  1. 转换思路，现在让我们从染色体臂 2L 提取一些信息：

    h5_2L = h5py.File('ag1000g.crosses.phase1.ar3sites.2L.h5', 'r')
samples_hdf5 = list(map(lambda sample: sample.decode('utf-8'), h5_2L['/2L/samples']))
calldata_DP = h5_2L['/2L/calldata/DP']
POS = h5_2L['/2L/variants/POS']

  2. 在这里，我们只关心父母：

    def get_parent_indexes(samples_hdf5, parents_pd):
    parents = []
    for i, individual in parents_pd.T.iteritems():
        index = samples_hdf5.index(individual.id)
        parents.append(index)
    return parents
parents_pd = samples[samples['function'] == 'parent']
parent_indexes = get_parent_indexes(samples_hdf5, parents_pd)

  3. 我们提取每个父母的样本 DP：

    all_dps = []
for i, pos in enumerate(POS):
    if random.random() > 0.01:
        continue
    pos_dp = calldata_DP[i]
    parent_pos_dp = [pos_dp[parent_index] for parent_index in parent_indexes]
    all_dps.append(parent_pos_dp + [pos])
all_dps = np.array(all_dps)
np.save(gzip.open('DP_2L.npy.gz', 'wb'), all_dps, allow_pickle=False, fix_imports=False)

现在，我们已经为本章的分析准备好了数据集。

使用孟德尔误差信息进行质量控制

那么，如何使用孟德尔遗传规则推断调用质量呢？让我们看看父母不同基因型配置的预期结果：

  • 对于某个潜在的双等位基因 SNP，如果母亲是 AA 且父亲也是 AA，则所有后代将是 AA。

  • 如果母亲是 AA，父亲是 TT，则所有后代必须是杂合子（AT）。他们总是从母亲那里得到一个 A，总是从父亲那里得到一个 T。

  • 如果母亲是 AA，父亲是 AT，则后代可以是 AA 或 AT。他们总是从母亲那里得到一个 A，但可以从父亲那里得到一个 A 或一个 T。

  • 如果母亲和父亲都是杂合子（AT），则后代可以是任何基因型。从理论上讲，在这种情况下我们无法做太多。

实际上，我们可以忽略突变，这是在大多数真核生物中都可以安全做到的。从我们的角度来看，突变（噪声）的数量比我们正在寻找的信号低几个数量级。

在这个例子中，我们将进行一个小的理论研究，分析分布和孟德尔错误，并进一步处理数据以便下游分析。这相关的笔记本文件是 Chapter04/Mendel.py。

如何做…

  1. 我们需要几个导入：

    import random
import matplotlib.pyplot as plt

  2. 在进行任何经验分析之前，让我们尝试理解在母亲是 AA 且父亲是 AT 的情况下我们可以提取什么信息。我们来回答这个问题，如果我们有 20 个后代，所有后代为杂合子的概率是多少？：

    num_sims = 100000
num_ofs = 20
num_hets_AA_AT = []
for sim in range(num_sims):
    sim_hets = 0
    for ofs in range(20):
        sim_hets += 1 if random.choice([0, 1]) == 1 else 0
    num_hets_AA_AT.append(sim_hets)

fig, ax = plt.subplots(1,1, figsize=(16,9))
ax.hist(num_hets_AA_AT, bins=range(20))
print(len([num_hets for num_hets in num_hets_AA_AT if num_hets==20]))

我们得到以下输出：

图 4.1 - 来自 100,000 次模拟的结果：在母亲是 AA 且父亲是杂合子的情况下，后代在某些基因座上为杂合子的数量

图 4.1 - 来自 100,000 次模拟的结果：在母亲是 AA 且父亲是杂合子的情况下，后代在某些基因座上为杂合子的数量

在这里，我们进行了 100,000 次模拟。就我而言（这是随机的，所以你的结果可能不同），我得到了零次模拟结果，其中所有的后代都是杂合子的。事实上，这些是带重复的排列，因此所有都是杂合子的概率是 或 9.5367431640625e-07——并不太可能。所以，即使对于单个后代，我们可能得到 AT 或 AA，对于 20 个后代来说，它们全部都是同一种类型的概率也非常小。这就是我们可以用来更深入解释孟德尔错误的信息。

  1. 让我们重复一下母亲和父亲都为 AT 的分析：

    num_AAs_AT_AT = []
num_hets_AT_AT = []
for sim in range(num_sims):
    sim_AAs = 0
    sim_hets = 0
    for ofs in range(20):
        derived_cnt = sum(random.choices([0, 1], k=2))
        sim_AAs += 1 if derived_cnt == 0 else 0
        sim_hets += 1 if derived_cnt == 1 else 0
    num_AAs_AT_AT.append(sim_AAs)
    num_hets_AT_AT.append(sim_hets)
fig, ax = plt.subplots(1,1, figsize=(16,9))
ax.hist([num_hets_AT_AT, num_AAs_AT_AT], histtype='step', fill=False, bins=range(20), label=['het', 'AA'])
plt.legend()

输出如下：

图 4.2 - 来自 100,000 次模拟的结果：在母亲和父亲都是杂合子的情况下，后代在某个基因座上为 AA 或杂合子的数量

图 4.2 - 来自 100,000 次模拟的结果：在母亲和父亲都是杂合子的情况下，后代在某个基因座上为 AA 或杂合子的数量

在这种情况下，我们也有带重复的排列，但我们有四个可能的值，而不是两个：AA、AT、TA 和 TT。结果是所有个体为 AT 的概率相同：9.5367431640625e-07。对于所有个体都是同型合子的情况（所有为 TT 或者所有为 AA），情况更糟（实际上是两倍糟糕）。

  1. 好的，在这个概率性的序言之后，让我们开始更多数据处理的工作。我们首先要做的是检查我们有多少个错误。让我们从前一个例子中加载数据：

    import gzip
import pickle
import random
import numpy as np
mendelian_errors = pickle.load(gzip.open('mendelian_errors.pickle.gz', 'rb'))
feature_fit = np.load(gzip.open('feature_fit.npy.gz', 'rb'))
ordered_features = np.load(open('ordered_features', 'rb'))
num_features = len(ordered_features)

  2. 让我们看看我们有多少个错误：

    print(len(mendelian_errors), len(list(filter(lambda x: x[0] > 0,mendelian_errors.values()))))

输出如下：

(10905732, 541688)

并不是所有的调用都有孟德尔错误——只有大约 5%，很好。

  1. 让我们创建一个平衡集，其中大约一半的集合有错误。为此，我们将随机丢弃大量正常调用。首先，我们计算错误的比例：

    total_observations = len(mendelian_errors)
error_observations = len(list(filter(lambda x: x[0] > 0,mendelian_errors.values())))
ok_observations = total_observations - error_observations
fraction_errors = error_observations/total_observations
print (total_observations, ok_observations, error_observations, 100*fraction_errors)
del mendelian_errors

  2. 我们使用这些信息来获取一组被接受的条目：所有错误和一个大致相等数量的正常调用。我们在最后打印条目数量（这将随着 OK 列表的随机性而变化）：

    prob_ok_choice = error_observations / ok_observations
def accept_entry(row):
    if row[-1] == 1:
        return True
    return random.random() <= prob_ok_choice
accept_entry_v = np.vectorize(accept_entry, signature='(i)->()')
accepted_entries = accept_entry_v(feature_fit)
balanced_fit = feature_fit[accepted_entries]
del feature_fit
balanced_fit.shape
len([x for x in balanced_fit if x[-1] == 1]), len([x for x in balanced_fit if x[-1] == 0])

  3. 最后，我们保存它：

    np.save(gzip.open('balanced_fit.npy.gz', 'wb'), balanced_fit, allow_pickle=False, fix_imports=False)

还有更多……

关于孟德尔错误及其对成本函数的影响，让我们思考以下情况：母亲是 AA，父亲是 AT，所有后代都是 AA。这是否意味着父亲的基因型判断错误，或者我们未能检测到一些杂合的后代？从这个推理来看，可能是父亲的基因型判断错误。这在一些更精细的孟德尔错误估计函数中有影响：让几个后代错误比仅仅一个样本（父亲）错误可能更有成本。在这种情况下，你可能会认为这很简单（没有杂合子后代的概率很低，所以可能是父亲的错误），但如果有 18 个后代是 AA，2 个是 AT，那是否还能算“简单”呢？这不仅仅是一个理论问题，因为它会严重影响成本函数的设计。

我们在之前的食谱中的函数，为分析准备数据集，虽然很简单，但足够满足我们进一步获得有趣结果所需的精度。

使用标准统计方法探索数据

现在我们有了孟德尔错误分析的见解，让我们探索数据，以便获得更多可能帮助我们更好地过滤数据的见解。你可以在Chapter04/Exploration.py中找到此内容。

如何做到……

  1. 我们像往常一样，先导入必要的库：

    import gzip
import pickle
import random
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pandas.plotting import scatter_matrix

  2. 然后我们加载数据。我们将使用 pandas 来导航：

    fit = np.load(gzip.open('balanced_fit.npy.gz', 'rb'))
ordered_features = np.load(open('ordered_features', 'rb'))
num_features = len(ordered_features)
fit_df = pd.DataFrame(fit, columns=ordered_features + ['pos', 'error'])
num_samples = 80
del fit

  3. 让我们让 pandas 显示所有注释的直方图：

    fig,ax = plt.subplots(figsize=(16,9))
fit_df.hist(column=ordered_features, ax=ax)

生成的直方图如下：

图 4.3 - 数据集中所有注释的直方图，错误约占 50%

图 4.3 - 数据集中所有注释的直方图，错误约占 50%

  1. 对于某些注释，我们没有得到有趣的信息。我们可以尝试放大，举个例子，使用 DP：

    fit_df['MeanDP'] = fit_df['DP'] / 80
fig, ax = plt.subplots()
_ = ax.hist(fit_df[fit_df['MeanDP']<50]['MeanDP'], bins=100)

图 4.4 - 放大显示 DP 相关兴趣区域的直方图

图 4.4 - 放大显示 DP 相关兴趣区域的直方图

实际上，我们将 DP 除以样本数，以便得到一个更有意义的数字。

  1. 我们将把数据集分为两部分，一部分用于错误，另一部分用于没有孟德尔错误的位置：

    errors_df = fit_df[fit_df['error'] == 1]
ok_df = fit_df[fit_df['error'] == 0]

  2. 让我们看一下 QUAL 并以 0.005 为分割点，检查我们如何得到错误和正确调用的分割：

    ok_qual_above_df = ok_df[ok_df['QUAL']>0.005]
errors_qual_above_df = errors_df[errors_df['QUAL']>0.005]
print(ok_df.size, errors_df.size, ok_qual_above_df.size, errors_qual_above_df.size)
print(ok_qual_above_df.size / ok_df.size, errors_qual_above_df.size / errors_df.size)

结果如下：

6507972 6500256 484932 6114096
0.07451353509203788 0.9405931089483245

显然，['QUAL']>0.005产生了很多错误，而没有产生很多正确的位置。这是积极的，因为我们有一些希望通过过滤来处理这些错误。

  1. 对 QD 做同样的处理：

    ok_qd_above_df = ok_df[ok_df['QD']>0.05]
errors_qd_above_df = errors_df[errors_df['QD']>0.05]
print(ok_df.size, errors_df.size, ok_qd_above_df.size, errors_qd_above_df.size)
print(ok_qd_above_df.size / ok_df.size, errors_qd_above_df.size / errors_df.size)

再次，我们得到了有趣的结果：

6507972 6500256 460296 5760288
0.07072802402960554 0.8861632526472804

  1. 让我们选择一个错误较少的区域，研究注释之间的关系。我们将成对绘制注释：

    not_bad_area_errors_df = errors_df[(errors_df['QUAL']<0.005)&(errors_df['QD']<0.05)]
_ = scatter_matrix(not_bad_area_errors_df[['FS', 'ReadPosRankSum', 'MQ', 'HRun']], diagonal='kde', figsize=(16, 9), alpha=0.02)

前面的代码生成了以下输出：

图 4.5 - 搜索空间区域的错误注释散点矩阵

图 4.5 - 搜索空间区域的错误注释散点矩阵

  1. 现在对正确的调用做相同的处理：

    not_bad_area_ok_df = ok_df[(ok_df['QUAL']<0.005)&(ok_df['QD']<0.05)]
_ = scatter_matrix(not_bad_area_ok_df[['FS', 'ReadPosRankSum', 'MQ', 'HRun']], diagonal='kde', figsize=(16, 9), alpha=0.02)

输出结果如下：

图 4.6 - 搜索空间区域内良好标记的注释散点矩阵

图 4.6 - 搜索空间区域内良好标记的注释散点矩阵

  1. 最后，让我们看看我们的规则在完整数据集上的表现（记住，我们使用的数据集大约由 50%的错误和 50%的正确标记组成）：

    all_fit_df = pd.DataFrame(np.load(gzip.open('feature_fit.npy.gz', 'rb')), columns=ordered_features + ['pos', 'error'])
potentially_good_corner_df = all_fit_df[(all_fit_df['QUAL']<0.005)&(all_fit_df['QD']<0.05)]
all_errors_df=all_fit_df[all_fit_df['error'] == 1]
print(len(all_fit_df), len(all_errors_df), len(all_errors_df) / len(all_fit_df))

我们得到如下结果：

10905732 541688 0.04967002673456491

让我们记住，我们的完整数据集中大约有 1090 万个标记，误差大约为 5%。

  1. 让我们获取一些关于我们good_corner的统计数据：

    potentially_good_corner_errors_df = potentially_good_corner_df[potentially_good_corner_df['error'] == 1]
print(len(potentially_good_corner_df), len(potentially_good_corner_errors_df), len(potentially_good_corner_errors_df) / len(potentially_good_corner_df))
print(len(potentially_good_corner_df)/len(all_fit_df))

输出结果如下：

9625754 32180 0.0033431147315836243
0.8826325458942141

所以，我们将误差率从 5%降低到了 0.33%，同时标记数量仅减少到了 960 万个。

还有更多……

从 5%的误差减少到 0.3%，同时失去 12%的标记，这样的变化是好还是坏呢？嗯，这取决于你接下来想要做什么样的分析。也许你的方法能抵御标记丢失，但不太容忍错误，如果是这种情况，这样的变化可能会有帮助。但如果情况相反，也许即使数据集错误更多，你也更倾向于保留完整的数据集。如果你使用不同的方法，可能会根据方法的不同而使用不同的数据集。在这个疟蚊数据集的具体案例中，数据量非常大，因此减少数据集大小对几乎所有情况都没有问题。但如果标记数量较少，你需要根据标记和质量来评估你的需求。

从测序注释中寻找基因组特征

我们将以一个简单的步骤总结这一章及本书内容，表明有时你可以从简单的、意外的结果中学到重要的东西，而表面上的质量问题可能掩盖了重要的生物学问题。

我们将绘制读取深度——DP——在染色体臂 2L 上所有交叉父本的分布情况。此步骤的代码可以在 Chapter04/2L.py 中找到。

如何做……

我们将从以下步骤开始：

  1. 让我们从常规导入开始：

    from collections import defaultdict
import gzip
import numpy as np
import matplotlib.pylab as plt

  2. 让我们加载在第一步中保存的数据：

    num_parents = 8
dp_2L = np.load(gzip.open('DP_2L.npy.gz', 'rb'))
print(dp_2L.shape)

  3. 现在让我们打印整个染色体臂的中位 DP，以及其中部的一部分数据，针对所有父本：

    for i in range(num_parents):
    print(np.median(dp_2L[:,i]), np.median(dp_2L[50000:150000,i]))

输出结果如下：

17.0 14.0
23.0 22.0
31.0 29.0
28.0 24.0
32.0 27.0
31.0 31.0
25.0 24.0
24.0 20.0

有趣的是，整个染色体的中位数有时并不适用于中间的那个大区域，所以我们需要进一步挖掘。

  1. 我们将打印染色体臂上 200,000 kb 窗口的中位 DP。让我们从窗口代码开始：

    window_size = 200000
parent_DP_windows = [defaultdict(list) for i in range(num_parents)]
def insert_in_window(row):
    for parent in range(num_parents):
        parent_DP_windows[parent][row[-1] // window_size].append(row[parent])
insert_in_window_v = np.vectorize(insert_in_window, signature='(n)->()')
_ = insert_in_window_v(dp_2L)

  2. 让我们绘制它：

    fig, axs = plt.subplots(2, num_parents // 2, figsize=(16, 9), sharex=True, sharey=True, squeeze=True)
for parent in range(num_parents):
    ax = axs[parent // 4][parent % 4]
    parent_data = parent_DP_windows[parent]
    ax.set_ylim(10, 40)
    ax.plot(*zip(*[(win*window_size, np.mean(lst)) for win, lst in parent_data.items()]), '.')

  3. 以下是绘制的结果：

图 4.7 - 所有父本数据集在染色体臂 2L 上每个窗口的中位 DP

图 4.7 - 所有父本数据集在染色体臂 2L 上每个窗口的中位 DP

你会注意到，在某些蚊子样本中，例如第一列和最后一列，染色体臂中间有明显的 DP 下降。在某些样本中，比如第三列的样本，下降程度较轻——不那么明显。对于第二列的底部父本样本，根本没有下降。

还有更多……

前述模式有一个生物学原因，最终对测序产生影响：按蚊可能在 2L 臂中间有一个大的染色体倒位。与用于做调用的参考基因组不同的核型，由于进化分化，较难进行调用。这导致该区域的测序读取数量较少。这在这种物种中特别明显，但你可能会期望在其他生物中也出现类似特征。

一个更广为人知的案例是 n，此时你可以预期在整个基因组中看到 n 倍于中位数的 DP。

但在一般情况下，保持警惕并留意分析中的 异常 结果是个好主意。有时候，这正是一个有趣生物学特征的标志，就像这里一样。要么这就是指向一个错误的信号：例如，主成分分析（PCA）可以用来发现标签错误的样本（因为它们可能会聚集在错误的组中）。

使用 QIIME 2 Python API 进行宏基因组学分析

Wikipedia 表示，宏基因组学是直接从环境样本中回收遗传物质的研究。注意，这里的“环境”应广泛理解：在我们的例子中，我们将处理肠道微生物组，研究的是儿童肠道问题中的粪便微生物组移植。该研究是 QIIME 2 的其中一个教程，QIIME 2 是最广泛使用的宏基因组数据分析应用之一。QIIME 2 有多个接口：图形用户界面（GUI）、命令行界面和一个称为 Artifact API 的 Python API。

Tomasz Kościółek 提供了一个出色的 Artifact API 使用教程，基于 QIIME 2 最成熟的（客户端基础的，而非基于 Artifact 的）教程，即 “Moving Pictures” 教程（nbviewer.jupyter.org/gist/tkosciol/29de5198a4be81559a075756c2490fde）。在这里，我们将创建一个 Python 版本的粪便微生物群移植研究，正如客户端接口一样，详细内容可参见 docs.qiime2.org/2022.2/tutorials/fmt/。你应该熟悉它，因为我们不会在这里深入探讨生物学的细节。我走的路线比 Tomasz 更为复杂：这将帮助你更好地了解 QIIME 2 Python 的内部结构。获得这些经验后，你可能更倾向于按照 Tomasz 的路线，而非我的路线。然而，你在这里获得的经验将使你在使用 QIIME 的内部功能时更加舒适和自信。

准备开始

这个设置稍微复杂一些。我们将需要创建一个 conda 环境，将 QIIME 2 的软件包与其他应用程序的软件包分开。你需要遵循的步骤很简单。

在 OS X 上，使用以下代码创建一个新的 conda 环境：

wget wget https://data.qiime2.org/distro/core/qiime2-2022.2-py38-osx-conda.yml
conda env create -n qiime2-2022.2 --file qiime2-2022.2-py38-osx-conda.yml

在 Linux 上，使用以下代码创建环境：

wget wget https://data.qiime2.org/distro/core/qiime2-2022.2-py38-linux-conda.yml
conda env create -n qiime2-2022.2 --file qiime2-2022.2-py38-linux-conda.yml

如果这些指令不起作用，请查看 QIIME 2 网站上的更新版本（docs.qiime2.org/2022.2/install/native）。QIIME 2 会定期更新。

在这个阶段，你需要通过使用source activate qiime2-2022.2进入 QIIME 2 的conda环境。如果你想回到标准的conda环境，可以使用source deactivate。我们将安装jupyter lab和jupytext：

conda install jupyterlab jupytext

你可能想要在 QIIME 2 的环境中使用conda install安装其他你想要的包。

为了准备 Jupyter 执行，你应该安装 QIIME 2 扩展，方法如下：

jupyter serverextension enable --py qiime2 --sys-prefix

提示

该扩展具有高度交互性，允许你从不同的视角查看数据，这些视角在本书中无法捕捉。缺点是它不能在nbviewer中工作（某些单元格的输出在静态查看器中无法显示）。记得与扩展中的输出进行交互，因为许多输出是动态的。

你现在可以启动 Jupyter。Notebook 可以在Chapter4/QIIME2_Metagenomics.py文件中找到。

警告

由于 QIIME 的包安装具有流动性，我们没有为其提供 Docker 环境。这意味着如果你是通过我们的 Docker 安装工作，你将需要下载食谱并手动安装这些包。

你可以找到获取 Notebook 文件和 QIIME 2 教程数据的说明。

如何操作...

让我们看一下接下来的步骤：

  1. 让我们首先检查一下有哪些插件可用：

    import pandas as pd
from qiime2.metadata.metadata import Metadata
from qiime2.metadata.metadata import CategoricalMetadataColumn
from qiime2.sdk import Artifact
from qiime2.sdk import PluginManager
from qiime2.sdk import Result
pm = PluginManager()
demux_plugin = pm.plugins['demux']
#demux_emp_single = demux_plugin.actions['emp_single']
demux_summarize = demux_plugin.actions['summarize']
print(pm.plugins)

我们还在访问解混插件及其摘要操作。

  1. 让我们来看看摘要操作，即inputs、outputs和parameters：

    print(demux_summarize.description)
demux_summarize_signature = demux_summarize.signature
print(demux_summarize_signature.inputs)
print(demux_summarize_signature.parameters)
print(demux_summarize_signature.outputs)

输出将如下所示：

Summarize counts per sample for all samples, and generate interactive positional quality plots based on `n` randomly selected sequences.
 OrderedDict([('data', ParameterSpec(qiime_type=SampleData[JoinedSequencesWithQuality | PairedEndSequencesWithQuality | SequencesWithQuality], view_type=<class 'q2_demux._summarize._visualizer._PlotQualView'>, default=NOVALUE, description='The demultiplexed sequences to be summarized.'))])
 OrderedDict([('n', ParameterSpec(qiime_type=Int, view_type=<class 'int'>, default=10000, description='The number of sequences that should be selected at random for quality score plots. The quality plots will present the average positional qualities across all of the sequences selected. If input sequences are paired end, plots will be generated for both forward and reverse reads for the same `n` sequences.'))])
 OrderedDict([('visualization', ParameterSpec(qiime_type=Visualization, view_type=None, default=NOVALUE, description=NOVALUE))])

  1. 现在，我们将加载第一个数据集，对其进行解混，并可视化一些解混统计数据：

    seqs1 = Result.load('fmt-tutorial-demux-1-10p.qza')
sum_data1 = demux_summarize(seqs1)
sum_data1.visualization

这是来自 QIIME 扩展为 Jupyter 提供的输出的一部分：

图 4.8 - QIIME 2 扩展为 Jupyter 提供的输出部分

图 4.8 - QIIME 2 扩展为 Jupyter 提供的输出部分

记住，扩展是迭代的，并且提供的信息比仅仅这张图表多得多。

提示

本食谱的原始数据是以 QIIME 2 格式提供的。显然，你将拥有自己原始数据的其他格式（可能是 FASTQ 格式）—请参见*还有更多...*部分，了解如何加载标准格式。

QIIME 2 的.qza和.qzv格式只是压缩文件。你可以通过unzip查看其内容。

图表将类似于 QIIME CLI 教程中的图表，但务必检查我们输出的交互质量图。

  1. 让我们对第二个数据集做相同的操作：

    seqs2 = Result.load('fmt-tutorial-demux-2-10p.qza')
sum_data2 = demux_summarize(seqs2)
sum_data2.visualization

  2. 让我们使用 DADA2（github.com/benjjneb/dada2）插件进行质量控制：

    dada2_plugin = pm.plugins['dada2']
dada2_denoise_single = dada2_plugin.actions['denoise_single']
qual_control1 = dada2_denoise_single(demultiplexed_seqs=seqs1,
                                    trunc_len=150, trim_left=13)
qual_control2 = dada2_denoise_single(demultiplexed_seqs=seqs2,
                                    trunc_len=150, trim_left=13)

  3. 让我们从去噪（第一组）中提取一些统计数据：

    metadata_plugin = pm.plugins['metadata']
metadata_tabulate = metadata_plugin.actions['tabulate']
stats_meta1 = metadata_tabulate(input=qual_control1.denoising_stats.view(Metadata))
stats_meta1.visualization

同样，结果可以在 QIIME 2 命令行版本的教程中找到。

  1. 现在，让我们对第二组数据做相同的操作：

    stats_meta2 = metadata_tabulate(input=qual_control2.denoising_stats.view(Metadata))
stats_meta2.visualization

  2. 现在，合并去噪后的数据：

    ft_plugin = pm.plugins['feature-table']
ft_merge = ft_plugin.actions['merge']
ft_merge_seqs = ft_plugin.actions['merge_seqs']
ft_summarize = ft_plugin.actions['summarize']
ft_tab_seqs = ft_plugin.actions['tabulate_seqs']
table_merge = ft_merge(tables=[qual_control1.table, qual_control2.table])
seqs_merge = ft_merge_seqs(data=[qual_control1.representative_sequences, qual_control2.representative_sequences])

  3. 然后，从合并结果中收集一些质量统计数据：

    ft_sum = ft_summarize(table=table_merge.merged_table)
ft_sum.visualization

  4. 最后，让我们获取一些关于合并序列的信息：

    tab_seqs = ft_tab_seqs(data=seqs_merge.merged_data)
tab_seqs.visualization

还有更多内容...

上面的代码没有展示如何导入数据。实际代码会根据情况有所不同（单端数据、双端数据，或已经解多重条形码的数据），但对于主要的 QIIME 2 教程《移动图像》，假设你已经将单端、未解多重条形码的数据和条形码下载到名为data的目录中，你可以执行以下操作：

data_type = 'EMPSingleEndSequences'
conv = Artifact.import_data(data_type, 'data')
conv.save('out.qza')

如上面的代码所述，如果你在 GitHub 上查看这个笔记本，静态的nbviewer系统将无法正确渲染笔记本（你需要自己运行它）。这远非完美；它不具备交互性，质量也不是很好，但至少它能让你在不运行代码的情况下大致了解输出结果。

第六章：处理基因组

计算生物学中的许多任务依赖于参考基因组的存在。如果你正在进行序列比对、寻找基因或研究种群遗传学，你将直接或间接使用参考基因组。在本章中，我们将开发一些处理参考基因组的技术，解决不同质量的参考基因组问题，这些基因组的质量可能从高质量（在这里，高质量仅指基因组组装的状态，这是本章的重点），如人类基因组，到存在问题的非模式物种基因组。我们还将学习如何处理基因组注释（处理能够指引我们发现基因组中有趣特征的数据库），并使用注释信息提取序列数据。我们还将尝试跨物种寻找基因同源物。最后，我们将访问基因本体论（GO）数据库。

在本章中，我们将涵盖以下内容：

  • 处理高质量参考基因组

  • 处理低质量参考基因组

  • 遍历基因组注释

  • 使用注释从参考基因组中提取基因

  • 使用 Ensembl REST API 查找同源基因

  • 从 Ensembl 获取基因本体信息

技术要求

如果你通过 Docker 运行本章内容，你可以使用tiagoantao/bioinformatics_genomes镜像。如果你使用 Anaconda，本章所需的软件将在每个相关部分介绍。

处理高质量参考基因组

在本节中，你将学习一些操作参考基因组的通用技术。作为一个示例，我们将研究恶性疟原虫的 GC 含量——即基因组中以鸟嘌呤-胞嘧啶为基础的部分，这是导致疟疾的最重要寄生虫物种。参考基因组通常以 FASTA 文件的形式提供。

准备工作

生物体基因组的大小差异很大，从像 HIV 这样的病毒（其基因组为 9.7 kbp）到像大肠杆菌这样的细菌，再到像恶性疟原虫这样的原生动物（其基因组跨越 14 条染色体、线粒体和顶体，大小为 22 Mbp），再到果蝇，具有三对常染色体、线粒体和 X/Y 性染色体，再到人类，其基因组由 22 对常染色体、X/Y 染色体和线粒体组成，总大小为 3 Gbp，一直到日本巴黎，一种植物，基因组大小为 150 Gbp。在这个过程中，你会遇到不同的倍性和性染色体组织。

提示

正如你所看到的，不同的生物体有非常不同的基因组大小。这个差异可以达到几个数量级。这对你的编程风格有重大影响。处理一个大型基因组将要求你更加节省内存。不幸的是，更大的基因组将从更高效的编程技术中受益（因为你需要分析的数据更多）；这些是相互矛盾的要求。一般来说，对于较大的基因组，你必须更加小心地处理效率（无论是速度还是内存）。

为了让这个方法不那么繁琐，我们将使用恶性疟原虫的一个小型真核基因组。这个基因组仍然具有大基因组的一些典型特征（例如，多个染色体）。因此，它在复杂性和大小之间是一个很好的折衷。请注意，对于像恶性疟原虫这样大小的基因组，可以通过将整个基因组加载到内存中来执行许多操作。然而，我们选择了一种可以应用于更大基因组（例如，哺乳动物）的编程风格，这样你可以以更通用的方式使用这个方法，但如果是像这种小型基因组，也可以使用更依赖内存的方式。

我们将使用 Biopython，这是你在第一章，Python 与相关软件生态中安装的。如往常一样，这个方法可以在本书的 Jupyter 笔记本中找到，路径为Chapter05/Reference_Genome.py，在本书的代码包中。我们需要下载参考基因组——你可以在上述笔记本中找到最新的下载地址。为了生成本方法最后的图表，我们还需要reportlab：

conda install -c bioconda reportlab

现在，我们准备开始。

如何操作...

按照以下步骤进行：

  1. 我们将首先检查参考基因组 FASTA 文件中所有序列的描述：

    from Bio import SeqIO
genome_name = 'PlasmoDB-9.3_Pfalciparum3D7_Genome.fasta'
recs = SeqIO.parse(genome_name, 'fasta')
for rec in recs:
    print(rec.description)

这段代码应该很熟悉，来自上一章，第三章，下一代测序。让我们来看一下部分输出：

图 5.1 – 显示恶性疟原虫参考基因组的 FASTA 描述的输出

不同的基因组参考将有不同的描述行，但它们通常包含重要信息。在这个示例中，你可以看到我们有染色体、线粒体和顶体。我们还可以查看染色体的大小，但我们将使用序列长度中的值。

  1. 让我们解析描述行，以提取染色体编号。我们将从序列中获取染色体大小，并基于窗口计算每个染色体的GC含量：

    from Bio import SeqUtils
recs = SeqIO.parse(genome_name, 'fasta')
chrom_sizes = {}
chrom_GC = {}
block_size = 50000
min_GC = 100.0
max_GC = 0.0
for rec in recs:
    if rec.description.find('SO=chromosome') == -1:
        continue
    chrom = int(rec.description.split('_')[1])
    chrom_GC[chrom] = []
    size = len(rec.seq)
    chrom_sizes[chrom] = size
    num_blocks = size // block_size + 1
    for block in range(num_blocks):
        start = block_size * block
        if block == num_blocks - 1:
            end = size
        else:
            end = block_size + start + 1
        block_seq = rec.seq[start:end]
        block_GC = SeqUtils.GC(block_seq)
        if block_GC < min_GC:
            min_GC = block_GC
        if block_GC > max_GC:
            max_GC = block_GC
        chrom_GC[chrom].append(block_GC)
print(min_GC, max_GC)

在这里，我们对所有染色体进行了窗口分析，类似于我们在第三章中所做的，下一代测序。我们从定义一个 50 kbp 的窗口大小开始。这对于Plasmodium falciparum来说是合适的（你可以自由调整其大小），但对于那些染色体大小差异数量级较大的基因组，你会想考虑其他的数值。

注意，我们正在重新读取文件。由于基因组如此之小，实际上在步骤 1中可以将整个基因组加载到内存中。对于小型基因组来说，尝试这种编程风格是可行的——它更快！然而，我们的代码是为了可以在更大的基因组上复用而设计的。

  1. 注意，在for循环中，我们通过解析SO条目来忽略线粒体和顶体。chrom_sizes字典将维护染色体的大小。

chrom_GC字典是我们最有趣的数据结构，将包含每个 50 kbp 窗口中GC含量的一个分数的列表。所以，染色体 1 的大小为 640,851 bp，它会有 14 个条目，因为该染色体的大小为 14 个 50 kbp 的块。

注意疟原虫（Plasmodium falciparum）基因组的两个不寻常特征：该基因组非常富含 AT，即 GC 贫乏。因此，你得到的数字会非常低。此外，染色体是按大小排序的（这很常见），但从最小的大小开始。通常的约定是从最大的大小开始（比如在人类基因组中）。

  1. 现在，让我们创建一个GC分布的基因组图。我们将使用蓝色的不同色调表示GC含量。然而，对于高离群值，我们将使用红色的不同色调。对于低离群值，我们将使用黄色的不同色调：

    from reportlab.lib import colors
from reportlab.lib.units import cm
from Bio.Graphics import BasicChromosome
chroms = list(chrom_sizes.keys())
chroms.sort()
biggest_chrom = max(chrom_sizes.values())
my_genome = BasicChromosome.Organism(output_format="png")
my_genome.page_size = (29.7*cm, 21*cm)
telomere_length = 10
bottom_GC = 17.5
top_GC = 22.0
for chrom in chroms:
    chrom_size = chrom_sizes[chrom]
    chrom_representation = BasicChromosome.Chromosome ('Cr %d' % chrom)
    chrom_representation.scale_num = biggest_chrom
    tel = BasicChromosome.TelomereSegment()
    tel.scale = telomere_length
    chrom_representation.add(tel)
    num_blocks = len(chrom_GC[chrom])
    for block, gc in enumerate(chrom_GC[chrom]):
        my_GC = chrom_GC[chrom][block]
        body = BasicChromosome.ChromosomeSegment()
        if my_GC > top_GC:
            body.fill_color = colors.Color(1, 0, 0)
        elif my_GC < bottom_GC:
            body.fill_color = colors.Color(1, 1, 0)
        else:
            my_color = (my_GC - bottom_GC) / (top_GC -bottom_GC)
            body.fill_color = colors.Color(my_color,my_color, 1)
        if block < num_blocks - 1:
            body.scale = block_size
        else:
            body.scale = chrom_size % block_size
        chrom_representation.add(body)
    tel = BasicChromosome.TelomereSegment(inverted=True)
    tel.scale = telomere_length
    chrom_representation.add(tel)
    my_genome.add(chrom_representation)
my_genome.draw('falciparum.png', 'Plasmodium falciparum')

第一行将keys方法的返回值转换为一个列表。在 Python 2 中这是多余的，但在 Python 3 中并非如此，因为keys方法的返回类型是特定的dict_keys类型。

我们按顺序绘制染色体（因此需要排序）。我们需要最大的染色体的大小（在Plasmodium falciparum中为 14）来确保染色体的大小以正确的比例打印（即biggest_chrom变量）。

然后，我们创建一个 A4 大小的有机体表示，并输出为 PNG 文件。注意，我们绘制了非常小的端粒（10 bp）。这将产生一个类似矩形的染色体。你可以将端粒做得更大，给它们一个圆形的表示，或者你也可以选择使用适合你物种的端粒尺寸。

我们声明，任何GC含量低于 17.5%或高于 22.0%的都会被视为离群值。记住，对于大多数其他物种来说，这个值会更高。

然后，我们打印这些染色体：它们被端粒限制，且由 50 kbp 的染色体片段组成（最后一个片段的大小为剩余部分）。每个片段将用蓝色表示，并基于两种极值之间的线性归一化，呈现红绿成分。每个染色体片段将为 50 kbp，或者如果它是染色体的最后一个片段，可能会更小。输出结果如下图所示：

图 5.2 – 疟原虫的 14 条染色体，用 GC 含量进行颜色编码（红色表示超过 22%，黄色表示少于 17%，蓝色阴影代表这两个数值之间的线性渐变）

提示

Biopython 代码的演变发生在 Python 成为流行语言之前。过去，库的可用性非常有限。reportlab的使用大多可以视为遗留问题。我建议你只需学到足够的知识，能与 Biopython 配合使用。如果你计划学习现代的 Python 绘图库，那么标准的代表是 Matplotlib，正如我们在*第二章**中所学的，《了解 NumPy、pandas、Arrow 和 Matplotlib》。替代方案包括 Bokeh、HoloViews，或 Python 版本的 ggplot（甚至更复杂的可视化替代方案，如 Mayavi、可视化工具包（VTK）甚至 Blender API）。

  1. 最后，你可以在笔记本中内嵌打印图像：

    from IPython.core.display import Image
Image("falciparum.png")

就这样完成了这个方案！

还有更多……

疟原虫是一个合理的例子，它是一个基因组较小的真核生物，能够让你进行一个具有足够特征的小型数据练习，同时对于大多数真核生物仍然具有参考价值。当然，它没有性染色体（如人类的 X/Y 染色体），但这些应该容易处理，因为参考基因组并不涉及倍性问题。

疟原虫确实有线粒体，但由于篇幅限制，我们在此不讨论它。Biopython 确实具有打印圆形基因组的功能，你也可以将其用于细菌。关于细菌和病毒，这些基因组更容易处理，因为它们的大小非常小。

另见

这里有一些你可以深入了解的资源：

  • 你可以在 Ensembl 网站上找到许多模式生物的参考基因组，网址是www.ensembl.org/info/data/ftp/index.xhtml。

  • 像往常一样，国家生物技术信息中心（NCBI）也提供了一个庞大的基因组列表，网址是www.ncbi.nlm.nih.gov/genome/browse/。

  • 有许多网站专注于单一生物体（或一组相关生物体）。除了你从Plasmodium falciparum基因组下载的 PlasmoDB（plasmodb.org/plasmo/），你还会在下一个关于病媒生物的食谱中找到 VectorBase（www.vectorbase.org/）。用于*果蝇（Drosophila melanogaster）*的 FlyBase（flybase.org/）也值得一提，但不要忘记搜索你感兴趣的生物体。

处理低质量的基因组参考

不幸的是，并非所有参考基因组都具备Plasmodium falciparum那样的质量。除了某些模型物种（例如人类，或常见的果蝇Drosophila melanogaster）和少数其他物种外，大多数参考基因组仍有待改进。在本食谱中，我们将学习如何处理低质量的参考基因组。

准备工作

继续围绕疟疾主题，我们将使用两种疟疾传播蚊子的参考基因组：Anopheles gambiae（这是最重要的疟疾传播者，存在于撒哈拉以南的非洲）和Anopheles atroparvus，一种欧洲的疟疾传播者（尽管欧洲已经消除了该病，但该传播者依然存在）。Anopheles gambiae基因组质量较好，大多数染色体已经被映射，尽管 Y 染色体仍需要一些工作。还有一个相当大的未知染色体，可能由 X 和 Y 染色体的部分片段以及中肠微生物群组成。该基因组中有一些位置未被注释（即，你会看到N而不是 ACTG）。Anopheles atroparvus基因组仍处于框架格式。遗憾的是，许多非模型物种的基因组都是这种情况。

请注意，我们将稍微提高难度。Anopheles基因组比Plasmodium falciparum基因组大一个数量级（但仍然比大多数哺乳动物小一个数量级）。

我们将使用你在第一章中安装的 Biopython，Python 和周边软件生态。像往常一样，本食谱可在本书的 Jupyter 笔记本中找到，路径为Chapter05/Low_Quality.py，在本书的代码包中。在笔记本的开始部分，你可以找到两个基因组的最新位置，以及下载它们的代码。

如何操作...

请按照以下步骤操作：

  1. 让我们从列出Anopheles gambiae基因组的染色体开始：

    import gzip
from Bio import SeqIO
gambiae_name = 'gambiae.fa.gz'
atroparvus_name = 'atroparvus.fa.gz'
recs = SeqIO.parse(gzip.open(gambiae_name, 'rt', encoding='utf-8'), 'fasta')
for rec in recs:
    print(rec.description)

这将产生一个输出，其中包括生物体的染色体（以及一些未映射的超级重组片段，未显示）：

AgamP4_2L | organism=Anopheles_gambiae_PEST | version=AgamP4 | length=49364325 | SO=chromosome
AgamP4_2R | organism=Anopheles_gambiae_PEST | version=AgamP4 | length=61545105 | SO=chromosome
AgamP4_3L | organism=Anopheles_gambiae_PEST | version=AgamP4 | length=41963435 | SO=chromosome
AgamP4_3R | organism=Anopheles_gambiae_PEST | version=AgamP4 | length=53200684 | SO=chromosome
AgamP4_X | organism=Anopheles_gambiae_PEST | version=AgamP4 | length=24393108 | SO=chromosome
AgamP4_Y_unplaced | organism=Anopheles_gambiae_PEST | version=AgamP4 | length=237045 | SO=chromosome
AgamP4_Mt | organism=Anopheles_gambiae_PEST | version=AgamP4 | length=15363 | SO=mitochondrial_chromosome

代码非常简单。我们使用gzip模块，因为较大基因组的文件通常是压缩的。我们可以看到四个染色体臂（2L、2R、3L 和 3R）、线粒体（Mt）、X 染色体和 Y 染色体，后者非常小，名字几乎表明它的状态可能不佳。此外，未知的（UNKN）染色体占参考基因组的较大比例，几乎相当于一个染色体臂的大小。

不要在Anopheles atroparvus上执行此操作；否则，您将得到超过一千个条目，这归功于 scaffolding 状态。

  1. 现在，让我们检查一下未调用的位置（Ns）及其在*按蚊（Anopheles gambiae）*基因组中的分布：

    recs = SeqIO.parse(gzip.open(gambiae_name, 'rt', encoding='utf-8'), 'fasta')
chrom_Ns = {}
chrom_sizes = {}
for rec in recs:
    if rec.description.find('supercontig') > -1:
        continue
    print(rec.description, rec.id, rec)
    chrom = rec.id.split('_')[1]
    if chrom in ['UNKN']:
        continue
    chrom_Ns[chrom] = []
    on_N = False
    curr_size = 0
    for pos, nuc in enumerate(rec.seq):
        if nuc in ['N', 'n']:
            curr_size += 1
            on_N = True
        else:
            if on_N:
                chrom_Ns[chrom].append(curr_size)
                curr_size = 0
            on_N = False
    if on_N:
        chrom_Ns[chrom].append(curr_size)
    chrom_sizes[chrom] = len(rec.seq)
for chrom, Ns in chrom_Ns.items():
    size = chrom_sizes[chrom]
    if len(Ns) > 0:
        max_Ns = max(Ns)
    else:
        max_Ns = 'NA'
    print(f'{chrom} ({size}): %Ns ({round(100 * sum(Ns) / size, 1)}), num Ns: {len(Ns)}, max N: {max_Ns}')

上面的代码将需要一些时间来运行，请耐心等待；我们将检查所有常染色体的每个碱基对。和往常一样，我们将重新打开并重新读取文件以节省内存。

我们有两个字典：一个包含染色体大小，另一个包含 Ns 运行的大小分布。为了计算 Ns 运行，我们必须遍历所有常染色体（注意何时 N 位置开始和结束）。然后，我们必须打印 Ns 分布的基本统计信息：

2L (49364325): %Ns (1.7), num Ns: 957, max N: 28884
2R (61545105): %Ns (2.3), num Ns: 1658, max N: 36427
3L (41963435): %Ns (2.9), num Ns: 1272, max N: 31063
3R (53200684): %Ns (1.8), num Ns: 1128, max N: 24292
X (24393108): %Ns (4.1), num Ns: 1287, max N: 21132
Y (237045): %Ns (43.0), num Ns: 63, max N: 7957
Mt (15363): %Ns (0.0), num Ns: 0, max N: NA

因此，对于 2L 染色体臂（大小为 49 Mbp），1.7% 是 N 调用，并且分布在 957 个运行中。最大的运行是 28884 bp。请注意，X 染色体具有最高的 Ns 位置比例。

  1. 现在，让我们将注意力转向Anopheles Atroparvus基因组。让我们计算一下 scaffolds 的数量，并且查看 scaffold 大小的分布：

    import numpy as np
recs = SeqIO.parse(gzip.open(atroparvus_name, 'rt', encoding='utf-8'), 'fasta')
sizes = []
size_N = []
for rec in recs:
    size = len(rec.seq)
    sizes.append(size)
    count_N = 0
    for nuc in rec.seq:
        if nuc in ['n', 'N']:
            count_N += 1
    size_N.append((size, count_N / size))
print(len(sizes), np.median(sizes), np.mean(sizes),
      max(sizes), min(sizes),
      np.percentile(sizes, 10), np.percentile(sizes, 90))

这段代码与我们之前看到的相似，但我们使用 NumPy 打印了更详细的统计信息，因此我们得到了以下结果：

1320 7811.5 170678.2 58369459 1004 1537.1 39644.7

因此，我们有 1371 个 scaffolds（与Anopheles gambiae基因组中的七个条目相比），中位大小为 7811.5（平均值为 17,0678.2）。最大的 scaffold 为 5.8 Mbp，最小的 scaffold 为 1004 bp。大小的第十百分位为 1537.1，而第九十百分位为 39644.7。

  1. 最后，让我们绘制 scaffold 的比例 —— 即 N —— 作为其大小的函数：

    import matplotlib.pyplot as plt
small_split = 4800
large_split = 540000
fig, axs = plt.subplots(1, 3, figsize=(16, 9), squeeze=False, sharey=True)
xs, ys = zip(*[(x, 100 * y) for x, y in size_N if x <= small_split])
axs[0, 0].plot(xs, ys, '.')
xs, ys = zip(*[(x, 100 * y) for x, y in size_N if x > small_split and x <= large_split])
axs[0, 1].plot(xs, ys, '.')
axs[0, 1].set_xlim(small_split, large_split)
xs, ys = zip(*[(x, 100 * y) for x, y in size_N if x > large_split])
axs[0, 2].plot(xs, ys, '.')
axs[0, 0].set_ylabel('Fraction of Ns', fontsize=12)
axs[0, 1].set_xlabel('Contig size', fontsize=12)
fig.suptitle('Fraction of Ns per contig size', fontsize=26)

上面的代码将生成如下图所示的输出，在该图中，我们将图表按 scaffold 大小分成三部分：一部分用于小于 4,800 bp 的 scaffolds，一部分用于介于 4,800 和 540,000 bp 之间的 scaffolds，另一部分用于更大的 scaffolds。小型 scaffolds 的 Ns 比例非常低（始终低于 3.5%）；中型 scaffolds 的 Ns 比例有较大变异（大小范围为 0% 至超过 90%）；而对于最大的 scaffolds，Ns 的变异性较小（在 0% 至 25% 之间）。

图 5.3 – 作为其大小函数的 N 片段比例

还有更多内容...

有时，参考基因组携带额外的信息。例如，按蚊基因组是软屏蔽的。这意味着对基因组进行了某些操作，以识别低复杂度区域（这些区域通常更难分析）。这种情况可以通过大写字母注释：ACTG 表示高复杂度，而 actg 表示低复杂度。

拥有大量 scaffolds 的参考基因组不仅仅是麻烦。例如，非常小的 scaffold（比如小于 2000 bp）在使用比对工具时（如Burrows-Wheeler 比对器（BWA））可能会出现比对问题，尤其是在极端位置（大多数 scaffold 在极端位置会有比对问题，但如果 scaffold 较小，这些问题将占据 scaffold 更大的比例）。如果你使用这样的参考基因组进行比对，建议在比对到小 scaffold 时考虑忽略配对信息（假设你有双端读取），或者至少衡量 scaffold 大小对比对工具性能的影响。无论如何，关键在于你应该小心，因为 scaffold 的大小和数量会时不时地带来麻烦。

对于这些基因组，仅识别出完全模糊（N）。请注意，其他基因组组装可能会给出一个介于模糊和确定性（ACTG）之间的中间代码。

另见

以下是一些你可以从中了解更多信息的资源：

遍历基因组注释

拥有一个基因组序列很有趣，但我们还需要从中提取特征，例如基因、外显子和编码序列。这类注释信息通常以通用特征格式（GFF）和通用转移格式（GTF）文件的形式提供。在本教程中，我们将以按蚊基因组的注释为例，学习如何解析和分析 GFF 文件。

准备工作

使用提供的Chapter05/Annotations.py笔记本文件，该文件包含在本书的代码包中。我们将使用的 GFF 文件的最新位置可以在笔记本顶部找到。

你需要安装gffutils：

conda install -c bioconda gffutils

现在，我们准备开始了。

如何做...

按照以下步骤操作：

  1. 让我们从创建一个基于 GFF 文件的注释数据库开始，使用gffutils：

    import gffutils
import sqlite3
try:
    db = gffutils.create_db('gambiae.gff.gz', 'ag.db')
except sqlite3.OperationalError:
    db = gffutils.FeatureDB('ag.db')

gffutils库创建一个 SQLite 数据库来高效存储注释信息。在这里，我们将尝试创建数据库，如果数据库已存在，则使用现有的。这一步可能需要一些时间。

  1. 现在，让我们列出所有可用的特征类型并统计它们：

    print(list(db.featuretypes()))
for feat_type in db.featuretypes():
    print(feat_type, db.count_features_of_type(feat_type))

这些特征包括 contigs、基因、外显子、转录本等等。请注意，我们将使用 gffutils 包的 featuretypes 函数。它将返回一个生成器，但我们会将其转换为列表（在这里这样做是安全的）。

  1. 让我们列出所有的 seqids：

    seqids = set()
for e in db.all_features():
    seqids.add(e.seqid)
for seqid in seqids:
    print(seqid)

这将显示所有染色体臂和性染色体、线粒体以及未知染色体的注释信息。

  1. 现在，让我们按染色体提取大量有用的信息，比如基因数量、每个基因的转录本数量、外显子数量等等：

    from collections import defaultdict
num_mRNAs = defaultdict(int)
num_exons = defaultdict(int)
max_exons = 0
max_span = 0
for seqid in seqids:
    cnt = 0
    for gene in db.region(seqid=seqid, featuretype='protein_coding_gene'):
        cnt += 1
        span = abs(gene.start - gene.end) # strand
        if span > max_span:
            max_span = span
            max_span_gene = gene
        my_mRNAs = list(db.children(gene, featuretype='mRNA'))
        num_mRNAs[len(my_mRNAs)] += 1
        if len(my_mRNAs) == 0:
            exon_check = [gene]
        else:
            exon_check = my_mRNAs
        for check in exon_check:
            my_exons = list(db.children(check, featuretype='exon'))
            num_exons[len(my_exons)] += 1
            if len(my_exons) > max_exons:
                max_exons = len(my_exons)
                max_exons_gene = gene
    print(f'seqid {seqid}, number of genes {cnt}')
print('Max number of exons: %s (%d)' % (max_exons_gene.id, max_exons))
print('Max span: %s (%d)' % (max_span_gene.id, max_span))
print(num_mRNAs)
print(num_exons)

我们将在提取所有蛋白质编码基因的同时遍历所有 seqids（使用 region）。在每个基因中，我们统计可变转录本的数量。如果没有（注意这可能是注释问题，而不是生物学问题），我们统计外显子数（children）。如果有多个转录本，我们统计每个转录本的外显子数。我们还会考虑跨度大小，以检查跨越最大区域的基因。

我们遵循类似的步骤来查找基因和最大的外显子数。最后，我们打印一个字典，包含每个基因的可变转录本数量分布（num_mRNAs）和每个转录本的外显子数量分布（num_exons）。

还有更多...

GFF/GTF 格式有许多变种。不同的 GFF 版本和许多非官方的变体。如果可能的话，选择 GFF 版本 3。然而，残酷的事实是，你会发现处理这些文件非常困难。gffutils 库尽可能地适应了这一点。事实上，许多关于这个库的文档都是帮助你处理各种尴尬变体的（参考 pythonhosted.org/gffutils/examples.xhtml）。

使用 gffutils 也有替代方案（无论是因为你的 GFF 文件有问题，还是你不喜欢这个库的接口或它依赖 SQL 后端）。自己手动解析文件。如果你看一下格式，你会发现它并不复杂。如果你只进行一次性操作，手动解析或许足够。当然，长期来看，一次性操作往往并不是最好的选择。

另外，请注意，注释的质量往往差异很大。随着质量的提高，复杂性也会增加。只要看看人类注释，就能看到这一点的例子。可以预见，随着我们对生物体的认识不断深入，注释的质量和复杂性也会逐渐提升。

另请参见

以下是一些你可以学习更多资源：

从参考基因组中提取基因信息

在这个步骤中，我们将学习如何借助注释文件提取基因序列，并将其与参考 FASTA 文件的坐标对齐。我们将使用 Anopheles gambiae 基因组及其注释文件（如前两个步骤所示）。首先，我们将提取 电压门控钠通道（VGSC）基因，该基因与抗虫剂的抗性相关。

准备开始

如果你已经按照前两个步骤操作，你就准备好了。如果没有，下载 Anopheles gambiae 的 FASTA 文件和 GTF 文件。你还需要准备 gffutils 数据库：

import gffutils
import sqlite3
try:
    db = gffutils.create_db('gambiae.gff.gz', 'ag.db')
except sqlite3.OperationalError:
    db = gffutils.FeatureDB('ag.db')

和往常一样，你将会在 Chapter05/Getting_Gene.py 笔记本文件中找到所有这些内容。

如何操作...

按照以下步骤操作：

  1. 让我们从获取基因的注释信息开始：

    import gzip
from Bio import Seq, SeqIO
gene_id = 'AGAP004707'
gene = db[gene_id]
print(gene)
print(gene.seqid, gene.strand)

gene_id 是从 VectorBase 获取的，它是一个专门用于疾病传播媒介基因组学的在线数据库。对于其他特定情况，你需要知道你的基因 ID（它将依赖于物种和数据库）。输出将如下所示：

AgamP4_2L       VEuPathDB       protein_coding_gene     2358158 2431617 .       +       .       ID=AGAP004707;Name=para;description=voltage-gated sodium channel
AgamP4_2L +

注意该基因位于 2L 染色体臂，并且是以正向方向编码的（+ 链）。

  1. 让我们将 2L 染色体臂的序列保存在内存中（它只有一个染色体，所以我们可以稍微放松一些）：

    recs = SeqIO.parse(gzip.open('gambiae.fa.gz', 'rt', encoding='utf-8'), 'fasta')
for rec in recs:
    print(rec.description)
    if rec.id == gene.seqid:
        my_seq = rec.seq
        break

输出将如下所示：

AgamP4_2L | organism=Anopheles_gambiae_PEST | version=AgamP4 | length=49364325 | SO=chromosome

  1. 让我们创建一个函数来为一系列 CDS 构建基因序列：

    def get_sequence(chrom_seq, CDSs, strand):
    seq = Seq.Seq('')
    for CDS in CDSs:
        my_cds = Seq.Seq(str(my_seq[CDS.start - 1:CDS.end]))
        seq += my_cds
    return seq if strand == '+' else seq.reverse_complement()

这个函数将接收一个染色体序列（在我们的案例中是 2L 臂），一个编码序列的列表（从注释文件中提取），以及链的信息。

我们必须非常小心序列的起始和结束（注意 GFF 文件是基于 1 的，而 Python 数组是基于 0 的）。最后，如果链是负向的，我们将返回反向互补序列。

  1. 尽管我们已经得到了 gene_id，但我们只需要选择这个基因的三个转录本中的一个，因此需要选择一个：

    mRNAs = db.children(gene, featuretype='mRNA')
for mRNA in mRNAs:
    print(mRNA.id)
    if mRNA.id.endswith('RA'):
        break

  2. 现在，让我们获取转录本的编码序列，然后获取基因序列，并进行翻译：

    CDSs = db.children(mRNA, featuretype='CDS', order_by='start')
gene_seq = get_sequence(my_seq, CDSs, gene.strand)
print(len(gene_seq), gene_seq)
prot = gene_seq.translate()
print(len(prot), prot)

  3. 让我们获取以负链方向编码的基因。我们将提取 VGSC 旁边的基因（恰好是负链）。

    reverse_transcript_id = 'AGAP004708-RA'
reverse_CDSs = db.children(reverse_transcript_id, featuretype='CDS', order_by='start')
reverse_seq = get_sequence(my_seq, reverse_CDSs, '-')
print(len(reverse_seq), reverse_seq)
reverse_prot = reverse_seq.translate()
print(len(reverse_prot), reverse_prot)

在这里，我避免了获取基因的所有信息，只是硬编码了转录本 ID。关键是你需要确保你的代码无论在什么链上都能正常工作。

还有更多...

这是一个简单的步骤，涵盖了本章和 第三章，《下一代测序》中介绍的几个概念。虽然它在概念上很简单，但不幸的是充满了陷阱。

提示

使用不同的数据库时，确保基因组组装版本是同步的。使用不同版本可能会导致严重且潜在的隐性错误。请记住，不同版本（至少在主版本号上）有不同的坐标。例如，人在基因组 36 版本中 3 号染色体上的位置 1,234 可能与基因组 38 版本中的 1,234 不同，可能指向不同的 SNP。在人类数据中，你可能会发现很多芯片使用的是基因组 36 版本，而整个基因组序列使用的是基因组 37 版本，而最新的人类组装版本是基因组 38 版本。对于我们的Anopheles示例，你将会看到 3 和 4 版本。大多数物种都会遇到这种情况，所以请注意！

还有一个问题是 Python 中的 0 索引数组与 1 索引的基因组数据库之间的差异。不过，需要注意的是，一些基因组数据库可能也使用 0 索引。

这里还有两个容易混淆的点：转录本与基因选择，就像在更丰富的注释数据库中一样。在这里，你将有几个备选的转录本（如果你想查看一个复杂到让人迷惑的数据库，可以参考人类注释数据库）。另外，标记为exon的字段包含的信息比编码序列要多。为了这个目的，你将需要 CDS 字段。

最后，还有一个链条问题，你将需要基于反向互补进行翻译。

另见

以下是一些资源，你可以从中获取更多信息：

使用 Ensembl REST API 查找同源基因

在这个食谱中，我们将学习如何为某个基因寻找同源基因。这个简单的食谱不仅介绍了同源性检索，还教你如何使用网页上的 REST API 来访问生物数据。最后，虽然不是最重要的，它将作为如何使用编程 API 访问 Ensembl 数据库的入门教程。

在我们的示例中，我们将尝试为人类horse基因组寻找任何同源基因。

准备工作

这个食谱不需要任何预先下载的数据，但由于我们使用的是 Web API，因此需要互联网访问。传输的数据量将受到限制。

我们还将使用requests库来访问 Ensembl。请求 API 是一个易于使用的 Web 请求封装。自然，你也可以使用标准的 Python 库，但这些要麻烦得多。

像往常一样，你可以在 Chapter05/Orthology.py 笔记本文件中找到这些内容。

如何操作...

按照以下步骤操作：

  1. 我们将从创建一个支持函数开始，以执行网络请求：

    import requests
ensembl_server = 'http://rest.ensembl.org'
def do_request(server, service, *args, **kwargs):
    url_params = ''
    for a in args:
        if a is not None:
            url_params += '/' + a
    req = requests.get('%s/%s%s' % (server, service, url_params), params=kwargs, headers={'Content-Type': 'application/json'})
    if not req.ok:
        req.raise_for_status()
    return req.json()

我们首先导入 requests 库并指定根 URL。然后，我们创建一个简单的函数，传入要调用的功能（参见以下示例），并生成完整的 URL。它还会添加可选参数，并指定负载类型为 JSON（这样就能获取默认的 JSON 响应）。它将返回 JSON 格式的响应。通常这是一个嵌套的 Python 数据结构，包含列表和字典。

  1. 接下来，我们将检查服务器上所有可用的物种，写这本书时大约有 110 种物种：

    answer = do_request(ensembl_server, 'info/species')
for i, sp in enumerate(answer['species']):
    print(i, sp['name'])

请注意，这将构造一个以 http://rest.ensembl.org/info/species 为前缀的 URL，用于 REST 请求。顺便说一下，前面的链接在你的浏览器中无法使用，它应该仅通过 REST API 使用。

  1. 现在，让我们尝试查找与人类数据相关的 HGNC 数据库：

    ext_dbs = do_request(ensembl_server, 'info/external_dbs', 'homo_sapiens', filter='HGNC%')
print(ext_dbs)

我们将搜索限制为与人类相关的数据库（homo_sapiens）。我们还会筛选以 HGNC 开头的数据库（这个筛选使用 SQL 表达式）。HGNC 是 HUGO 数据库。我们要确保它可用，因为 HUGO 数据库负责整理人类基因名称并维护我们的 LCT 标识符。

  1. 现在我们知道 LCT 标识符可能可用，我们希望检索该基因的 Ensembl ID，如以下代码所示：

    answer = do_request(ensembl_server, 'lookup/symbol', 'homo_sapiens', 'LCT')
print(answer)
lct_id = answer['id']

提示

正如你现在可能知道的，不同的数据库会为相同的对象分配不同的 ID。我们需要将我们的 LCT 标识符解析为 Ensembl ID。当你处理与相同对象相关的外部数据库时，数据库之间的 ID 转换可能是你需要完成的第一项任务。

  1. 仅供参考，我们现在可以获取包含基因的区域的序列。请注意，这可能是整个区间，因此如果你想恢复基因，你需要使用类似于我们在前一个步骤中使用的方法：

    lct_seq = do_request(ensembl_server, 'sequence/id', lct_id)
print(lct_seq)

  2. 我们还可以检查 Ensembl 已知的其他数据库；参见以下基因：

    lct_xrefs = do_request(ensembl_server, 'xrefs/id', lct_id)
for xref in lct_xrefs:
    print(xref['db_display_name'])
    print(xref)

你会发现不同种类的数据库，比如 脊椎动物基因组注释 (Vega) 项目、UniProt（参见 第八章，使用蛋白质数据库）和 WikiGene。

  1. 让我们获取这个基因在 horse 基因组上的同源基因：

    hom_response = do_request(ensembl_server, 'homology/id', lct_id, type='orthologues', sequence='none')
homologies = hom_response['data'][0]['homologies']
for homology in homologies:
    print(homology['target']['species'])
    if homology['target']['species'] != 'equus_caballus':
        continue
    print(homology)
    print(homology['taxonomy_level'])
    horse_id = homology['target']['id']

我们本可以通过在 do_request 中指定 target_species 参数，直接获取 horse 基因组的同源基因。然而，这段代码允许你检查所有可用的同源基因。

你将获得关于同源基因的许多信息，例如同源性分类的分类学级别（Boreoeutheria—有胎盘哺乳动物是人类与马匹之间最近的系统发育级别）、同源基因的 Ensembl ID、dN/dS 比率（非同义突变与同义突变的比例），以及 CIGAR 字符串（请参见前一章，第三章，下一代测序）来表示序列之间的差异。默认情况下，你还会得到同源序列的比对结果，但为了简化输出，我已经将其移除。

  1. 最后，让我们查找horse_id的 Ensembl 记录：

    horse_req = do_request(ensembl_server, 'lookup/id', horse_id)
print(horse_req)

从这一点开始，你可以使用之前配方中的方法来探索 LCT horse同源基因。

还有更多…

你可以在rest.ensembl.org/找到所有可用功能的详细解释。这包括所有接口和 Python 代码片段等语言。

如果你对同源基因感兴趣，可以通过前面的步骤轻松地从之前的配方中获取这些信息。在调用homology/id时，只需将类型替换为paralogues。

如果你听说过 Ensembl，那么你可能也听说过 UCSC 的一个替代服务：基因组浏览器(genome.ucsc.edu/)。从用户界面的角度来看，它们在同一层级。从编程的角度来看，Ensembl 可能更加成熟。访问 NCBI Entrez 数据库在第三章，下一代测序中有介绍。

另一种完全不同的编程接口方式是下载原始数据表并将其导入到本地 MySQL 数据库中。请注意，这本身会是一个相当大的工程（你可能只想加载非常小的一部分数据表）。然而，如果你打算进行非常密集的使用，可能需要考虑创建部分数据库的本地版本。如果是这种情况，你可能需要重新考虑 UCSC 的替代方案，因为从本地数据库的角度来看，它和 Ensembl 一样优秀。

从 Ensembl 获取基因本体论信息

在本步骤中，你将再次学习如何通过查询 Ensembl REST API 来使用基因本体论信息。基因本体论是用于注释基因及基因产物的受控词汇。这些词汇以概念树的形式提供（越通用的概念在层次结构的顶部）。基因本体论有三个领域：细胞成分、分子功能和生物过程。

准备工作

与之前的步骤一样，我们不需要任何预下载的数据，但由于我们使用的是 Web API，因此需要互联网访问。传输的数据量将是有限的。

如常，你可以在Chapter05/Gene_Ontology.py笔记本文件中找到这些内容。我们将使用在前一部分(使用 Ensembl REST API 查找直系同源基因)中定义的do_request函数。为了绘制 GO 树，我们将使用pygraphviz，一个图形绘制库：

conda install pygraphviz

好的，我们准备好了。

如何操作...

按照以下步骤操作：

  1. 让我们从检索与 LCT 基因相关的所有 GO 术语开始（你已经在前一部分学会了如何检索 Ensembl ID）。记住，你需要使用前一部分中的do_request函数：

    lct_id = 'ENSG00000115850'
refs = do_request(ensembl_server, 'xrefs/id', lct_id,external_db='GO', all_levels='1')
print(len(refs))
print(refs[0].keys())
for ref in refs:
    go_id = ref['primary_id']
    details = do_request(ensembl_server, 'ontology/id', go_id)
    print('%s %s %s' % (go_id, details['namespace'], ref['description']))
    print('%s\n' % details['definition'])

注意自由格式的定义和每个术语的不同命名空间。在循环中报告的前两个项目如下（当你运行时，它可能会有所变化，因为数据库可能已经更新）：

GO:0000016 molecular_function lactase activity
 "Catalysis of the reaction: lactose + H2O = D-glucose + D-galactose." [EC:3.2.1.108]

 GO:0004553 molecular_function hydrolase activity, hydrolyzing O-glycosyl compounds
 "Catalysis of the hydrolysis of any O-glycosyl bond." [GOC:mah]

  1. 让我们集中注意力在乳糖酶活性分子功能上，并获取更多关于它的详细信息（以下go_id来自前一步）：

    go_id = 'GO:0000016'
my_data = do_request(ensembl_server, 'ontology/id', go_id)
for k, v in my_data.items():
    if k == 'parents':
        for parent in v:
            print(parent)
            parent_id = parent['accession']
    else:
        print('%s: %s' % (k, str(v)))
parent_data = do_request(ensembl_server, 'ontology/id', parent_id)
print(parent_id, len(parent_data['children']))

我们打印出乳糖酶活性记录（它目前是 GO 树分子功能的一个节点），并检索潜在的父节点列表。此记录只有一个父节点。我们获取该父节点并打印出它的子节点数量。

  1. 让我们检索所有与乳糖酶活性分子功能相关的一般术语（再次强调，父术语及所有其他祖先术语）：

    refs = do_request(ensembl_server, 'ontology/ancestors/chart', go_id)
for go, entry in refs.items():
    print(go)
    term = entry['term']
    print('%s %s' % (term['name'], term['definition']))
    is_a = entry.get('is_a', [])
    print('\t is a: %s\n' % ', '.join([x['accession'] for x in is_a]))

我们通过遵循is_a关系来获取祖先列表（更多关于可能关系类型的详细信息请参考另见部分中的 GO 网站）。

  1. 让我们定义一个函数，该函数将创建一个包含术语的祖先关系字典，并返回每个术语的摘要信息，作为一个对：

    def get_upper(go_id):
    parents = {}
    node_data = {}
    refs = do_request(ensembl_server, 'ontology/ancestors/chart', go_id)
    for ref, entry in refs.items():
        my_data = do_request(ensembl_server, 'ontology/id', ref)
        node_data[ref] = {'name': entry['term']['name'], 'children': my_data['children']}
        try:
            parents[ref] = [x['accession'] for x in entry['is_a']]
        except KeyError:
            pass  # Top of hierarchy
    return parents, node_data

  2. 最后，我们将打印出乳糖酶活性术语的关系树。为此，我们将使用pygraphivz库：

    parents, node_data = get_upper(go_id)
import pygraphviz as pgv
g = pgv.AGraph(directed=True)
for ofs, ofs_parents in parents.items():
    ofs_text = '%s\n(%s)' % (node_data[ofs]['name'].replace(', ', '\n'), ofs)
    for parent in ofs_parents:
        parent_text = '%s\n(%s)' % (node_data[parent]['name'].replace(', ', '\n'), parent)
        children = node_data[parent]['children']
        if len(children) < 3:
            for child in children:
                if child['accession'] in node_data:
                    continue
                g.add_edge(parent_text, child['accession'])
        else:
            g.add_edge(parent_text, '...%d...' % (len(children) - 1))
        g.add_edge(parent_text, ofs_text)
print(g)
g.graph_attr['label']='Ontology tree for Lactase activity'
g.node_attr['shape']='rectangle'
g.layout(prog='dot')
g.draw('graph.png')

以下输出显示了乳糖酶活性术语的本体树：

图 5.4 – “乳糖酶活性”术语的本体树（顶部的术语更为一般）；树的顶部是 molecular_function；对于所有祖先节点，还标注了额外后代的数量（如果少于三个，则进行列举）

图 5.4 – “乳糖酶活性”术语的本体树（顶部的术语更为一般）；树的顶部是 molecular_function；对于所有祖先节点，还标注了额外后代的数量（如果少于三个，则进行列举）

还有更多内容...

如果你对基因本体论感兴趣，主要参考网站是geneontology.org，在这里你会找到更多关于这个主题的信息。除了molecular_function，基因本体还包括生物过程和细胞组分。在我们的配方中，我们遵循了is a的分层关系，但也存在其他部分关系。例如，“线粒体核糖体”（GO:0005761）是一个细胞组分，是“线粒体基质”的一部分（参考amigo.geneontology.org/amigo/term/GO:0005761#display-lineage-tab，点击图形视图）。

与前面的配方一样，你可以下载基因本体数据库的 MySQL 转储（你可能更喜欢以这种方式与数据交互）。有关详细信息，请参阅geneontology.org/page/download-go-annotations。同样，请准备一些时间来理解关系数据库架构。此外，请注意，绘制树和图形有许多替代方案可用于 Graphviz。我们将在本书后续章节回顾这个主题。

另见

这里有一些可以进一步了解的资源：

第七章：群体遗传学

群体遗传学是研究群体中等位基因频率变化的学科，这些变化基于选择、漂变、突变和迁徙。前几章主要集中在数据处理和清理上；这是我们第一次真正推断有趣的生物学结果。

基于序列数据的群体遗传学分析有很多有趣的内容，但由于我们已经有了不少处理序列数据的操作，我们将把注意力转向其他地方。此外，我们不会在此讨论基因组结构变异，如拷贝数变异（CNVs）或倒位。我们将集中分析 SNP 数据，这是最常见的数据类型之一。我们将使用 Python 执行许多标准的群体遗传学分析，例如使用固定指数（FST）计算 F 统计量、主成分分析（PCA），并研究群体结构。

我们将主要使用 Python 作为脚本语言，来将执行必要计算的应用程序连接在一起，这是一种传统的做法。话虽如此，由于 Python 软件生态系统仍在不断发展，你至少可以使用 scikit-learn 在 Python 中执行 PCA，如我们在第十一章中将看到的那样。

在群体遗传学数据中并没有所谓的默认文件格式。这个领域的严峻现实是，存在大量的文件格式，其中大多数是为特定应用而开发的；因此，没有任何一种是通用的。尽管有一些尝试创建更通用的格式（或者只是开发一个支持多种格式的文件转换器），但这些努力的成功是有限的。更重要的是，随着我们对基因组学理解的不断深入，我们将无论如何需要新的格式（例如，支持某种之前未知的基因组结构变异）。在这里，我们将使用 PLINK（www.cog-genomics.org/plink/2.0/），该工具最初是为在人类数据上执行全基因组关联研究（GWAS）而开发的，但它有更多的应用。如果你有下一代测序（NGS）数据，可能会问，为什么不使用变异调用格式（VCF）？嗯，VCF 文件通常会进行注释，以帮助测序分析，而在这个阶段你并不需要这些（此时你应该已经拥有一个经过过滤的数据集）。如果你将单核苷酸多态性（SNP）的调用从 VCF 转换为 PLINK，你会大约节省 95%的存储空间（这是与压缩后的 VCF 相比的结果）。更重要的是，处理 VCF 文件的计算成本要远高于其他两种格式的处理成本（想象一下处理所有这些高度结构化的文本）。如果你使用 Docker，请使用镜像 tiagoantao/bioinformatics_popgen。

本章我们将覆盖以下内容：

  • 使用 PLINK 管理数据集

  • 使用 sgkit 和 xarray 进行群体遗传学分析

  • 使用 sgkit 探索数据集

  • 分析人口结构

  • 执行 PCA 分析

  • 通过混合分析调查人口结构

首先，让我们从文件格式问题的讨论开始，然后继续讨论有趣的数据分析。

使用 PLINK 管理数据集

在这里，我们将使用 PLINK 来管理我们的数据集。我们将从 HapMap 项目中的主数据集中创建适合以下食谱分析的子集。

警告

请注意，PLINK 及其他类似程序并不是为了他们的文件格式而开发的。创建人口遗传学数据的默认文件标准可能并不是一个目标。在这个领域，你需要做好格式转换的准备（为此，Python 非常适合），因为你将使用的每个应用程序可能都有自己的独特需求。从这个食谱中你要学到的最重要的一点是，使用的不是文件格式（尽管这些是相关的），而是一种“文件转换思维”。除此之外，本食谱中的一些步骤还传达了你可能希望使用的真正的分析技巧，例如，子抽样或连锁不平衡（LD-）修剪。

准备工作

在本章中，我们将使用国际 HapMap 项目的数据。你可能记得我们在第三章中使用了 1,000 基因组项目的数据，下一代测序，而 HapMap 项目在许多方面是 1,000 基因组项目的前身；它使用的是基因分型，而非全基因组测序。HapMap 项目的大多数样本都用于 1,000 基因组项目，因此如果你已经阅读了第三章中的内容，下一代测序，你就已经对该数据集（包括可用的人群）有所了解。我不会再对数据集做更多介绍，但你可以参考第三章中的内容，下一代测序，以及 HapMap 官网(www.genome.gov/10001688/international-hapmap-project)获取更多信息。请记住，我们有来自世界各地不同人群的基因分型数据。我们将按人群的缩写来引用这些人群。以下是从www.sanger.ac.uk/resources/downloads/human/hapmap3.xhtml获取的列表：

缩写人群
ASW美国西南部的非洲血统人群
CEU来自 CEPH 收藏的北欧和西欧血统的犹他州居民
CHB中国北京的汉族人
CHD科罗拉多州丹佛市的华裔居民
GIH德克萨斯州休斯顿的古吉拉特印度人
JPT日本东京的日本人
LWK肯尼亚韦布耶的卢希亚人
MXL加利福尼亚州洛杉矶的墨西哥裔人
MKK肯尼亚金亚瓦的马赛人
TSI意大利托斯卡纳地区的人
YRI尼日利亚伊巴丹的约鲁巴人

表 6.1 - 基因组计划中的群体

注意

我们将使用 HapMap 项目的数据，该项目实际上已被 1,000 基因组计划所取代。为了教学目的，教授 Python 中的群体遗传学编程技术，HapMap 项目的数据比 1,000 基因组项目更易于处理，因为数据要小得多。HapMap 样本是 1,000 基因组样本的子集。如果你从事人类群体遗传学研究，强烈建议使用 1,000 基因组计划作为基础数据集。

这将需要一个相当大的下载（大约 1GB），并且需要解压。确保你有大约 20GB 的磁盘空间用于本章。文件可以在ftp.ncbi.nlm.nih.gov/hapmap/genotypes/hapmap3_r3/plink_format/找到。

使用以下命令解压 PLINK 文件：

bunzip2 hapmap3_r3_b36_fwd.consensus.qc.poly.map.gz
bunzip2 hapmap3_r3_b36_fwd.consensus.qc.poly.ped.gz

现在，我们有了 PLINK 文件；MAP 文件包含了整个基因组中标记的位置，而 PED 文件包含了每个个体的实际标记，以及一些家谱信息。我们还下载了一个元数据文件，其中包含了每个个体的信息。查看这些文件并熟悉它们。像往常一样，这些内容也可以在Chapter06/Data_Formats.py Notebook 文件中找到，所有内容都已处理好。

最终，这个教程的大部分内容将会大量使用 PLINK（www.cog-genomics.org/plink/2.0/）。Python 主要作为连接语言来调用 PLINK。

如何操作...

请查看以下步骤：

  1. 让我们获取我们样本的元数据。我们将加载每个样本的人口信息，并记录数据集中所有其他个体的后代：

    from collections import defaultdict
f = open('relationships_w_pops_041510.txt')
pop_ind = defaultdict(list)
f.readline() # header
offspring = []
for l in f:
    toks = l.rstrip().split('\t')
    fam_id = toks[0]
    ind_id = toks[1]
    mom = toks[2]
    dad = toks[3]
    if mom != '0' or dad != '0':
        offspring.append((fam_id, ind_id))
    pop = toks[-1]
pop_ind[pop].append((fam_id, ind_id))
f.close()

这将加载一个字典，其中人口是键（CEU，YRI等），而其值是该人口中个体的列表。该字典还将存储个体是否为其他人的后代信息。每个个体通过家族和个体 ID 进行标识（这些信息可以在 PLINK 文件中找到）。HapMap 项目提供的文件是一个简单的制表符分隔文件，不难处理。虽然我们使用标准的 Python 文本处理读取文件，但这是一个典型的例子，pandas 会有所帮助。

这里有一个重要的点：之所以将此信息提供在一个单独的临时文件中，是因为 PLINK 格式没有考虑到种群结构（该格式仅为 PLINK 设计时所用的病例与对照信息提供了支持）。这并不是格式的缺陷，因为它从未被设计用来支持标准的种群遗传学研究（它是一个 GWAS 工具）。然而，这是种群遗传学中数据格式的普遍特征：无论你最终使用什么格式，都会有所缺失。

我们将在本章的其他例子中使用这些元数据。我们还将进行元数据与 PLINK 文件之间的一致性分析，但我们将把这部分推迟到下一个例子。

  1. 现在，让我们以 10% 和 1% 的标记数对数据集进行子抽样，具体如下：

    import os
os.system('plink2 --pedmap hapmap3_r3_b36_fwd.consensus.qc.poly --out hapmap10 --thin 0.1 --geno 0.1 --export ped')
os.system('plink2 --pedmap hapmap3_r3_b36_fwd.consensus.qc.poly --out hapmap1 --thin 0.01 --geno 0.1 --export ped')

在 Jupyter Notebook 中，你只需要这样做：

!plink2 --pedmap hapmap3_r3_b36_fwd.consensus.qc.poly --out hapmap10 --thin 0.1 --geno 0.1 --export ped
!plink2 --pedmap hapmap3_r3_b36_fwd.consensus.qc.poly --out hapmap1 --thin 0.01 --geno 0.1 --export ped

注意微妙之处，你实际上并不会得到 1% 或 10% 的数据；每个标记有 1% 或 10% 的机会被选中，因此你将得到大约 1% 或 10% 的标记。

显然，由于过程是随机的，不同的运行会产生不同的标记子集。这将在后续的分析中产生重要的影响。如果你想复现完全相同的结果，你仍然可以使用 --seed 选项。

我们还将去除所有基因型率低于 90% 的 SNP（使用 --geno 0.1 参数）。

注意

这段代码与 Python 本身没有什么特殊之处，但你可能有两个原因希望对数据进行子抽样。首先，如果你正在对自己的数据集进行探索性分析，你可能希望从一个较小的版本开始，因为这样更易处理。而且，你将能更全面地查看你的数据。第二，一些分析方法可能不需要全部数据（实际上，一些方法甚至可能无法使用全部数据）。不过，对于最后一点要非常小心；也就是说，对于每种分析方法，确保你理解要回答的科学问题对数据的要求。通常提供过多的数据是可以的（即便你需要支付时间和内存的代价），但提供过少的数据将导致不可靠的结果。

  1. 现在，让我们只使用常染色体生成子集（也就是说，我们将去除性染色体和线粒体），具体如下：

    def get_non_auto_SNPs(map_file, exclude_file):
    f = open(map_file)
    w = open(exclude_file, 'w')
    for l in f:
        toks = l.rstrip().split('\t')
        try:
            chrom = int(toks[0])
        except ValueError:
            rs = toks[1]
            w.write('%s\n' % rs)
    w.close()
get_non_auto_SNPs('hapmap1.map', 'exclude1.txt')
get_non_auto_SNPs('hapmap10.map', 'exclude10.txt')
os.system('plink2 –-pedmap hapmap1 --out hapmap1_auto --exclude exclude1.txt --export ped')
os.system('plink2 --pedmap hapmap10 --out hapmap10_auto --exclude exclude10.txt --export ped')

  2. 让我们创建一个函数，生成一个包含所有不属于常染色体的 SNP 列表。对于人类数据，这意味着所有非数字染色体。如果你使用的是其他物种，要小心染色体编码，因为 PLINK 是面向人类数据的。如果你的物种是二倍体，具有少于 23 条常染色体，并且有性别决定系统，即 X/Y，这将很简单；如果不是，请参考 www.cog-genomics.org/plink2/input#allow_extra_chr 以获取一些替代方案（例如 --allow-extra-chr 标志）。

  3. 然后，我们为 10%和 1%的子样本数据集创建仅包含常染色体的 PLINK 文件（分别以hapmap10_auto和hapmap1_auto为前缀）。

  4. 让我们创建一些没有后代的数据集。这些数据集将用于大多数种群遗传学分析，这些分析要求个体之间在一定程度上没有亲缘关系：

    os.system('plink2 --pedmap hapmap10_auto --filter-founders --out hapmap10_auto_noofs --export ped')

注意

这一步代表了大多数种群遗传学分析的事实，即这些分析要求样本之间有一定程度的非亲缘关系。显然，由于我们知道一些后代存在于 HapMap 中，因此我们需要将其去除。

然而，请注意，使用你的数据集时，你需要比这更加精细。例如，运行plink --genome或者使用其他程序来检测相关个体。这里的关键是，你必须花费一些精力去检测样本中的相关个体；这并不是一项微不足道的任务。

  1. 我们还将生成一个 LD 修剪后的数据集，这是许多 PCA 和混合分析算法所要求的，具体如下：

    os.system('plink2 --pedmap hapmap10_auto_noofs --indep-pairwise 50 10 0.1 --out keep --export ped')
os.system('plink2 --pedmap hapmap10_auto_noofs --extract keep.prune.in --recode --out hapmap10_auto_noofs_ld --export ped')

第一步生成一个标记列表，供数据集进行 LD 修剪时使用。这使用一个50个 SNP 的滑动窗口，每次推进10个 SNP，切割值为0.1。第二步从之前生成的列表中提取 SNP。

  1. 让我们将几个案例以不同格式重新编码：

    os.system('plink2 --file hapmap10_auto_noofs_ld --recode12 tab --out hapmap10_auto_noofs_ld_12 --export ped 12')
os.system('plink2 --make-bed --file hapmap10_auto_noofs_ld --out hapmap10_auto_noofs_ld')

第一个操作将把一个使用 ACTG 核苷酸字母的 PLINK 格式转换为另一种格式，这种格式将等位基因重新编码为1和2。我们将在执行 PCA的配方中稍后使用这个。

第二个操作将文件重新编码为二进制格式。如果你在 PLINK 中工作（使用 PLINK 提供的许多有用操作），二进制格式可能是最合适的格式（例如，提供更小的文件大小）。我们将在混合分析配方中使用此格式。

  1. 我们还将提取单一的染色体（2）进行分析。我们将从 10%子样本的常染色体数据集开始：

    os.system('plink2 --pedmap hapmap10_auto_noofs --chr 2 --out hapmap10_auto_noofs_2 --export ped')

还有更多……

你可能有许多理由想要创建不同的数据集进行分析。你可能想对数据进行一些快速的初步探索——例如，如果你计划使用的分析算法对数据格式有要求，或对输入有某些约束，例如标记数量或个体之间的关系。很可能你会有许多子集需要分析（除非你的数据集本身就非常小，例如微卫星数据集）。

这可能看起来是一个小细节，但实际上并非如此：在命名文件时要非常小心（请注意，我在生成文件名时遵循了一些简单的约定）。确保文件名能提供有关子集选项的一些信息。在进行后续分析时，你需要确保选择正确的数据集；你希望你的数据集管理是灵活且可靠的，最重要的是。最糟糕的情况是，你创建了一个包含错误数据集的分析，而这个数据集不符合软件要求的约束条件。

我们使用的 LD 剪枝对于人类分析来说是标准的，但如果你使用非人类数据，请务必检查参数。

我们下载的 HapMap 文件是基于旧版参考基因组（构建 36）。如前一章 第五章 所述，与基因组合作，如果你计划使用此文件进行更多分析，请确保使用构建 36 的注释。

本示例为接下来的示例做了准备，其结果将被广泛使用。

另请参见

使用 sgkit 进行群体遗传学分析与 xarray

Sgkit 是进行群体遗传学分析的最先进的 Python 库。它是一个现代实现，利用了几乎所有 Python 中的基础数据科学库。当我说几乎所有时，我并没有夸张；它使用了 NumPy、pandas、xarray、Zarr 和 Dask。NumPy 和 pandas 在 第二章 中已介绍。在这里，我们将介绍 xarray 作为 sgkit 的主要数据容器。因为我觉得不能要求你对数据工程库有极为深入的了解，所以我会略过 Dask 部分（主要是将 Dask 结构当作等价的 NumPy 结构来处理）。你可以在 第十一章 中找到有关超大内存 Dask 数据结构的更详细信息。

准备工作

你需要运行之前的示例，因为它的输出是本示例所需的：我们将使用其中一个 PLINK 数据集。你需要安装 sgkit。

与往常一样，这可以在 Chapter06/Sgkit.py 笔记本文件中找到，但你仍然需要运行之前的笔记本文件以生成所需的文件。

如何操作...

看一下以下步骤：

  1. 让我们加载在之前的示例中生成的 hapmap10_auto_noofs_ld 数据集：

    import numpy as np
from sgkit.io import plink
data = plink.read_plink(path='hapmap10_auto_noofs_ld', fam_sep='\t')

记住我们正在加载一组 PLINK 文件。事实证明，sgkit 为该数据创建了非常丰富且结构化的表示。这种表示基于 xarray 数据集。

  1. 让我们检查一下数据的结构 —— 如果你在笔记本中，只需输入以下内容：

    data

sgkit – 如果在笔记本中 – 将生成以下表示：

图 6.1 - sgkit 加载的 xarray 数据概览，适用于我们的 PLINK 文件

图 6.1 - sgkit 加载的 xarray 数据概览，适用于我们的 PLINK 文件

data 是一个 xarray 数据集。xarray 数据集本质上是一个字典，其中每个值都是一个 Dask 数组。就我们而言，可以假设它是一个 NumPy 数组。在这种情况下，我们可以看到我们有 56241 个变异，涵盖了 1198 个样本。每个变异有 2 个等位基因，且倍性为 2。

在笔记本中，我们可以展开每一项。在我们的案例中，我们展开了 call_genotype。这是一个三维数组，包含 variants、samples 和 ploidy 维度。该数组的类型是 int8。接下来，我们可以找到一些与条目相关的元数据，如 mixed_ploidy 和注释。最后，你将看到 Dask 实现的摘要。Array 列展示了数组的大小和形状的详细信息。对于 Chunk 列，请参阅 第十一章——但现在可以安全忽略它。

  1. 获取摘要信息的另一种方法，尤其在你没有使用笔记本时，便是检查 dims 字段：

    print(data.dims)

输出应该不言自明：

Frozen({'variants': 56241, 'alleles': 2, 'samples': 1198, 'ploidy': 2})

  1. 让我们提取一些关于样本的信息：

    print(len(data.sample_id.values))
print(data.sample_id.values)
print(data.sample_family_id.values)
print(data.sample_sex.values)

输出如下：

1198
['NA19916' 'NA19835' 'NA20282' ... 'NA18915' 'NA19250' 'NA19124']
['2431' '2424' '2469' ... 'Y029' 'Y113' 'Y076']
[1 2 2 ... 1 2 1]

我们有 1198 个样本。第一个样本的样本 ID 是 NA19916，家庭 ID 是 2431，性别为 1（男性）。请记住，考虑到 PLINK 作为数据源，样本 ID 并不足以作为主键（你可能会有多个样本具有相同的样本 ID）。主键是样本 ID 和样本家庭 ID 的复合键。

提示

你可能已经注意到我们在所有数据字段后加了 .values：这实际上是将一个懒加载的 Dask 数组渲染成一个具象的 NumPy 数组。现在，我建议你忽略它，但如果你在阅读完 第十一章 后重新回到这一章，.values 类似于 Dask 中的 compute 方法。

.values 调用并不麻烦——我们的代码之所以能工作，是因为我们的数据集足够小，可以适应内存，这对于我们的教学示例来说非常好。但如果你有一个非常大的数据集，前面的代码就过于简单了。再次强调， 第十一章 会帮助你解决这个问题。现在，简单性是为了教学目的。

  1. 在查看变异数据之前，我们必须了解 sgkit 如何存储 contigs：

    print(data.contigs)

输出如下：

['1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12', '13', '14', '15', '16', '17', '18', '19', '20', '21', '22']

这里的 contigs 是人类的常染色体（如果你的数据来自于大多数其他物种，可能就不那么幸运了——你可能会看到一些丑陋的标识符）。

  1. 现在，让我们来看看变异数据：

    print(len(data.variant_contig.values))
print(data.variant_contig.values)
print(data.variant_position.values)
print(data.variant_allele.values)
print(data.variant_id.values)

这是输出的简化版本：

56241
[ 0  0  0 ... 21 21 21]
[  557616   782343   908247 ... 49528105 49531259 49559741]
[[b'G' b'A']
 ...
 [b'C' b'A']]
['rs11510103' 'rs2905036' 'rs13303118' ... 'rs11705587' 'rs7284680'
 'rs2238837']

我们有 56241 个变异。contig 索引是 0，如果你查看前一步的步骤，它代表染色体 1。变异位于位置 557616（参照人类基因组版本 36），可能的等位基因是 G 和 A。它有一个 SNP ID，rs11510103。

  1. 最后，让我们看看 genotype 数据：

    call_genotype = data.call_genotype.values
print(call_genotype.shape)
first_individual = call_genotype[:,0,:]
first_variant = call_genotype[0,:,:]
first_variant_of_first_individual = call_genotype[0,0,:]
print(first_variant_of_first_individual)
print(data.sample_family_id.values[0], data.sample_id.values[0])
print(data.variant_allele.values[0])

call_genotype 的形状为 56,241 x 1,1198,2，这是它的变异、样本和倍性维度。

要获取第一个个体的所有变异数据，你需要固定第二个维度。要获取第一个变异的所有样本数据，你需要固定第一个维度。

如果你打印出第一个个体的详细信息（样本和家庭 ID），你会得到 2431 和 NA19916 ——如预期的那样，正如在上一个样本探索中的第一个案例。

还有更多内容...

这个食谱主要是对 xarray 的介绍，伪装成 sgkit 教程。关于 xarray 还有很多内容要说——一定要查看 docs.xarray.dev/。值得重申的是，xarray 依赖于大量的 Python 数据科学库，而我们现在暂时忽略了 Dask。

使用 sgkit 探索数据集

在这个食谱中，我们将对其中一个生成的数据集进行初步的探索性分析。现在我们对 xarray 有了一些基本了解，可以实际尝试进行一些数据分析。在这个食谱中，我们将忽略群体结构问题，这一问题将在下一个食谱中回到。

准备工作

你需要先运行第一个食谱，并确保你有 hapmap10_auto_noofs_ld 文件。此食谱中有一个 Notebook 文件，名为 Chapter06/Exploratory_Analysis.py。你还需要为上一个食谱安装的软件。

如何操作...

看一下以下步骤：

  1. 我们首先使用 sgkit 加载 PLINK 数据，方法与上一个食谱完全相同：

    import numpy as np
import xarray as xr
import sgkit as sg
from sgkit.io import plink

data = plink.read_plink(path='hapmap10_auto_noofs_ld', fam_sep='\t')

  2. 让我们请求 sgkit 提供 variant_stats：

    variant_stats = sg.variant_stats(data)
variant_stats

输出结果如下：

图 6.2 - sgkit 提供的变异统计数据

图 6.2 - sgkit 提供的变异统计数据

  1. 现在我们来看一下统计数据，variant_call_rate：

    variant_stats.variant_call_rate.to_series().describe()

这里有更多需要解释的内容，可能看起来不太明显。关键部分是 to_series() 调用。Sgkit 返回给你的是一个 Pandas 序列——记住，sgkit 与 Python 数据科学库高度集成。获得 Series 对象后，你可以调用 Pandas 的 describe 函数并得到以下结果：

count    56241.000000
mean         0.997198
std          0.003922
min          0.964107
25%          0.996661
50%          0.998331
75%          1.000000
max          1.000000
Name: variant_call_rate, dtype: float64

我们的变异呼叫率相当不错，这并不令人惊讶，因为我们在查看的是阵列数据——如果你有一个基于 NGS 的数据集，数据质量可能会更差。

  1. 现在让我们来看一下样本统计数据：

    sample_stats = sg.sample_stats(data)
sample_stats

同样，sgkit 提供了许多现成的样本统计数据：

图 6.3 - 通过调用 sample_stats 获得的样本统计数据

图 6.3 - 通过调用 sample_stats 获得的样本统计数据

  1. 现在我们来看一下样本的呼叫率：

    sample_stats.sample_call_rate.to_series().hist()

这次，我们绘制了样本呼叫率的直方图。同样，sgkit 通过利用 Pandas 自动实现这一功能：

图 6.4 - 样本呼叫率的直方图

图 6.4 - 样本呼叫率的直方图

还有更多内容...

事实上，对于人口遗传学分析，R 是无可比拟的；我们强烈建议你查看现有的 R 人口遗传学库。不要忘记，Python 和 R 之间有一个桥接，在 第一章《Python 和周边软件生态》一章中已作讨论。

如果在更大的数据集上执行，大多数在此展示的分析将会消耗大量计算资源。实际上，sgkit 已经准备好应对这个问题，因为它利用了 Dask。虽然在这个阶段介绍 Dask 过于复杂，但对于大数据集，第十一章 会讨论如何解决这些问题。

另见

  • 统计遗传学的 R 包列表可以在 cran.r-project.org/web/views/Genetics.xhtml 找到。

  • 如果你需要了解更多人口遗传学的内容，我推荐阅读 Daniel L. Hartl 和 Andrew G. Clark 合著的《人口遗传学原理》一书，出版商是 Sinauer Associates。

分析人群结构

之前，我们介绍了不考虑人群结构的 sgkit 数据分析。事实上，大多数数据集，包括我们正在使用的这个数据集，都有一定的人群结构。sgkit 提供了分析具有群体结构的基因组数据集的功能，我们将在这里进行探讨。

准备工作

你需要先运行第一个配方，并且应该已经下载了我们生成的 hapmap10_auto_noofs_ld 数据和原始的人群元数据 relationships_w_pops_041510.txt 文件。配方文件中有一个 Notebook 文件，包含了 06_PopGen/Pop_Stats.py。

如何操作...

请看以下步骤：

  1. 首先，让我们用 sgkit 加载 PLINK 数据：

    from collections import defaultdict
from pprint import pprint
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
import xarray as xr
import sgkit as sg
from sgkit.io import plink

data = plink.read_plink(path='hapmap10_auto_noofs_ld', fam_sep='\t')

  2. 现在，让我们加载数据并将个体分配到各个群体：

    f = open('relationships_w_pops_041510.txt')
pop_ind = defaultdict(list)
f.readline()  # header
for line in f:
    toks = line.rstrip().split('\t')
    fam_id = toks[0]
    ind_id = toks[1]
    pop = toks[-1]
    pop_ind[pop].append((fam_id, ind_id))
pops = list(pop_ind.keys())

我们最终得到一个字典 pop_ind，其中键是人群代码，值是样本列表。请记住，样本的主键是家庭 ID 和样本 ID。

我们还在 pops 变量中列出了人群列表。

  1. 我们现在需要告知 sgkit 每个样本属于哪个人群或队列：

    def assign_cohort(pops, pop_ind, sample_family_id, sample_id):
    cohort = []
    for fid, sid in zip(sample_family_id, sample_id):
        processed = False
        for i, pop in enumerate(pops):
            if (fid, sid) in pop_ind[pop]:
                processed = True
                cohort.append(i)
                break
        if not processed:
            raise Exception(f'Not processed {fid}, {sid}')
    return cohort
cohort = assign_cohort(pops, pop_ind, data.sample_family_id.values, data.sample_id.values)
data['sample_cohort'] = xr.DataArray(
    cohort, dims='samples')

请记住，sgkit 中的每个样本都有一个在数组中的位置。因此，我们需要创建一个数组，其中每个元素都指向样本内的特定人群或队列。assign_cohort 函数正是做了这件事：它获取我们从 relationships 文件加载的元数据和来自 sgkit 文件的样本列表，并为每个样本获取人群索引。

  1. 现在我们已经将人群信息结构加载到 sgkit 数据集中，接下来可以开始在人口或队列层级计算统计数据。首先，让我们计算每个人群的单一等位基因位点数量：

    cohort_allele_frequency = sg.cohort_allele_frequencies(data)['cohort_allele_frequency'].values
monom = {}
for i, pop in enumerate(pops):
    monom[pop] = len(list(filter(lambda x: x, np.isin(cohort_allele_frequency[:, i, 0], [0, 1]))))
pprint(monom)

我们首先请求 sgkit 计算每个队列或人群的等位基因频率。然后，我们筛选每个人群中的所有位点，其中第一个等位基因的频率为0或1（即其中一个等位基因已经固定）。最后，我们打印出来。顺便提一下，我们使用pprint.pprint函数让输出看起来更清晰（如果你希望以可读的方式渲染输出，且结构较复杂时，这个函数非常有用）：

{'ASW': 3332,
 'CEU': 8910,
 'CHB': 11130,
 'CHD': 12321,
 'GIH': 8960,
 'JPT': 13043,
 'LWK': 3979,
 'MEX': 6502,
 'MKK': 3490,
 'TSI': 8601,
 'YRI': 5172}

  1. 让我们获取每个人群所有位点的最小等位基因频率。这仍然是基于cohort_allele_frequency——因此不需要再次调用 sgkit：

    mafs = {}
for i, pop in enumerate(pops):
    min_freqs = map(
        lambda x: x if x < 0.5 else 1 - x,
        filter(
            lambda x: x not in [0, 1],
            cohort_allele_frequency[:, i, 0]))
    mafs[pop] = pd.Series(min_freqs)

我们为每个人群创建 Pandas Series对象，因为这样可以使用许多有用的功能，例如绘图。

  1. 现在，我们将打印YRI和JPT人群的 MAF 直方图。我们将利用 Pandas 和 Matplotlib 来完成此操作：

    maf_plot, maf_ax = plt.subplots(nrows=2, sharey=True)
mafs['YRI'].hist(ax=maf_ax[0], bins=50)
maf_ax[0].set_title('*YRI*')
mafs['JPT'].hist(ax=maf_ax[1], bins=50)
maf_ax[1].set_title('*JPT*')
maf_ax[1].set_xlabel('MAF')

我们让 Pandas 生成直方图，并将结果放入 Matplotlib 图中。结果如下所示：

图 6.5 - YRI 和 JPT 人群的 MAF 直方图

图 6.5 - YRI 和 JPT 人群的 MAF 直方图

  1. 现在我们将集中计算 FST。FST 是一个广泛使用的统计量，旨在表示由人群结构产生的遗传变异。让我们使用 sgkit 来计算它：

    fst = sg.Fst(data)
fst = fst.assign_coords({"cohorts_0": pops, "cohorts_1": pops})

第一行计算fst，在这种情况下，它将是队列或人群之间的配对fst。第二行通过使用 xarray 坐标功能为每个队列分配名称。这使得代码更简洁，更具声明性。

  1. 让我们比较CEU和CHB人群与CHB和CHD人群之间的fst：

    remove_nan = lambda data: filter(lambda x: not np.isnan(x), data)
ceu_chb = pd.Series(remove_nan(fst.stat_Fst.sel(cohorts_0='CEU', cohorts_1='CHB').values))
chb_chd = pd.Series(remove_nan(fst.stat_Fst.sel(cohorts_0='CHB', cohorts_1='CHD').values))
ceu_chb.describe()
chb_chd.describe()

我们将stat_FST中sel函数返回的配对结果用于比较，并创建一个 Pandas Series。请注意，我们可以按名称引用人群，因为我们已经在前一步中准备好了坐标。

  1. 让我们基于多位点配对 FST 绘制人群间的距离矩阵。在此之前，我们将准备计算：

    mean_fst = {}
for i, pop_i in enumerate(pops):
    for j, pop_j in enumerate(pops):
        if j <= i:
            continue
        pair_fst = pd.Series(remove_nan(fst.stat_Fst.sel(cohorts_0=pop_i, cohorts_1=pop_j).values))
        mean = pair_fst.mean()
        mean_fst[(pop_i, pop_j)] = mean
min_pair = min(mean_fst.values())
max_pair = max(mean_fst.values())

我们计算所有人群对之间的 FST 值。执行这段代码时将消耗大量时间和内存，因为我们实际上要求 Dask 进行大量计算，以呈现我们的 NumPy 数组。

  1. 现在我们可以对所有人群的平均 FST 进行配对绘图：

    sns.set_style("white")
num_pops = len(pops)
arr = np.ones((num_pops - 1, num_pops - 1, 3), dtype=float)
fig = plt.figure(figsize=(16, 9))
ax = fig.add_subplot(111)
for row in range(num_pops - 1):
    pop_i = pops[row]
    for col in range(row + 1, num_pops):
        pop_j = pops[col]
        val = mean_fst[(pop_i, pop_j)]
        norm_val = (val - min_pair) / (max_pair - min_pair)
        ax.text(col - 1, row, '%.3f' % val, ha='center')
        if norm_val == 0.0:
            arr[row, col - 1, 0] = 1
            arr[row, col - 1, 1] = 1
            arr[row, col - 1, 2] = 0
        elif norm_val == 1.0:
            arr[row, col - 1, 0] = 1
            arr[row, col - 1, 1] = 0
            arr[row, col - 1, 2] = 1
        else:
            arr[row, col - 1, 0] = 1 - norm_val
            arr[row, col - 1, 1] = 1
            arr[row, col - 1, 2] = 1
ax.imshow(arr, interpolation='none')
ax.set_title('Multilocus Pairwise FST')
ax.set_xticks(range(num_pops - 1))
ax.set_xticklabels(pops[1:])
ax.set_yticks(range(num_pops - 1))
ax.set_yticklabels(pops[:-1])

在下面的图表中，我们将绘制一个上三角矩阵，其中单元格的背景颜色表示分化的度量；白色表示差异较小（较低的 FST），蓝色表示差异较大（较高的 FST）。CHB和CHD之间的最小值用黄色表示，而JPT和YRI之间的最大值用洋红色表示。每个单元格中的值是这两个人群之间的平均配对 FST：

图 6.6 - HapMap 项目中 11 个人群所有常染色体的平均配对 FST

图 6.6 - HapMap 项目中 11 个人群所有常染色体的平均配对 FST

另见

  • F 统计量是一个非常复杂的话题，所以我首先将引导你到维基百科页面：en.wikipedia.org/wiki/F-statistics。

  • Holsinger 和 Weir 的论文中提供了一个很好的解释，论文标题为《Genetics in geographically structured populations: defining, estimating, and interpreting FST》，刊登在《Nature Reviews Genetics》杂志上，链接：www.nature.com/nrg/journal/v10/n9/abs/nrg2611.xhtml。

执行 PCA

PCA 是一种统计方法，用于将多个变量的维度降低到一个较小的子集，这些子集是线性无关的。它在群体遗传学中的实际应用是帮助可视化被研究个体之间的关系。

本章中的大多数食谱都使用 Python 作为粘合语言（Python 调用外部应用程序来完成大部分工作），而在 PCA 中，我们有一个选择：我们可以使用外部应用程序（例如，EIGENSOFT SmartPCA），也可以使用 scikit-learn 并在 Python 中执行所有操作。在本食谱中，我们将使用 SmartPCA——如果你想体验使用 scikit-learn 的原生机器学习方法，可以参考第十章。

提示

事实上，你还有第三个选择：使用 sgkit。然而，我想向你展示如何执行计算的替代方案。这样做有两个很好的理由。首先，你可能不想使用 sgkit——尽管我推荐它，但我不想强迫你；其次，你可能需要使用一个在 sgkit 中没有实现的替代方法。PCA 实际上就是一个很好的例子：论文的审稿人可能要求你运行一个已发布且广泛使用的方法，例如 EIGENSOFT SmartPCA。

准备工作

你需要先运行第一个食谱，才能使用hapmap10_auto_noofs_ld_12 PLINK 文件（其中等位基因已重新编码为1和2）。PCA 需要 LD 修剪后的标记；我们不会在这里使用后代数据，因为这可能会偏向结果。我们将使用重新编码后的 PLINK 文件，其中等位基因为1和2，因为这样可以使 SmartPCA 和 scikit-learn 的处理更加方便。

我有一个简单的库来帮助进行一些基因组处理。你可以在github.com/tiagoantao/pygenomics找到这个代码。你可以使用以下命令进行安装：

pip install pygenomics

对于本食谱，你需要下载 EIGENSOFT（www.hsph.harvard.edu/alkes-price/software/），其中包含我们将使用的 SmartPCA 应用程序。

在Chapter06/PCA.py食谱中有一个 Notebook 文件，但你仍然需要先运行第一个食谱。

如何做...

请查看以下步骤：

  1. 让我们加载元数据，步骤如下：

    f = open('relationships_w_pops_041510.txt')
ind_pop = {}
f.readline() # header
for l in f:
    toks = l.rstrip().split('\t')
    fam_id = toks[0]
    ind_id = toks[1]
    pop = toks[-1]
    ind_pop['/'.join([fam_id, ind_id])] = pop
f.close()
ind_pop['2469/NA20281'] = ind_pop['2805/NA20281']

在这种情况下，我们将添加一个与 PLINK 文件中已有内容一致的条目。

  1. 让我们将 PLINK 文件转换为 EIGENSOFT 格式：

    from genomics.popgen.plink.convert import to_eigen
to_eigen('hapmap10_auto_noofs_ld_12', 'hapmap10_auto_noofs_ld_12')

这使用了我编写的一个函数，将 PLINK 数据转换为 EIGENSOFT 格式。这主要是文本处理——并不是最激动人心的代码。

  1. 现在，我们将运行 SmartPCA 并解析其结果，如下所示：

    from genomics.popgen.pca import smart
ctrl = smart.SmartPCAController('hapmap10_auto_noofs_ld_12')
ctrl.run()
wei, wei_perc, ind_comp = smart.parse_evec('hapmap10_auto_noofs_ld_12.evec', 'hapmap10_auto_noofs_ld_12.eval')

同样，这将使用 pygenomics 中的几个函数来控制 SmartPCA，然后解析输出。代码是此类操作的典型代码，尽管你可以查看它，但实际上它非常直接。

parse 函数将返回 PCA 权重（我们不会使用这些权重，但你应该检查一下）、归一化权重，以及每个个体的主成分（通常最多到 PC 10）。

  1. 然后，我们绘制 PC 1 和 PC 2，如下所示的代码：

    from genomics.popgen.pca import plot
plot.render_pca(ind_comp, 1, 2, cluster=ind_pop)

这将生成以下图示。我们将提供绘图函数和从元数据中检索的种群信息，这样可以用不同的颜色绘制每个种群。结果与已发布的结果非常相似；我们将找到四个群体。大多数亚洲种群位于顶部，非洲种群位于右侧，欧洲种群位于底部。还有两个混合种群（GIH 和 MEX）位于中间：

图 6.7 - 由 SmartPCA 生成的 HapMap 数据的 PC 1 和 PC 2

图 6.7 - 由 SmartPCA 生成的 HapMap 数据的 PC 1 和 PC 2

注意

请注意，PCA 图可以在任何轴上对称，因为信号并不重要。重要的是簇应该相同，并且个体之间（以及这些簇之间）的距离应该相似。

还有更多...

这里有一个有趣的问题，你应该使用哪种方法——SmartPCA 还是 scikit-learn？我们将在第十章中使用 scikit-learn。结果是相似的，所以如果你正在进行自己的分析，可以自由选择。然而，如果你将结果发表在科学期刊上，SmartPCA 可能是更安全的选择，因为它基于遗传学领域已发布的软件；审稿人可能更倾向于选择这一方法。

另见

使用混合分析研究种群结构

在种群遗传学中，一个典型的分析是由结构程序（web.stanford.edu/group/pritchardlab/structure.xhtml）推广的，这个程序用于研究种群结构。此类软件用于推断存在多少个种群（或有多少个祖先种群生成了当前的种群），并识别潜在的迁徙者和混合个体。结构程序是在很久以前开发的，当时基因型标记数量较少（当时主要是少量微卫星标记），之后开发出了更快的版本，包括同一实验室开发的fastStructure（rajanil.github.io/fastStructure/）。在这里，我们将使用 Python 与在 UCLA 开发的同类程序——admixture（dalexander.github.io/admixture/download.xhtml）进行接口。

准备工作

你需要运行第一个步骤才能使用hapmap10_auto_noofs_ld二进制 PLINK 文件。同样，我们将使用经过 LD 修剪且没有后代的 10%自动体样本。

如同前面的步骤，你将使用pygenomics库来协助；你可以在github.com/tiagoantao/pygenomics找到这些代码文件。你可以使用以下命令安装它：

pip install pygenomics

理论上，对于本步骤，你需要下载 admixture（www.genetics.ucla.edu/software/admixture/）。然而，在本案例中，我将提供我们将使用的 HapMap 数据上运行 admixture 的输出，因为运行 admixture 需要很长时间。你可以使用已提供的结果，或者自己运行 admixture。此过程的 Notebook 文件位于Chapter06/Admixture.py中，但你仍然需要先运行该步骤。

如何做到这一点...

看一下以下步骤：

  1. 首先，让我们定义我们感兴趣的k（祖先种群的数量）范围，如下所示：

    k_range = range(2, 10)  # 2..9

  2. 让我们对所有的k进行 admixture 分析（或者，你也可以跳过此步骤，使用提供的示例数据）：

    for k in k_range:
    os.system('admixture --cv=10 hapmap10_auto_noofs_ld.bed %d > admix.%d' % (k, k))

注意

这是进行混合分析的最糟糕方式，如果你按这种方式进行，可能需要超过 3 小时。这是因为它会按顺序运行所有的k值，从2到9。有两种方法可以加速这一过程：使用混合工具提供的多线程选项（-j），或者并行运行多个应用程序。在这里，我假设最坏的情况，你只有一个核心和线程可用，但你应该能通过并行化更高效地运行。我们将在第十一章中详细讨论这个问题。

  1. 我们将需要 PLINK 文件中个体的顺序，因为混合工具以此顺序输出个体结果：

    f = open('hapmap10_auto_noofs_ld.fam')
ind_order = []
for l in f:
    toks = l.rstrip().replace(' ', '\t').split('\t')
    fam_id = toks[0]
    ind_id = toks[1]
    ind_order.append((fam_id, ind_id))
f.close()

  2. 交叉验证误差给出了“最佳”k的衡量标准，如下所示：

    import matplotlib.pyplot as plt
CVs = []
for k in k_range:
    f = open('admix.%d' % k)
    for l in f:
        if l.find('CV error') > -1:
            CVs.append(float(l.rstrip().split(' ')[-1]))
            break
    f.close()
fig = plt.figure(figsize=(16, 9))
ax = fig.add_subplot(111)
ax.set_title('Cross-Validation error')
ax.set_xlabel('K')
ax.plot(k_range, CVs)

以下图表绘制了K值从2到9之间的交叉验证误差，越低越好。从这个图表中应该可以看出，我们可能需要运行更多的K值（事实上，我们有 11 个种群；如果不是更多的话，至少应该运行到 11），但由于计算成本问题，我们停在了9。

关于是否存在“最佳”K，这将是一个非常技术性的讨论。现代科学文献表明，可能没有“最佳”K；这些结果值得一些解释。我认为在解释K结果之前，你应该意识到这一点：

图 6.8 - K 的误差

图 6.8 - K 的误差

  1. 我们需要人口信息的元数据：

    f = open('relationships_w_pops_041510.txt')
pop_ind = defaultdict(list)
f.readline() # header
for l in f:
   toks = l.rstrip().split('\t')
   fam_id = toks[0]
   ind_id = toks[1]
   if (fam_id, ind_id) not in ind_order:
      continue
   mom = toks[2]
   dad = toks[3]
   if mom != '0' or dad != '0':
      continue
 pop = toks[-1]
 pop_ind[pop].append((fam_id, ind_id))
f.close()

我们将忽略 PLINK 文件中没有的个体。

  1. 让我们加载个体组成部分，如下所示：

    def load_Q(fname, ind_order):
    ind_comps = {}
    f = open(fname)
    for i, l in enumerate(f):
        comps = [float(x) for x in l.rstrip().split(' ')]
        ind_comps[ind_order[i]] = comps
    f.close()
    return ind_comps
comps = {}
for k in k_range:
    comps[k] = load_Q('hapmap10_auto_noofs_ld.%d.Q' % k, ind_order)

混合工具会生成一个文件，包含每个个体的祖先组成部分（例如，查看任何生成的Q文件）；将会有与您选择研究的k值相等的组成部分数量。在这里，我们将加载我们研究的所有k的Q文件，并将其存储在一个字典中，个体 ID 作为键。

  1. 然后，我们将对个体进行聚类，如下所示：

    from genomics.popgen.admix import cluster
ordering = {}
for k in k_range:
    ordering[k] = cluster(comps[k], pop_ind)

请记住，个体是通过混合获得祖先种群的组成部分；我们希望根据它们在祖先组成部分上的相似性对它们进行排序（而不是按 PLINK 文件中的顺序）。这并不是一项简单的任务，需要使用聚类算法。

此外，我们并不希望对所有个体进行排序；我们希望按每个种群排序，然后再对每个种群进行排序。

为此，我在github.com/tiagoantao/pygenomics/blob/master/genomics/popgen/admix/__init__.py提供了一些聚类代码。这并不完美，但允许你执行一些看起来合理的绘图。我的代码使用了 SciPy 的聚类代码。我建议你看看（顺便说一句，改进它并不难）。

  1. 在合理的个体顺序下，我们现在可以绘制混合图：

    from genomics.popgen.admix import plot
plot.single(comps[4], ordering[4])
fig = plt.figure(figsize=(16, 9))
plot.stacked(comps, ordering[7], fig)

这将产生两个图表；第二个图表显示在下图中（第一个图表实际上是从顶部的第三个混合图的变体）。

第一个图 K = 4 需要每个个体的组分及其顺序。它将按种群排序并分割所有个体。

第二个图将绘制从 K = 2 到 9 的一组堆叠混合图。它需要一个 figure 对象（由于堆叠混合物的数量可能很多，图的尺寸会有所不同）。个体的顺序通常会遵循其中一个 K（这里我们选择了 K = 7）。

请注意，所有的 K 都值得一些解释（例如，K = 2 将非洲人口与其他人群分离开来，而 K = 3 则将欧洲人口分离，并展示了 GIH 和 MEX 的混合）：

图 6.9 - HapMap 示例的堆叠混合图（K 值介于 2 和 9 之间）

图 6.9 - HapMap 示例的堆叠混合图（K 值介于 2 和 9 之间）

还有更多内容...

不幸的是，你不能仅运行单个混合实例来获得结果。最佳做法实际上是运行 100 个实例，并获得具有最佳对数似然的结果（这在混合输出中报告）。显然，我无法要求你为此配方的每个 7 个不同的 K 运行 100 个实例（我们在谈论两周的计算时间），但如果你想要有可发布的结果，你可能必须执行这些步骤。需要集群（或至少非常好的机器）来运行这个。你可以使用 Python 处理输出并选择最佳对数似然的结果。选择每个 K 的最佳结果后，你可以轻松应用此配方来绘制输出。

avatar
文章
阅读
粉丝