Python-生物信息学秘籍-二-

Python 生物信息学秘籍（二）

原文：Bioinformatics with Python Cookbook

协议：CC BY-NC-SA 4.0

四、高级 NGS 数据处理

如果你处理下一代测序 ( NGS )数据，你就会知道质量分析和处理是获得结果的两大耗时环节。在本章的第一部分，我们将通过使用包含亲属信息的数据集来更深入地研究 NGS 分析——在我们的例子中，是母亲、父亲和大约 20 个后代。这是执行质量分析的常用技术，因为谱系信息将允许我们对我们的过滤规则可能产生的错误数量进行推断。我们还将利用这个机会，使用相同的数据集，根据现有的注释来寻找基因组特征。

本章的最后一个方法将利用 NGS 数据深入研究另一个高级主题:宏基因组学。我们将使用宏基因组学的 Python 包 QIIME2 来分析数据。

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

在本章中，有以下配方:

• 准备用于分析的数据集
• 利用孟德尔误差信息进行质量控制
• 使用标准统计数据探索数据
• 从测序注释中寻找基因组特征
• 用 QIIME2 做宏基因组学

准备用于分析的数据集

我们的起点将是一个 VCF 文件(或等效文件),其中包含由基因分型器(基因组分析工具包(在我们的例子中为 GATK )发出的调用，包括注释。由于我们将过滤 NGS 数据，我们需要可靠的决策标准来调用一个站点。那么，我们如何获得这些信息呢？一般情况下，我们不能，但如果我们需要这样做，有三种基本方法:

• 使用更强大的测序技术进行比较——例如，使用桑格测序来验证 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 而父亲是杂合的某些基因座，杂合的后代数量

在这里，我们做了 10 万次模拟。在我的例子中(这是随机的，所以你的结果可能会有所不同)，我得到了所有后代都是杂合的零个模拟。事实上，这些都是重复的排列，所以所有都是杂合的概率是或 9.5367431640625 e-07——不太可能。所以，即使对于单个后代，我们可以有 AT 或 AA；对于 20 来说，都是同一类型的可能性很小。这是我们可以用来对孟德尔错误进行不那么幼稚的解释的信息。

1. 让我们重复母亲和父亲都在的分析:

   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 或杂合的后代数量

在这种情况下，我们也有重复排列，但我们有四个可能的值，而不是两个:AA，AT，TA 和 TT。我们最终得到所有个体的相同概率为:9.53631640625 e-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 调用。我们在最后打印条目的数量(由于 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. 然后我们加载数据。我们将用熊猫来导航:

   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. 让我们请熊猫展示所有注释的直方图:

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

生成以下直方图:

图 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 感兴趣区域的直方图

我们实际上是将 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 -搜索空间某一区域的错误注释散布矩阵

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 -搜索空间的一个区域的良好调用的注释的散布矩阵

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 万个标记。

还有更多……

当丢失 12%的标记时，误差从 5%减少到 0.3%是好还是坏？嗯，那要看你接下来想做什么分析了。也许你的方法对丢失标记有弹性，但不会有太多错误，在这种情况下，这可能会有帮助。但是如果反过来，也许你更喜欢完整的数据集，即使它有更多的错误。如果你应用不同的方法，也许你会在不同的方法中使用不同的数据集。在这个按蚊数据集的具体例子中，数据太多了，减少数据量可能对任何事情都没问题。但是如果你有更少的标记，你将不得不根据标记和质量来评估你的需求。

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

我们将用一个简单的方法来结束本章和本书，这个方法表明，有时你可以从简单的意外结果中学到重要的东西，表面的质量问题可能掩盖重要的生物学问题。

我们将为我们杂交的所有亲本绘制穿过 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 kbp 窗口的中值 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

你会注意到，对于一些蚊子来说，例如，在第一列和最后一列的蚊子，在染色体臂的中间有一个明显的 DP 下降。其中一些，比如第三列，有一点下降——不那么明显。而对于第二列的底层父代，根本没有下降。

还有更多……

前面的模式有一个生物学原因，最终对测序产生影响:按蚊可能在 2L 臂中部有一个大的染色体倒位。由于进化上的差异，与用于判断的参考基因组不同的核型更难判断。这使得该区域的序列器读取次数减少。这是这个物种特有的，但你可能会认为其他种类的特征也会出现在其他生物身上。

一个更广为人知的案例是拷贝数变异 ( CNV ):如果一个参考基因组只有一个特征的拷贝，但你正在测序的个体有n，那么你可以预期看到一个n倍于整个基因组中值的 DP。

但是，在一般情况下，在整个分析过程中留意奇怪的结果是一个好主意。有时，这是一个有趣的生物特征的标志，就像这里的。或者，这是一个指向错误的指针:例如，主成分分析 ( PCA )可以用来找到错误标记的样本(因为它们可能会聚集在错误的组中)。

使用 QIIME 2 Python API 进行宏基因组学研究

维基百科称宏基因组学是对直接从环境样本中回收的遗传物质的研究。请注意，这里的“环境”应该广义地解释:在我们的例子中，我们将在患有胃肠道问题的儿童的粪便微生物群移植研究中处理胃肠道微生物群。这项研究是 QIIME 2 的教程之一，QIIME 2 是宏基因组学中最广泛使用的数据分析应用程序之一。QIIME 2 有几个接口:一个 GUI、一个命令行和一个称为工件 API 的 Python API。

Tomasz kocióek 有一个使用工件 API 的优秀教程，它基于 QIIME 2 上最完善的(基于客户端的，而不是基于工件的)教程，“移动图片”教程(nb viewer . jupyter . org/gist/tkosciol/29de 5198 a4be 81559 a 075756 c 2490 FDE)。在这里，我们将创建一个 Python 版本的粪便微生物群移植研究，可以在 docs.qiime2.org/2022.2/tuto…](docs.qiime2.org/2022.2/tuto… Tomasz 更复杂的路线:这将允许您对 QIIME 2 Python 内部有更多的了解。在你获得这种体验之后，你很可能会想走托马斯的路线，而不是我的。但是，您在这里获得的经验会让您对 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

您可能希望使用conda install在 QIIME 2 的环境中安装其他软件包。

要为 Jupyter 执行准备,您应该安装 QIIME 2 扩展，如下所示:

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

小费

该扩展具有高度的交互性，允许您从不同的角度查看本书中无法捕捉到的数据。缺点是它在nbviewer中不起作用(一些单元格输出在静态查看器中不可见)。记住与扩展的输出交互，因为许多输出是动态的。

你现在可以开始 Jupyter 了。笔记本可以在Chapter4/QIIME2_Metagenomics.py文件中找到。

警告

由于 QIIME 软件包安装的流动性，我们没有为它提供一个 Docker 环境。这意味着如果你从我们的 Docker 安装开始工作，你将不得不下载配方并手动安装软件包。

您可以找到获取笔记本文件和 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
   

下面是 Juypter 的 QIIME 扩展的一部分输出:

图 4.8-Jupyter 的 QIIME2 扩展的部分输出

记住扩展是迭代的，提供的信息远不止这张图表。

小费

该配方的原始数据以 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. 让我们使用 dada 2(【github.com/benjjneb/da…:

   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 CLI 版本中在线找到结果。

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 文件的形式提供。

准备就绪

生物基因组的大小差异很大，从 9.7 kbp 的艾滋病毒等病毒，到 9.7 kbp 的大肠杆菌等细菌，到分布在 14 条染色体、线粒体和顶体上的 22 Mbp 的原生动物，如分布在 22 条常染色体、X/Y 染色体和线粒体上的恶性疟原虫等原生动物，到具有三条常染色体、一条线粒体和 X/Y 性染色体的果蝇，再到分布在 22 条常染色体、X/Y 染色体和线粒体上的三对 Gbp 的人类一路上，你有不同的倍性和性染色体组织。

小费

如你所见，不同的生物有非常不同的基因组大小。这种差异可以有几个数量级。这可能会对您的编程风格产生重大影响。处理大基因组需要你在记忆方面更加保守。不幸的是，更大的基因组将受益于速度更高效的编程技术(因为你有更多的数据要分析)；这些是相互矛盾的要求。一般的规则是，对于更大的基因组，你必须更加注意效率(速度和内存)。

为了减轻这个食谱的负担，我们将使用来自恶性疟原虫的小真核基因组。这个基因组仍然具有许多较大基因组的典型特征(例如，多条染色体)。因此，这是复杂性和大小之间的一个很好的折衷。请注意，有了像恶性疟原虫那么大的基因组，将整个基因组加载到内存中就有可能进行许多操作。然而，我们选择了一种可以用于更大基因组(例如，哺乳动物)的编程风格，这样您就可以以更通用的方式使用这种方法，但也可以随意对像这样的小基因组使用更占用内存的方法。

我们会使用你在 第一章**Python 以及周边软件生态中安装的 Biopython。像往常一样，这份食谱可以在本书的 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 大小的窗口开始。这对于恶性疟原虫来说是合适的，但是对于染色体数量级与此不同的基因组来说，你需要考虑其他的值。

请注意，我们正在重新读取文件。对于如此小的基因组，对整个基因组进行内存加载是可行的(在步骤 1 )。无论如何，请随意为小基因组尝试这种编程风格——它更快！然而，我们的代码被设计用于更大的基因组。

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

chrom_GC字典是我们最感兴趣的数据结构，它将包含每个 50 kbp 窗口的GC内容的片段的列表。因此，对于大小为 640，851 bp 的染色体 1，将有 14 个条目，因为该染色体的大小为 14 个 50 kbp 的块。

请注意恶性疟原虫基因组的两个不同寻常的特征:基因组富含 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 中不是，在 Python 3 中,keys方法有一个特定的dict_keys返回类型。

我们按顺序画出染色体(因此排序)。我们需要最大染色体的大小(在恶性疟原虫中为 14)，以确保染色体的大小以正确的比例打印出来(变量biggest_chrom)。

然后，我们创建一个具有 PNG 输出的有机体的 A4 大小的表示。注意我们画的是非常小的 10 bp 的端粒。这将产生一个类似矩形的染色体。你可以让端粒变大，给它们一个圆形的代表，或者你可能有一个更好的想法，为你的物种使用正确的端粒大小。

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

然后，我们打印这些染色体:它们以端粒为界，由 50 kbp 的染色体片段组成(最后一个片段的大小与余数相同)。每个段将被涂成蓝色，红绿成分基于两个异常值之间的线性归一化。每个染色体片段要么是 50 kbp，要么可能更小，如果它是染色体的最后一个。输出如下图所示:

图 5.2-恶性疟原虫的 14 条染色体，用 GC 含量进行颜色编码(红色超过 22%，黄色低于 17%，蓝色阴影代表两个数字之间的线性梯度)

小费

在 python 成为如此流行的语言之前，Biopython 代码就已经发展了。过去，图书馆的可用性相当有限。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 确实有打印环形基因组的功能，你也可以用它来打印细菌。至于细菌和病毒，这些基因组更容易处理，因为它们非常小。

参见

您可以从以下资源中了解更多信息:

• 你可以在www.ensembl.org/info/data/ftp/index.xhtml的恩森布尔找到许多模式生物的参考基因组。
• 和往常一样，国家生物技术信息中心 ( NCBI )也在www.ncbi.nlm.nih.gov/genome/browse/提供了一份基因组的大名单。
• 有很多网站致力于一个单一的有机体(或一组相关的有机体)。除了你从其下载了恶性疟原虫基因组的 PlasmoDB(plasmodb.org/plasmo/)之外，你将在疾病载体的下一个配方中找到 vector base(www.vectorbase.org/)。果蝇的 fly base(flybase.org/)也值得一提，但是不要忘记搜索你感兴趣的生物。

处理低质量的基因组参考

不幸的是，并不是所有的参考基因组都具有恶性疟原虫 ?? 的品质。除了一些模式物种(例如，人类，或常见的果蝇黑腹果蝇和少数其他物种)，大多数参考基因组可以使用一些改进。在这个食谱中，我们将学习如何处理低质量的参考基因组。

准备就绪

为了与疟疾主题保持一致，我们将使用两种蚊子的参考基因组，这两种蚊子是疟疾的传播媒介:冈比亚按蚊(这是疟疾最重要的传播媒介，可以在撒哈拉以南非洲找到)和萎缩按蚊(欧洲的一种疟疾传播媒介，虽然这种疾病在欧洲已经被根除，但这种传播媒介仍然存在)。冈比亚按蚊基因组质量合理。大多数染色体已经被绘制出来，尽管 Y 染色体还需要一些工作。有一个相当大的未知染色体，可能由 X 和 Y 染色体以及中肠微生物群组成。这个基因组有合理数量的未被调用的位置(也就是说，你会发现 N s 而不是 ACTGs)。萎缩按蚊的基因组仍然是支架形式。不幸的是，这是你会发现许多非模式物种。

请注意，我们将增加一点赌注。按蚊的基因组比恶性疟原虫 ?? 的基因组大一个数量级(但仍然比大多数哺乳动物小一个数量级)。

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

怎么做...

请遵循以下步骤:

1. 让我们首先列出冈比亚按蚊基因组*的染色体:

   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)染色体是参考基因组的一大部分，相当于一个染色体臂。

不要用微小按蚊来做这个；否则，由于脚手架状态，您将获得一千多个条目。

1. 现在，让我们检查冈比亚按蚊基因组:

   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游程大小的分布。为了计算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 bps。请注意，X染色体具有最高比例的Ns位置。

1. 现在，让我们把注意力转向萎缩按蚊的基因组。让我们数一下脚手架的数量，以及脚手架尺寸的分布:

   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个支架(相对于冈比亚按蚊基因组上的 7 个条目),其中值为7811.5(?? 的平均值)。最大的支架是 5.8 Mbp，最小的支架是 1,004 bp。尺寸的第十百分位数是1537.1，而第九十百分位数是39644.7。

1. 最后，让我们绘制支架的分数——即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)
   

前面的代码将生成下图所示的输出，其中我们根据支架大小将图分成三部分:一部分用于小于 4,800 bp 的支架，一部分用于 4800 到 540,000 bp 的支架，一部分用于更大的支架。对于小型支架，Ns的分数非常低(总是低于 3.5%)；对于中型支架，差异较大(大小在 0%和 90%以上)，对于最大的支架，差异较小(在 0%和 25%之间):

图 5.3-支架的 N 分数与其尺寸的函数关系

还有更多...

有时，参考基因组携带额外的信息。例如，冈比亚按蚊基因组被软屏蔽。这意味着在基因组上运行了一些程序来识别低复杂性区域(这些区域通常更难以分析)。这可以通过大写来说明:ACTG 将是高复杂度的，而 actg 将是低复杂度的。

有许多支架的参考基因组不仅仅是一个不方便的争论。例如，非常小的支架(比如 2,000 bp 以下)在使用校准器时可能会有作图问题(例如布伦斯-惠勒校准器 ( BWA ))，特别是在极端情况下(大多数支架在其极端情况下都会有作图问题，但如果支架很小，这些问题在支架中的比例会大得多)。如果您正在使用这样的参考基因组进行比对，您将需要考虑在绘制小支架时忽略配对信息(假设您有配对末端读数)，或者至少测量支架大小对您的比对器性能的影响。在任何情况下，总的想法是，你应该小心，因为脚手架的大小和数量会不时地露出它们丑陋的头。

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

参见

以下是一些资源，您可以从中了解更多信息:

• 像 RepeatMasker 这样的工具可以用来寻找低复杂度的基因组区域。查看www.repeatmasker.org/了解更多信息。
• IUPAC 歧义代码在处理其他基因组时可能是有用的。查看www.bioinformatics.org/sms/iupac.xhtml了解更多信息。

遍历基因组注释

拥有一个基因组序列是有趣的，但是我们想要从中提取特征，比如基因、外显子和编码序列。这种类型的注释信息在通用特征格式 ( GFF )和通用传输格式 ( GTF )文件中可用。在这个菜谱中，我们将学习如何解析和分析 GFF 文件，同时以冈比亚按蚊基因组的注释为例。

准备就绪

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

您需要安装gffutils:

conda install -c bioconda gffutils

现在，我们准备开始了。

怎么做...

请遵循以下步骤:

1. 让我们首先基于我们的 GFF 文件:

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

   用gffutils创建一个注释数据库

gffutils库创建了一个 SQLite 数据库来高效地存储注释。这里，我们将尝试使用来创建数据库，但是如果它已经存在，我们将使用现有的数据库。这一步可能很耗时。

1. 现在，让我们列出所有可用的功能类型并进行计数:

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

这些特征将包括重叠群、基因、外显子、转录本等等。注意，我们将使用gffutils包的featuretypes函数。它将返回一个生成器，但是我们将把它转换成一个列表(在这里这样做是安全的)。

1. 让我们列出所有的序列号:

   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。然而，丑陋的事实是，你会发现处理文件非常困难。图书馆尽最大努力去适应这一点。事实上，这个库的大部分文档都与帮助你处理各种尴尬的变化有关(参考 pythonhosted.org/gffutils/ex…](pythonhosted.org/gffutils/ex…)

使用gffutils还有一个替代方法(要么是因为你的 GFF 文件很奇怪，要么是因为你不喜欢库接口或它对 SQL 后端的依赖)。手动解析文件。如果你看看它的格式，你会发现它并不复杂。如果您只是执行一次性操作，那么手动解析可能就足够了。当然，从长远来看，一次性操作往往没那么好。

另外，请注意注释的质量往往变化很大。随着质量的提高，复杂性也在增加。只需查看人工注释中的一个例子。你可以预期，随着时间的推移，随着我们对有机体知识的发展，注释的质量和复杂性将会增加。

参见

以下是一些资源，您可以从中了解更多信息:

• 在 www.sanger.ac.uk/resources/s… 的可以找到 GFF 的规格。
• 在 gmod.org/wiki/GFF3 的可以找到对 GFF 格式最好的解释，以及最常见的版本和 GTF。

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

在这份食谱中，我们将学习如何在注释文件的帮助下提取基因序列，以获得其相对于参考 FASTA 的坐标。我们将使用冈比亚按蚊的基因组，以及它的注释文件(按照前两个食谱)。首先，我们将提取电压门控钠通道 ( VGSC )基因，该基因涉及对杀虫剂的抗性。

准备就绪

如果你已经遵循了前两个食谱，你将准备好了。如果没有，下载冈比亚按蚊 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 中检索到的，vector base 是一个疾病媒介基因组学的在线数据库。对于其他特定的情况，你需要知道你的基因的 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. 让我们创建一个函数来为一个列表CDSs :

   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 号染色体上的 1，234 位可能指的是与第 38 号染色体上的 1，234 位不同的 SNP。对于人类数据，你可能会在 build 36 上找到很多芯片，在 build 37 上找到大量全基因组序列，而最近的人类组装是 build 38。以我们的按蚊为例，你会有第 3 和第 4 个版本。大多数物种都会这样。所以，要注意！

Python 中的 0 索引数组与 1 索引基因组数据库之间也存在问题。尽管如此，请注意一些基因组数据库也可能是 0 索引的。

还有两个混淆的来源:转录本和基因选择，就像在更丰富的注释数据库中一样。这里，您将有几个可供选择的抄本(如果您想查看一个丰富到令人困惑的数据库，请参考人工注释数据库)。此外，与编码序列相比，标记有exon的字段将包含更多信息。为此，您将需要 CDS 字段。

最后，还有一个链的问题，在这里你要根据反向互补来翻译。

参见

以下是一些资源，您可以从中了解更多信息:

• 你可以在 www.ensembl.org/info/data/m… 的下载 Ensembl 的 MySQL 表格。
• genome.ucsc.edu/的[可以找到 UCSC 基因组浏览器](genome.ucsc.edu/)。请务必查看 hgdownload.soe.ucsc.edu/downloads.x… 的下载区。
• 参照基因组，你可以在 www.ensembl.org/info/data/f… 的恩森布尔找到模式生物的 GTF。
• 关于 CDSs 和外显子的简单解释可以在www.biostars.org/p/65162/找到。

使用 Ensembl REST API 查找直系同源词

在这份食谱中，我们将学习如何寻找某个基因的直系同源物。这个简单的食谱不仅将介绍正字法检索，还将介绍如何在 web 上使用 REST APIs 来访问生物学数据。最后，但同样重要的是，它将介绍如何使用编程 API 访问 Ensembl 数据库。

在我们的例子中，我们将试图找到人类乳糖酶 ( LCT )基因在horse基因组上的任何直系同源物。

准备就绪

这个食谱不需要任何预先下载的数据，但是由于我们使用网络应用编程接口，将需要互联网接入。可以传输的数据量将会受到限制。

我们还将利用requests库来访问 Ensembl。request API 是一个易于使用的 web 请求包装器。当然，您可以使用标准的 Python 库，但是这些库要麻烦得多。

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

怎么做...

请遵循以下步骤:

1. 我们将首先创建一个支持函数来执行一个 web 请求:

   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'])
   

注意，这将为 REST 请求构建一个以前缀http://rest.ensembl.org/info/species开头的 URL。顺便说一下，前面的链接在你的浏览器上不起作用；它应该只通过 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 数据库负责管理人类基因名称并维护我们的 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)
   

你会发现不同种类的数据库，比如脊椎动物基因组注释 ( 织女星)项目、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基因组的直系同源物。但是，此代码允许您检查所有可用的正字法。

你会得到关于一个直向同源物的相当多的信息，比如直向同源物的分类级别(boreoutheria——胎盘哺乳动物是人类和马之间最接近的系统发育级别)，直向同源物的 Ensembl ID，dN/dS 比率(非同义突变到同义突变)，以及序列间差异的雪茄串(参考上一章， 第三章 ，下一代测序)。默认情况下，您还将获得直向同源序列的对齐，但是我已经将其移除以清除输出。

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 数据库在 第三章 、下一代测序中有所介绍。

以编程方式与 Ensembl 交互的另一个完全不同的策略是下载原始表，并将它们注入本地 MySQL 数据库。请注意，这本身就是一项艰巨的任务(您可能只想加载非常小的一部分表)。但是，如果您打算大量使用，您可能需要考虑创建数据库的一部分的本地版本。如果是这种情况，您可能需要重新考虑 UCSC 的替代方案，因为从本地数据库的角度来看，它和 Ensembl 一样好。

从 Ensembl 中检索基因本体信息

在这个菜谱中，您将通过查询 Ensembl REST API 再次学习如何使用基因本体信息。基因本体是受控的词汇，用于注释基因和基因产物。这些以概念树的形式提供(更一般的概念位于层次结构的顶端)。基因本体论有三个领域:细胞成分、分子功能和生物过程。

准备就绪

和前面的方法一样，我们不需要任何预先下载的数据，但是因为我们使用 web APIs，所以需要互联网访问。将要传输的数据量将是有限的。

和往常一样，你可以在Chapter05/Gene_Ontology.py笔记本文件中找到这个内容。我们将使用do_request函数，它是在前面的配方的步骤 1 中定义的(使用 Ensembl REST API 查找直系同源)。为了绘制围棋树，我们将使用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. 让我们专注于lactase activity的分子功能，并检索关于它的更详细的信息(下面的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']))
   

我们打印lactase activity记录(当前是 GO 树分子函数的一个节点)并检索一个潜在父节点列表。此记录有一个单亲。我们检索它并打印孩子的数量。

1. 让我们检索lactase activity分子函数的所有通用术语(同样，父代和所有其他祖先):

   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关系来检索ancestor列表(请参考中的 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. 最后，我们将打印出lactase activity术语的关系树。为此，我们将使用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')
   

下面的输出显示了lactase activity术语的本体树:

图 5.4–术语“乳糖酶活性”的本体树(顶部的术语更通用)；树的顶端是分子 _ 功能；对于所有的祖先节点，还会记录额外后代的数量(如果少于三个，则进行枚举)

还有更多...

如果你对基因本体论感兴趣，你的主要停靠站将是geneontology.org，在那里你将找到更多关于这个主题的信息。除了molecular_function，基因本体还有一个生物过程和一个细胞成分。在我们的食谱中，我们遵循了等级关系是一个，但是其他的确实部分存在。例如，“线粒体核糖体”(GO:0005761)是一种细胞成分，是“线粒体基质”的一部分(参见amigo . gene ontology . org/amigo/term/GO:0005761 # display-lineage-tab并点击图形视图)。

与前面的方法一样，您可以下载基因本体数据库的 MySQL 转储(您可能更喜欢以那种方式与数据交互)。关于这一点，见geneontology.org/page/download-go-annotations。同样，希望分配一些时间来理解关系数据库模式。此外，请注意，Graphviz 有许多绘制树和图的替代方法。我们将在本书的后面回到这个话题。

参见

以下是一些资源，您可以从中了解更多信息:

• 如前所述，基因本体的主要资源是 geneontology.org。
• 为了可视化，我们使用了pygraphviz库，它是 Graphviz(【www.graphviz.org】??)之上的一个包装器。
• 围棋数据有非常好的用户界面，比如 AmiGO(amigo.geneontology.org)和 quick GO(www.ebi.ac.uk/QuickGO/)。
• 用 GO 进行的最常见的分析之一是基因富集分析，以检查某些 GO 术语在某个基因组中是过表达还是低表达。geneontology.org服务器使用黑豹(go.pantherdb.org/)，但也有其他替代方案(比如大卫，在david.abcc.ncifcrf.gov/)。

六、群体遗传学

群体遗传学是在选择、漂移、突变和迁移的基础上研究群体中等位基因频率的变化。前几章主要关注数据处理和清理；这是第一章，我们将实际推断有趣的生物学结果。

有许多基于序列数据的有趣的群体遗传学分析，但由于我们已经有了相当多的处理序列数据的方法，我们将把注意力转移到其他地方。此外，我们将不涵盖基因组结构变异，如拷贝数变异 ( 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/bio informatics _ pop gen。

在本章中，我们将介绍以下配方:

• 使用 PLINK 管理数据集
• 利用 sgkit 和 xarray 进行群体遗传学分析
• 使用 sgkit 探索数据集
• 分析人口结构
• 执行 PCA
• 混合研究种群结构

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

使用 PLINK 管理数据集

这里，我们将使用 PLINK 管理我们的数据集。我们将创建我们的主数据集(来自 HapMap 项目)的子集，这些子集适合于在下面的食谱中进行分析。

警告

请注意，PLINK 和任何类似的程序都不是为它们的文件格式开发的。可能没有为群体遗传学数据创建默认文件标准的目标。在这个领域中，您需要准备好从一种格式到另一种格式的转换(对于这一点，Python 是非常合适的)，因为您将使用的每个应用程序可能都有自己古怪的需求。从这份食谱中学到的最重要的一点是，使用的不是格式，尽管这些是相关的，而是一种“文件转换心态”。除此之外，这个方法中的一些步骤也传达了真正的分析技术，你可能想考虑使用，例如，子采样或连锁不平衡- ( LD- )修剪。

准备就绪

在本章中，我们将使用来自国际单体型图项目的数据。你可能还记得，我们在第三章 、下一代测序中使用了来自 1000 个基因组项目的数据，并且 HapMap 项目在许多方面是 1000 个基因组项目的前身；基因分型取代了全基因组测序。HapMap 项目的大多数样本在 1，000 个基因组项目中使用，所以如果你已经阅读了 第三章 、下一代测序中的配方，你将已经对数据集(包括可用群体)有了一个概念。我就不多介绍数据集了，不过你可以参考 第三章 ，下一代测序，还有课程的 HapMap 网站(www . genome . gov/10001688/international-hap map-project)了解更多信息。请记住，我们有全球不同人群中许多个体的基因分型数据。我们将用这些人群的首字母缩写来称呼他们。以下是摘自www . Sanger . AC . uk/resources/downloads/human/hapmap 3 . XHTML的列表:

| **缩写** | **人口** | | 科学工作者协会(Association of Scientific Workers) | 美国西南部的非洲血统 | | CEU | 犹他州居民与北欧和西欧血统从 CEPH 收集 | | CHB | 中国北京的汉族人 | | 冠心病 | 科罗拉多州丹佛市的华人 | | GIH | 德克萨斯州休斯顿的古吉拉特印第安人 | | 发动机尾喷管温度(jet pipe temperature) | 日本东京的日本人 | | LWK | 肯尼亚韦布耶的卢赫亚 | | MXL | 加利福尼亚洛杉矶的墨西哥血统 | | MKK | 肯尼亚金亚瓦的马赛人 | | tonspersquareinch 每平方英寸吨数 | 意大利的托斯卡尼 | | YRI | 尼日利亚伊巴丹的约鲁巴语 |

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

注意

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

这将需要相当大的下载量(大约 1 GB ),并且必须解压缩。确保您有大约 20 GB 的磁盘空间用于本章。这些文件可以在FTP . NCBI . NLM . NIH . gov/hap map/genetics/hap map 3 _ 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 文件；图谱文件包含基因组中标记位置的信息，而 PED 文件包含每个个体的实际标记，以及一些谱系信息。我们还下载了一个包含每个人信息的元数据文件。看看这些文件，熟悉一下。像往常一样，这也可以在Chapter06/Data_Formats.py笔记本文件中找到，在那里一切都已经处理好了。

最后，这个食谱的大部分将大量使用 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%的子样本数据集(前缀为hapmap10_auto和hapmap1_auto)创建仅自动染色体 PLINK 文件。

4. 让我们创建一些没有后代的数据集。大多数群体遗传分析都需要这些，这在一定程度上需要不相关的个体:

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

注意

这一步代表了大多数群体遗传分析需要样本在一定程度上不相关的事实。显然，当我们知道一些后代在 HapMap 中时，我们删除它们。

但是，请注意，对于您的数据集，您应该比这更精确。例如，运行plink --genome或使用另一个程序来检测相关的个人。这里的基本点是，你必须付出一些努力来检测你的样本中的相关个体；这不是一项微不足道的任务。

1. 我们还将生成一个 LD-pruned 数据集，如许多 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-pruned，则保留该列表。这使用了一个50 SNPs 的滑动窗口，每次前进10 SNPs，切割值为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-pruning 对于人类分析来说有些标准，但是一定要检查参数，特别是如果您使用非人类数据的话。

我们下载的 HapMap 文件基于参考基因组的旧版本(build 36)。如前一章所述， 第五章 ，使用基因组，如果您打算使用这个文件进行更多的分析，请务必使用 build 36 中的注释。

该配方为后续配方奠定了基础，其结果将被广泛使用。

参见

• LD 上的维基百科页面en.wikipedia.org/wiki/Linkage_disequilibrium是一个很好的起点。
• PLINKwww.cog-genomics.org/plink/2.0/的网站记录非常完整，这是许多遗传学软件所缺乏的。

利用 sgkit 对 xarray 进行群体遗传学分析

Sgkit 是做群体遗传学分析的最先进的 Python 库。这是一个现代的实现，利用了 Python 中几乎所有的基础数据科学库。当我说几乎所有的时候，我并没有夸大其词；它使用 NumPy，pandas，xarray，Zarr 和 Dask。NumPy 和熊猫是在 第二章 中介绍的。这里，我们将介绍 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 为我们的 PLINK 文件加载的 xarray 数据的概述

data是一个 xarray 数据集。xarray 数据集本质上是一个字典，其中每个值都是一个 Dask 数组。出于我们的目的，您可以假设它是一个 NumPy 数组。在这种情况下，我们可以看到我们有 1198 样本的 56241 变体。我们每个变体有 2 个等位基因，倍性为 2 。

在笔记本中，我们可以展开每个条目。在我们的例子中，我们扩展了call_genotype。这是一个三维数组，有variants、samples和ploidy个维度。数组的类型是int8。在此之后，我们可以找到一些与条目、mixed_ploidy和注释相关的元数据。最后，您对 Dask 实现进行了总结。数组列显示了数组的大小和形状的详细信息。对于块列，参见 第十一章——但是你现在可以放心地忽略它。

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方法。

这个调用并不麻烦——我们的代码工作的原因是我们的数据集足够小，可以放入内存，这对我们的教学例子来说非常好。但是如果你有一个非常大的数据集，前面的代码就太天真了。还是那句话， 第十一章 会帮你搞定。目前，这种简单是教学用的。

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 IDrs11510103。

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 教程。关于哈维还有很多要说的——一定要看看 docs.xarray.dev/。值得重申的是，xar… 依赖于大量的 Python 数据科学库，我们暂时忽略了 Dask。

使用 sgkit 探索数据集

在这个菜谱中，我们将对我们生成的数据集之一进行初步探索性分析。现在我们对 xarray 有了一些基本的了解，我们可以实际尝试做一些数据分析。在这个配方中，我们将忽略人口结构，这个问题我们将在下一个配方中再讨论。

准备就绪

您需要运行第一个配方，并且应该有可用的hapmap10_auto_noofs_ld文件。有一个笔记本文件里有这个食谱，叫做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-SG kit 的 variant_stats 提供的变量统计数据

1. 现在让我们看看统计数据，variant_call_rate :

   variant_stats.variant_call_rate.to_series().describe()
   

这里要解开的比看起来的要多。最基本的部分是to_series()调用。Sgkit 将向您返回一个熊猫系列——请记住，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 获得的样本统计数据

1. 现在我们将让来看一下样本看涨期权价格:

   sample_stats.sample_call_rate.to_series().hist()
   

这一次，我们绘制了样本调用率的直方图。同样，sgkit 通过利用 Pandas 免费获得此功能:

图 6.4 -样本呼叫率直方图

还有更多...

事实是，对于群体基因分析来说，没有什么能打败 R；绝对鼓励你去看看现有的群体遗传学 R 库。别忘了还有一个 Python-R 桥，在 第一章 、 Python 以及周边的软件生态中讨论过。

如果在更大的数据集上进行，这里介绍的大部分分析在计算上将是昂贵的。事实上，sgkit 已经准备好处理这个问题，因为它利用了 Dask。在这个阶段引入 Dask 会太复杂，但是对于大型数据集， 第十一章 将讨论解决这些问题的方法。

参见

• 统计遗传学的 R 包列表可在cran.r-project.org/web/views/Genetics.xhtml获得。
• 如果你需要了解更多关于群体遗传学的知识，我推荐《群体遗传学原理》这本书，作者是的丹尼尔·l·哈特尔和的安德鲁·g·克拉克。**

分析人口结构

之前，我们用 sgkit 引入了数据分析，忽略了群体结构。大多数数据集，包括我们正在使用的数据集，实际上都有一个群体结构。Sgkit 提供了用群体结构分析基因组数据集的功能，这就是我们在这里要研究的内容。

准备就绪

您需要运行第一个配方，并且应该下载我们生成的hapmap10_auto_noofs_ld数据和原始群体元数据relationships_w_pops_041510.txt文件。有一个笔记本文件，里面有06_PopGen/Pop_Stats.py的配方。

怎么做...

看看下面的步骤:

1. 首先，让我们用 sgkit:

   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')
   

   加载 PLINK 数据

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')
   

记住 SG kit 中的每个样本在一个数组中都有一个位置。因此，我们必须创建一个数组，其中的每个元素都指向样本中的特定人群或群体。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】

我们为每个种群创建 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 直方图

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 :

   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()
   

   人口之间的fst

我们从stat_FST中获取由sel函数返回的成对结果，用它来比较和创建熊猫系列。请注意，我们可以通过名称来引用群体，因为我们在上一步中已经准备好了坐标。

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-hap map 项目中所有常染色体的 11 个群体的平均成对 FST

参见

• F-statistics 是一个非常复杂的话题，所以我将首先带你去位于 en.wikipedia.org/wiki/F-stat… ?? 的维基百科页面。
• 在 www.nature.com/nrg/journal… FST* )中可以找到一个很好的解释。](www.nature.com/nrg/journal…)

执行 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-修剪标记；我们不会冒险在这里使用后代，因为它可能会使结果有偏差。我们将使用带有等位基因1和2的重新编码的 PLINK 文件，因为这使得用 SmartPCA 和 scikit-learn 进行处理更容易。

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

pip install pygenomics

对于这个食谱，你将需要下载 eigen soft(www.hsph.harvard.edu/alkes-price/software/)，其中包括我们将使用的 SmartPCA 应用程序。

在Chapter06/PCA.py配方中有一个笔记本文件，但是您仍然需要运行第一个配方。

怎么做...

看看下面的步骤:

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 和墨西哥)位于中间:

图 6.7-smart PCA 产生的 HapMap 数据的 PC 1 和 PC 2

注意

注意，PCA 图可以在运行的任何轴上对称，因为信号无关紧要。重要的是聚类应该是相同的，并且个体(和这些聚类)之间的距离应该是相似的。

还有更多...

这里一个有趣的问题是你应该使用哪种方法——smart PCA 还是 scikit-learn，我们会在 第十章 中用到。结果是相似的，所以如果您正在执行自己的分析，您可以自由选择。然而，如果你在科学杂志上发表你的结果，SmartPCA 可能是一个更安全的选择，因为它是基于遗传学领域的软件发表的；评论家可能会喜欢这个。

参见

• 可能普及了 PCA 在遗传学中的应用的论文是 11 月 bre 等人的Genes mirror geography within EuropeonNature，其中欧洲人的 PCA 几乎完美地映射到欧洲地图上。这个可以在www . nature . com/nature/journal/v 456/n 7218/ABS/nature 07331 . XHTML找到。请注意，没有任何关于 PCA 的内容可以保证它会映射到地理特征(只需提前检查我们的 PCA)。
• journals.plos.org/plosgenetic… Patterson 等人的种群结构和特征分析、 PLoS 遗传学中描述了 SmartPCAid = 10.1371/journal . pgen . 0020190](journals.plos.org/plosgenetic…
• 关于 PCA 含义的讨论可以在麦克维恩在 journals.plos.org/plosgenetic… PLoS Genetics 的论文中找到 id = 10.1371/journal . pgen . 1000686](journals.plos.org/plosgenetic…

用混合物研究种群结构

人口遗传学中一个典型的分析是由程序结构(web.stanford.edu/group/pritchardlab/structure.xhtml)推广的，用于研究人口结构。这种类型的软件用于推断存在多少种群(或多少祖先种群产生了当前种群)，并识别潜在的移民和混合个体。这个结构是在相当一段时间以前开发的，当时进行基因分型的标记要少得多(当时，这主要是少数微卫星)，并且开发了更快的版本，包括来自同一实验室的一个叫做fastStructure(rajanil.github.io/fastStructure/)。在这里，我们将使用 Python 与加州大学洛杉矶分校开发的相同类型的程序接口，该程序被称为混合程序(dalexander.github.io/admixture/download.xhtml)。

准备就绪

为了使用hapmap10_auto_noofs_ld二进制 PLINK 文件，您需要运行第一个配方。同样，我们将使用 10%的被 LD 修剪过的没有后代的常染色体的二次抽样。

正如在前面的食谱中，你将使用的pygenomics库来帮助；你可以在 github.com/tiagoantao/… 找到这些代码文件。您可以使用以下命令安装它:

pip install pygenomics

理论上，对于这个食谱，你需要下载混合物(www.genetics.ucla.edu/software/admixture/)。然而，在这种情况下，我将在我们将使用的 HapMap 数据上提供运行混合物的输出，因为运行混合物需要很多时间。您可以使用可用的结果或自己运行混合物。在Chapter06/Admixture.py配方中有一个笔记本文件，但是你仍然需要先运行配方。

怎么做...

看看下面的步骤:

1. 首先，让我们定义我们感兴趣的k(一些祖先种群)范围，如下:

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

2. 让我们为所有的k运行混合(或者，您可以跳过这一步，使用提供的示例数据):

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

注意

这是运行混合物最糟糕的方式，如果你这样做的话，可能要花 3 个多小时。这是因为它将按顺序运行从2到9的所有k。有两件事可以加快速度:使用混合提供的多线程选项(-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)
   

下面的图绘制了2的 a K和9之间的CV，越低越好。从这张图中可以清楚地看出，我们也许应该运行更多的K(实际上，我们有 11 个种群；如果不是更多，我们应该至少运行到 11)，但由于计算成本，我们停在9。

关于是否存在“最佳”这种东西，这将是一场非常技术性的辩论。现代科学文献表明，可能不存在一个“最佳”K；这些结果值得做一些解释。我认为在你继续解释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/py genomics/blob/master/genomics/pop gen/mixed/_ _ 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(我们在这里选择了7中的一个K)。

请注意，所有的K都值得进行一些解释(例如，K = 2将非洲人口与其他人分开，K = 3将欧洲人口分开，并显示了 GIH 和墨西哥的混合):

图 6.9-hap map 示例的堆积混合图(K 在 2 和 9 之间)

还有更多...

不幸的是，你不能运行一个混合物的实例来得到一个结果。最佳实践是实际运行 100 个实例，并获得具有最佳对数似然性的实例(在混合输出中报告)。显然，我不能要求你为这个食谱的 7 个不同的K中的每一个运行 100 个实例(我们正在谈论两周的计算)，但是如果你想要有可发布的结果，你可能需要来执行这个。运行这个需要一个集群(或者至少是一个非常好的机器)。您可以使用 Python 检查输出并选择最佳对数似然。在为每个K选择了具有最佳对数似然的结果后，您可以很容易地应用这个方法来绘制输出。