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

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

译者：飞龙

协议：CC BY-NC-SA 4.0

第十一章：生物信息学中的机器学习

机器学习在许多不同的领域中都有应用，计算生物学也不例外。机器学习在该领域有着无数的应用，最古老且最为人熟知的应用之一就是使用主成分分析（PCA）通过基因组学研究种群结构。随着该领域的蓬勃发展，还有许多其他潜在的应用。在本章中，我们将从生物信息学的角度介绍机器学习的概念。

鉴于机器学习是一个非常复杂的主题，足以填满一本书，在这里我们打算采用直观的方法，让你大致了解一些机器学习技术如何帮助解决生物学问题。如果你发现这些技术有用，你将理解基本概念，并可以进一步阅读更详细的文献。

如果你使用的是 Docker，由于本章中的所有库对于数据分析都是基础库，它们都可以在 Docker 镜像tiagoantao/bioinformatics_ml中找到。

在本章中，我们将覆盖以下配方：

• 用 PCA 示例介绍 scikit-learn

• 使用 PCA 聚类来分类样本

• 使用决策树探索乳腺癌特征

• 使用随机森林预测乳腺癌结果

用 PCA 示例介绍 scikit-learn

PCA 是一种统计方法，用于将多个变量的维度降低到一个较小的、线性无关的子集。在第六章中，我们已经看过基于使用外部应用程序的 PCA 实现。在这个配方中，我们将为种群遗传学实现相同的 PCA，但将使用scikit-learn库。Scikit-learn 是 Python 中用于机器学习的基础库之一，本配方是对该库的介绍。PCA 是一种无监督的机器学习方法——我们不提供样本的类别信息。我们将在本章的其他配方中讨论有监督的技术。

提醒一下，我们将计算来自 HapMap 项目的 11 个人类种群的 PCA。

准备工作

你需要运行第六章中的第一个配方，以生成hapmap10_auto_noofs_ld_12 PLINK 文件（等位基因记录为 1 和 2）。从种群遗传学的角度来看，我们需要 LD 修剪标记来生成可靠的 PCA。我们这里不使用后代样本，因为它可能会导致结果偏差。我们的配方将需要pygenomics库，可以使用以下命令进行安装：

pip install pygenomics

代码位于Chapter10/PCA.py笔记本中。

如何实现…

看一下以下步骤：

1. 我们首先加载样本的元数据。在我们的案例中，我们将加载每个样本所属的人类种群：

   import os
   from sklearn.decomposition import PCA
   import numpy as np
   from genomics.popgen.pca import plot
   f = open('../Chapter06/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()
   

2. 我们现在获得了个体的顺序以及我们将要处理的 SNP 数量：

   f = open('../Chapter06/hapmap10_auto_noofs_ld_12.ped')
   ninds = 0
   ind_order = []
   for line in f:
       ninds += 1
       toks = line[:100].replace(' ', '\t').split('\t')
       fam_id = toks[0]
       ind_id = toks[1]
       ind_order.append('%s/%s' % (fam_id, ind_id))
   nsnps = (len(line.replace(' ', '\t').split('\t')) - 6) // 2
   f.close()
   

3. 我们创建将要输入 PCA 的数组：

   pca_array = np.empty((ninds, nsnps), dtype=int)
   print(pca_array.shape)
   f = open('../Chapter06/hapmap10_auto_noofs_ld_12.ped')
   for ind, line in enumerate(f):
       snps = line.replace(' ', '\t').split('\t')[6:]
       for pos in range(len(snps) // 2):
           a1 = int(snps[2 * pos])
           a2 = int(snps[2 * pos])
           my_code = a1 + a2 - 2
           pca_array[ind, pos] = my_code
   f.close()
   

4. 最后，我们计算包含最多八个成分的 PCA。然后，使用transform方法获取所有样本的 8 维坐标。

   my_pca = PCA(n_components=8)
   my_pca.fit(pca_array)
   trans = my_pca.transform(pca_array)
   

5. 最后，我们绘制 PCA 图：

   sc_ind_comp = {}
   for i, ind_pca in enumerate(trans):
       sc_ind_comp[ind_order[i]] = ind_pca
   plot.render_pca_eight(sc_ind_comp, cluster=ind_pop)
   

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_10_01.jpg

图 10.1 - 由 scikit-learn 生成的我们数据集的 PC1 到 PC8

还有更多…

对于科学期刊中的发布，我建议使用第六章中的食谱，因为它基于一个已发布的、广受好评的方法。也就是说，这段代码的结果在定性上是相似的，并且以非常类似的方式对数据进行了聚类（如果你与第六章中的图进行比较，垂直轴方向的反转对于解读 PCA 图表而言是无关紧要的）。

使用 PCA 进行聚类来对样本进行分类

基因组学中的 PCA 可以让我们看到样本如何聚类。在许多情况下，同一群体的个体会聚集在图表的同一区域。但我们希望进一步预测新个体在群体中的位置。为此，我们将从 PCA 数据开始，因为它进行了降维处理—使得处理数据更为简便—然后应用 K-Means 聚类算法来预测新样本的位置。我们将使用与上述食谱相同的数据集。我们将使用除了一个样本外的所有样本来训练算法，然后预测剩下的样本位置。

K-Means 聚类可以是监督算法的一个例子。在这类算法中，我们需要一个训练数据集，以便算法能够学习。训练算法之后，它将能够对新样本预测某个结果。在我们的案例中，我们希望能够预测群体。

警告

当前的食谱旨在作为对监督算法及其背后概念的温和介绍。我们训练算法的方式远非最佳。关于如何正确训练一个监督算法的问题将在本章最后一个食谱中提到。

准备工作

我们将使用与之前食谱中相同的数据。此食谱的代码可以在Chapter10/Clustering.py中找到。

如何操作…

让我们来看一下：

1. 我们首先加载群体信息——这与我们在之前的食谱中所做的相似：

   import os
   import matplotlib.pyplot as plt
   from sklearn.cluster import KMeans
   from sklearn.decomposition import PCA
   import numpy as np
   from genomics.popgen.pca import plot
   f = open('../Chapter06/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()
   
   f = open('../Chapter06/hapmap10_auto_noofs_ld_12.ped')
   ninds = 0
   ind_order = []
   for line in f:
       ninds += 1
       toks = line[:100].replace(' ', '\t').split('\t') #  for speed
       fam_id = toks[0]
       ind_id = toks[1]
       ind_order.append('%s/%s' % (fam_id, ind_id))
   nsnps = (len(line.replace(' ', '\t').split('\t')) - 6) // 2
   print (nsnps)
   f.close()
   

2. 现在，我们将所有样本数据（SNPs）加载到一个 NumPy 数组中：

   all_array = np.empty((ninds, nsnps), dtype=int)
   f = open('../Chapter06/hapmap10_auto_noofs_ld_12.ped')
   for ind, line in enumerate(f):
       snps = line.replace(' ', '\t').split('\t')[6:]
       for pos in range(len(snps) // 2):
           a1 = int(snps[2 * pos])
           a2 = int(snps[2 * pos])
           my_code = a1 + a2 - 2
           all_array[ind, pos] = my_code
   f.close()
   

3. 我们将数组分成两个数据集，即包含所有个体（除一个外）作为训练集，和单个个体用于测试的案例：

   predict_case = all_array[-1, :]
   pca_array = all_array[:-1,:]
   
   last_ind = ind_order[-1]
   last_ind, ind_pop[last_ind]
   

我们的测试案例是个体 Y076/NA19124，我们知道他属于约鲁巴族群体。

1. 现在，我们计算用于 K-Means 聚类的训练集的 PCA：

   my_pca = PCA(n_components=2)
   my_pca.fit(pca_array)
   trans = my_pca.transform(pca_array)
   
   sc_ind_comp = {}
   for i, ind_pca in enumerate(trans):
       sc_ind_comp[ind_order[i]] = ind_pca
   plot.render_pca(sc_ind_comp, cluster=ind_pop)
   

这是输出结果，将有助于检查聚类结果：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_10_02.jpg

图 10.2 - PC1 和 PC2 与按种群编码的颜色

1. 在我们开始计算 K 均值聚类之前，先写一个函数来绘制运行算法后的聚类面：

   def plot_kmeans_pca(trans, kmeans):
       x_min, x_max = trans[:, 0].min() - 1, trans[:, 0].max() + 1
       y_min, y_max = trans[:, 1].min() - 1, trans[:, 1].max() + 1
       mesh_x, mesh_y = np.meshgrid(np.arange(x_min, x_max, 0.5), np.arange(y_min, y_max, 0.5))
   
       k_surface = kmeans.predict(np.c_[mesh_x.ravel(), mesh_y.ravel()]).reshape(mesh_x.shape)
       fig, ax = plt.subplots(1,1, dpi=300)
       ax.imshow(
           k_surface, origin="lower", cmap=plt.cm.Pastel1,
           extent=(mesh_x.min(), mesh_x.max(), mesh_y.min(), mesh_y.max()),
       )
       ax.plot(trans[:, 0], trans[:, 1], "k.", markersize=2)
       ax.set_title("KMeans clustering of PCA data")
       ax.set_xlim(x_min, x_max)
       ax.set_ylim(y_min, y_max)
       ax.set_xticks(())
       ax.set_yticks(())
       return ax
   

2. 现在让我们用我们的样本来拟合算法。因为我们有 11 个人群，我们将训练 11 个簇：

   kmeans11 = KMeans(n_clusters=11).fit(trans)
   plot_kmeans_pca(trans, kmeans11)
   

这里是输出结果：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_10_03.jpg

图 10.3 - 11 个簇的簇面

如果你与这里的图进行比较，你会直观地看到这个聚类没有多大意义：它并没有很好地映射到已知的人群上。有人可能会认为，使用 11 个簇的聚类算法并不是很有用。

提示

scikit-learn 中实现了许多其他聚类算法，在多种情况下，它们的表现可能优于 K 均值聚类。你可以在scikit-learn.org/stable/modules/clustering.xhtml找到它们。值得怀疑的是，在这个特定的案例中，任何替代方法都不太可能在 11 个簇的情况下表现得更好。

1. 尽管看起来 K 均值聚类无法解决 11 个人群的划分，但如果我们使用不同数量的簇，或许它仍然可以提供一些预测。仅通过查看图表，我们可以看到四个独立的块。如果我们使用四个簇，会得到什么结果呢？

   kmeans4 = KMeans(n_clusters=4).fit(trans)
   plot_kmeans_pca(trans, kmeans4)
   

这里是输出结果：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_10_04.jpg

图 10.4 - 四个簇的簇面

四个群体现在大致清晰了。但它们直观上有意义吗？如果有，那么我们可以利用这种聚类方法。事实上，它们确实有意义。左侧的簇由非洲人口组成，顶部的簇由欧洲人组成，底部的簇由东亚人组成。中间的簇最为难以理解，因为它包含了古吉拉特人和墨西哥后裔，但这种混合最初来自于主成分分析（PCA），而非聚类本身所致。

1. 让我们看看我们漏掉的那个单独个案的预测结果：

   pca_predict = my_pca.transform([predict_case])
   kmeans4.predict(pca_predict)
   

我们的样本被预测为属于簇 1。现在我们需要进一步挖掘一下。

1. 让我们来看看簇 1 代表什么。我们取训练集中的最后一个个体，他也是一名约鲁巴人，看看他被分配到了哪个簇：

   last_train = ind_order[-2]
   last_train, ind_pop[last_train]
   kmeans4.predict(trans)[0]
   

确实是簇 1，因此预测是正确的。

还有更多…

值得重申的是，我们正在努力实现对机器学习的直观理解。在这个阶段，你应该对监督学习能够带来什么有所了解，并且已经掌握了一个聚类算法的使用示例。关于训练机器学习算法的过程，还有很多内容值得探讨，我们将在最后的食谱中部分揭示。

使用决策树探索乳腺癌特征

当我们接收到一个数据集时，第一个问题之一就是决定从哪里开始分析。一开始，往往会有一种迷茫的感觉，不知道该先做什么。这里，我们将展示基于决策树的探索性方法。决策树的最大优点是，它们会给我们提供构建决策树的规则，让我们初步了解数据的情况。

在这个示例中，我们将使用一个包含乳腺癌患者特征观察的数据集。该数据集包含 699 条数据，包含如肿块厚度、细胞大小的均匀性或染色质类型等信息。结果是良性或恶性肿瘤。特征值被编码为 0 到 10 之间的值。关于该项目的更多信息可以在archive.ics.uci.edu/ml/datasets/breast+cancer+wisconsin+%28diagnostic%29中找到。

准备工作

我们将下载数据及其文档：

wget http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.data
wget http://archive.ics.uci.edu/ml/machine-learning-databases/breast-cancer-wisconsin/breast-cancer-wisconsin.names

数据文件的格式为 CSV 文件。关于内容的信息可以在第二个下载的文件中找到。

这个代码可以在Chapter10/Decision_Tree.py中找到。

如何操作…

按照以下步骤操作：

1. 我们做的第一件事是移除一小部分数据不完整的个体：

   import numpy as np
   import matplotlib.pyplot as plt
   import pandas as pd
   from sklearn import tree
   f = open('breast-cancer-wisconsin.data')
   w = open('clean.data', 'w')
   for line in f:
       if line.find('?') > -1:
           continue
       w.write(line)
   f.close()
   w.close()
   

提示

移除数据不完整的个体在这种情况下是合适的，因为它们只是数据集的一小部分，而且我们只是进行探索性分析。如果数据缺失很多，或者我们需要做更严格的分析，你需要使用处理缺失数据的方法，但我们这里不会进行探讨。

1. 现在我们将读取数据，为所有列命名：

   column_names = [
       'sample_id', 'clump_thickness', 'uniformity_cell_size',
       'uniformity_cell shape', 'marginal_adhesion',
       'single_epithelial_cell_size', 'bare_nuclei',
       'bland_chromatin', 'normal_nucleoli', 'mitoses',
       'class'
   ]
   samples = pd.read_csv('clean.data', header=None, names=column_names, index_col=0)
   

2. 现在我们将特征与结果分离，并使用 0 和 1 重新编码结果：

   training_input = samples.iloc[:,:-1]
   target = samples.iloc[:,-1].apply(lambda x: 0 if x == 2 else 1)
   

3. 现在让我们基于这些数据创建一个最大深度为 3 的决策树：

   clf = tree.DecisionTreeClassifier(max_depth=3)
   clf.fit(training_input, target)
   

4. 让我们先看看哪些特征最重要：

   importances = pd.Series(
       clf.feature_importances_ * 100,
       index=training_input.columns).sort_values(ascending=False)
   importances
   

以下是按重要性排名的特征：

uniformity_cell_size           83.972870
uniformity_cell shape           7.592903
bare_nuclei                     4.310045
clump_thickness                 4.124183
marginal_adhesion               0.000000
single_epithelial_cell_size     0.000000
bland_chromatin                 0.000000
normal_nucleoli                 0.000000
mitoses                         0.000000

记住，这只是探索性分析。在下一个配方中，我们将尝试生成更可靠的排名。底部特征为零的原因是我们要求最大深度为 3，在这种情况下，可能并非所有特征都会被使用。

1. 我们可以对实现的准确性进行一些原生分析：

   100 * clf.score(training_input, target)
   

我们得到的性能为 96%。我们不应该用训练集测试算法，因为这会导致结果较为循环。我们将在下一个配方中重新审视这一点。

1. 最后，让我们绘制决策树：

   fig, ax = plt.subplots(1, dpi=300)
   tree.plot_tree(clf,ax=ax, feature_names=training_input.columns, class_names=['Benign', 'Malignant'])
   

这会产生以下输出：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_10_05.jpg

图 10.5 - 乳腺癌数据集的决策树

让我们从根节点开始看：它的标准是uniformity_cell_size < 2.5，分类为良性。分裂树的主要特征是细胞大小的一致性。根节点的良性分类仅仅是因为数据集中的大多数样本都是良性的。现在看一下根节点右侧的节点：它有 265 个样本，其中大部分是恶性的，并且标准为uniformity_cell_shape < 2.5。

这些规则帮助你初步理解可能驱动数据集的因素。决策树的精确度不是很高，所以将这些视为你的初步步骤。

使用随机森林预测乳腺癌结果

现在我们将使用随机森林预测一些患者的结果。随机森林是一种集成方法（它将使用多个机器学习算法的实例），通过多个决策树得出关于数据的稳健结论。我们将使用与前一个例子相同的例子：乳腺癌特征和结果。

这个例子有两个主要目标：介绍随机森林及机器学习算法训练中的问题。

准备开始

这个例子的代码可以在Chapter10/Random_Forest.py找到。

如何操作…

看一下代码：

1. 我们像前一个例子一样，首先去除缺失信息的样本：

   import pandas as pd
   import numpy as np
   import pandas as pd
   from sklearn.ensemble import RandomForestClassifier
   from sklearn.model_selection import train_test_split
   from sklearn.tree import export_graphviz
   f = open('breast-cancer-wisconsin.data')
   w = open('clean.data', 'w')
   for line in f:
       if line.find('?') > -1:
           continue
       w.write(line)
   f.close()
   w.close()
   

2. 现在我们加载清理过的数据：

   column_names = [
       'sample_id', 'clump_thickness', 'uniformity_cell_size',
       'uniformity_cell shape', 'marginal_adhesion',
       'single_epithelial_cell_size', 'bare_nuclei',
       'bland_chromatin', 'normal_nucleoli', 'mitoses',
       'class'
   ]
   samples = pd.read_csv('clean.data', header=None, names=column_names, index_col=0)
   samples
   
   

3. 我们将数据分为特征和结果：

   training_input = samples.iloc[:, :-1]
   target = samples.iloc[:, -1]
   

4. 我们创建一个分类器并将数据拟合到它：

   clf = RandomForestClassifier(max_depth=3, n_estimators=200)
   clf.fit(training_input, target)
   

这里最重要的参数是n_estimators：我们要求构建一个由 200 棵树组成的森林。

1. 现在我们按重要性对特征进行排序：

   importances = pd.Series(
       clf.feature_importances_ * 100,
       index=training_input.columns).sort_values(ascending=False)
   importances
   

以下是输出：

uniformity_cell_size           30.422515
uniformity_cell shape          21.522259
bare_nuclei                    18.410346
single_epithelial_cell_size    10.959655
bland_chromatin                 9.600714
clump_thickness                 3.619585
normal_nucleoli                 3.549669
marginal_adhesion               1.721133
mitoses                         0.194124

结果是非确定性的，意味着你可能会得到不同的结果。另外，请注意，随机森林与前一个例子中的决策树有很大的不同。这是预期之中的，因为决策树是一个单一的估计器，而随机森林权衡了 200 棵树，因此更加可靠。

1. 我们可以对这个案例进行评分：

   clf.score(training_input, target)
   

我得到的结果是 97.95%。你可能会得到稍微不同的值，因为算法是随机的。正如我们在前一个例子中所说，从训练集获取得分是一个循环过程，远非最佳实践。

1. 为了更真实地了解算法的准确性，我们需要将数据分为两部分——训练集和测试集：

   for test_size in [0.01, 0.1, 0.2, 0.5, 0.8, 0.9, 0.99]:
       X_train, X_test, y_train, y_test = train_test_split(
           trainning_input, target, test_size=test_size)
       tclf = RandomForestClassifier(max_depth=3)
       tclf.fit(X_train, y_train)
       score = tclf.score(X_test, y_test)
       print(f'{1 - test_size:.1%} {score:.2%}')
   

输出结果如下（记住你会得到不同的值）：

99.0% 71.43%
90.0% 94.20%
80.0% 97.81%
50.0% 97.66%
20.0% 96.89%
10.0% 94.80%
1.0% 92.02%

如果你仅用 1%的数据进行训练，你的准确率只有 71%，而如果你使用更多的数据，准确率会超过 90%。注意，准确率并不随着训练集的大小单调增加。决定训练集的大小是一个复杂的问题，涉及许多因素，可能会导致意想不到的副作用。

还有更多…

我们仅仅触及了训练和测试机器学习算法的表面。例如，监督数据集通常会被分成 3 部分，而不是 2 部分（训练、测试和交叉验证）。在训练算法时，你需要考虑许多其他问题，还有更多种类的算法。在本章中，我们试图培养理解机器学习的基本直觉，但如果你打算继续这条路线，这只是你的起点。

第十二章：使用 Dask 和 Zarr 进行并行处理

生物信息学数据集正在以指数速度增长。基于标准工具（如 Pandas）的数据分析策略假设数据集能够装入内存（尽管会有一些外部存储分析的处理），或者假设单台机器能够高效地处理所有数据。不幸的是，这对于许多现代数据集来说并不现实。

在本章中，我们将介绍两种能够处理非常大数据集和昂贵计算的库：

• Dask 是一个支持并行计算的库，可以扩展到从单台计算机到非常大的云环境和集群环境。Dask 提供了与 Pandas 和 NumPy 类似的接口，同时允许你处理分布在多台计算机上的大数据集。

• Zarr 是一个存储压缩和分块多维数组的库。正如我们将看到的，这些数组专为处理在大型计算机集群中处理的大数据集而设计，同时在需要时也能在单台计算机上处理数据。

我们的食谱将使用蚊子基因组学的数据介绍这些高级库。你应该将这段代码作为起点，帮助你走上处理大数据集的道路。大数据集的并行处理是一个复杂的话题，而这只是你旅程的开始——而非结束。

因为所有这些库对于数据分析都至关重要，如果你正在使用 Docker，它们都可以在 tiagoantao/bioinformatics_dask Docker 镜像中找到。

在本章中，我们将介绍以下食谱：

• 使用 Zarr 读取基因组学数据

• 使用 Python 多进程进行数据并行处理

• 使用 Dask 基于 NumPy 数组处理基因组数据

• 使用 dask.distributed 调度任务

使用 Zarr 读取基因组学数据

Zarr (zarr.readthedocs.io/en/stable/) 将基于数组的数据（如 NumPy）存储在磁盘和云存储的层次结构中。Zarr 用来表示数组的数据结构不仅非常紧凑，而且还支持并行读取和写入，这一点我们将在接下来的食谱中看到。在本食谱中，我们将读取并处理来自按蚊基因组 1000 基因组计划的数据（malariagen.github.io/vector-data/ag3/download.xhtml）。在这里，我们将仅进行顺序处理，以便引入 Zarr；在接下来的食谱中，我们将进行并行处理。我们的项目将计算单一染色体上所有基因位置的缺失数据。

准备工作

按照 cloud.google.com/storage/docs/gsutil_install 上提供的指导，从 gsutil 中获取按蚊 1000 基因组数据。在安装了 gsutil 后，使用以下代码行下载数据（约 2 千兆字节 (GB)）：

mkdir -p data/AG1000G-AO/
gsutil -m rsync -r \
         -x '.*/calldata/(AD|GQ|MQ)/.*' \
         gs://vo_agam_release/v3/snp_genotypes/all/AG1000G-AO/ \
         data/AG1000G-AO/ > /dev/null

我们从项目中下载了一个样本子集。下载数据后，处理它的代码可以在Chapter11/Zarr_Intro.py中找到。

如何操作…

查看以下步骤以开始：

1. 让我们首先检查一下 Zarr 文件中提供的结构：

   import numpy as np
   import zarr 
   mosquito = zarr.open('data/AG1000G-AO')
   print(mosquito.tree())
   

我们从打开 Zarr 文件开始（正如我们很快会看到的，这可能实际上并不是一个文件）。之后，我们会打印出里面可用的数据树：

/
├── 2L
│   └── calldata
│       └── GT (48525747, 81, 2) int8
├── 2R
│   └── calldata
│       └── GT (60132453, 81, 2) int8
├── 3L
│   └── calldata
│       └── GT (40758473, 81, 2) int8
├── 3R
│   └── calldata
│       └── GT (52226568, 81, 2) int8
├── X
│   └── calldata
│       └── GT (23385349, 81, 2) int8
└── samples (81,) |S24

Zarr 文件包含五个数组：四个对应蚊子的染色体——2L、2R、3L、3R和X（Y不包括在内）——另一个包含文件中包含的 81 个样本。最后一个数组包含样本名称——我们在这个文件中有 81 个样本。染色体数据由 8 位整数（int8）组成，样本名称则是字符串。

1. 现在，让我们探索2L染色体的数据。首先来看一些基本信息：

   gt_2l = mosquito['/2L/calldata/GT']
   gt_2l
   

这是输出：

<zarr.core.Array '/2L/calldata/GT' (48525747, 81, 2) int8>

我们有一个包含4852547个 SNP 和81个样本的数组。对于每个 SNP 和样本，我们有2个等位基因。

1. 现在我们来检查数据是如何存储的：

   gt_2l.info
   

输出看起来是这样的：

Name               : /2L/calldata/GT
Type               : zarr.core.Array
Data type          : int8
Shape              : (48525747, 81, 2)
Chunk shape        : (300000, 50, 2)
Order              : C
Read-only          : False
Compressor         : Blosc(cname='lz4', clevel=5, shuffle=SHUFFLE, blocksize=0)
Store type         : zarr.storage.DirectoryStore
No. bytes          : 7861171014 (7.3G)
No. bytes stored   : 446881559 (426.2M)
Storage ratio      : 17.6
Chunks initialized : 324/324

这里有很多内容需要解析，但现在我们将专注于存储类型、存储的字节数和存储比率。Store type的值是zarr.storage.DirectoryStore，所以数据并不在一个单独的文件中，而是存储在一个目录内。数据的原始大小是7.3 GB！但是 Zarr 使用压缩格式，将数据的大小压缩到426.2 17.6。

1. 让我们来看看数据是如何存储在目录中的。如果你列出AG1000G-AO目录的内容，你会发现以下结构：

   .
   ├── 2L
   │   └── calldata
   │       └── GT
   ├── 2R
   │   └── calldata
   │       └── GT
   ├── 3L
   │   └── calldata
   │       └── GT
   ├── 3R
   │   └── calldata
   │       └── GT
   ├── samples
   └── X
       └── calldata
           └── GT
   

2. 如果你列出2L/calldata/GT目录的内容，你会发现很多文件在编码该数组：

   0.0.0
   0.1.0
   1.0.0
   ...
   160.0.0
   160.1.0
   

2L/calldata/GT目录中有 324 个文件。记住，在前一步中我们有一个叫做Chunk shape的参数，它的值是(300000, 50, 2)。

Zarr 将数组拆分成多个块——这些块比加载整个数组更容易在内存中处理。每个块包含 30000x50x2 个元素。考虑到我们有 48525747 个 SNP，我们需要 162 个块来表示这些 SNP 的数量（48525747/300000 = 161.75），然后乘以 2 以表示样本的数量（81 个样本/每块 50 个 = 1.62）。因此，我们最终会得到 162*2 个块/文件。

提示

分块是一个广泛应用的技术，用于处理不能完全一次性加载到内存中的数据。这包括许多其他库，如 Pandas 或 Zarr。稍后我们将看到一个 Zarr 的例子。更大的观点是，你应该意识到分块的概念，因为它在许多需要大数据的场景中都有应用。

1. 在我们加载 Zarr 数据进行处理之前，先创建一个函数来计算一个块的基本基因组统计信息。我们将计算缺失值、祖先纯合子数量和异合子数量：

   def calc_stats(my_chunk):
       num_miss = np.sum(np.equal(my_chunk[:,:,0], -1), axis=1)
       num_anc_hom = np.sum(
           np.all([
               np.equal(my_chunk[:,:,0], 0),
               np.equal(my_chunk[:,:,0], my_chunk[:,:,1])], axis=0), axis=1)
       num_het = np.sum(
           np.not_equal(
               my_chunk[:,:,0],
               my_chunk[:,:,1]), axis=1)
       return num_miss, num_anc_hom, num_het
   

如果你查看前面的函数，你会注意到没有任何与 Zarr 相关的内容：它只是 NumPy 代码。Zarr 有一个非常轻量的应用程序接口（API），它将 NumPy 中的大多数数据暴露出来，使得如果你熟悉 NumPy，它非常容易使用。

1. 最后，让我们遍历数据——也就是遍历我们的数据块来计算统计信息：

   complete_data = 0
   more_anc_hom = 0
   total_pos = 0
   for chunk_pos in range(ceil(max_pos / chunk_pos_size)):
       start_pos = chunk_pos * chunk_pos_size
       end_pos = min(max_pos + 1, (chunk_pos + 1) * chunk_pos_size)
       my_chunk = gt_2l[start_pos:end_pos, :, :]
       num_samples = my_chunk.shape[1]
       num_miss, num_anc_hom, num_het = calc_stats(my_chunk)
       chunk_complete_data = np.sum(np.equal(num_miss, 0))
       chunk_more_anc_hom = np.sum(num_anc_hom > num_het)
       complete_data += chunk_complete_data
       more_anc_hom += chunk_more_anc_hom
       total_pos += (end_pos - start_pos)
   print(complete_data, more_anc_hom, total_pos)
   

大多数代码负责管理数据块，并涉及算术运算来决定访问数组的哪部分。就准备好的 Zarr 数据而言，重要的部分是my_chunk = gt_2l[start_pos:end_pos, :, :]这一行。当你切片 Zarr 数组时，它会自动返回一个 NumPy 数组。

提示

在将数据加载到内存时要非常小心。记住，大多数 Zarr 数组的大小将远大于你可用的内存，因此如果尝试加载，可能会导致应用程序甚至计算机崩溃。例如，如果你执行all_data = gt_2l[:, :, :]，你将需要大约 8 GB 的空闲内存来加载它——正如我们之前看到的，数据大小为 7.3 GB。

还有更多…

Zarr 具有比这里展示的更多功能，虽然我们将在接下来的示例中探索一些其他功能，但仍有一些可能性是你应该了解的。例如，Zarr 是少数几个允许并发写入数据的库之一。你还可以更改 Zarr 表示的内部格式。

正如我们在这里看到的，Zarr 能够以非常高效的方式压缩数据——这是通过使用 Blosc 库实现的（www.blosc.org/）。由于 Blosc 的灵活性，你可以更改 Zarr 数据的内部压缩算法。

另见

Zarr 有替代格式——例如，分层数据格式 5（HDF5）（en.wikipedia.org/wiki/Hierarchical_Data_Format）和网络公共数据格式（NetCDF）（en.wikipedia.org/wiki/NetCDF）。虽然这些格式在生物信息学领域之外更为常见，但它们功能较少——例如，缺乏并发写入功能。

使用 Python 多处理进行数据并行处理

处理大量数据时，一种策略是并行处理，以便利用所有可用的中央处理单元（CPU）的计算能力，因为现代计算机通常有多个核心。在理论上的最佳情况下，如果你的计算机有八个核心，你可以通过并行处理获得八倍的性能提升。

不幸的是，典型的 Python 代码只能使用一个核心。话虽如此，Python 具有内置功能来利用所有可用的 CPU 资源；事实上，Python 提供了几种方法来实现这一点。在本配方中，我们将使用内置的multiprocessing模块。这里提供的解决方案在单台计算机上运行良好，并且如果数据集能适应内存的话也没有问题，但如果你想要在集群或云端扩展，应该考虑使用 Dask，我们将在接下来的两篇配方中介绍它。

我们在这里的目标仍然是计算与缺失值和杂合度相关的统计信息。

准备工作

我们将使用与之前的配方相同的数据。该配方的代码可以在Chapter11/MP_Intro.py中找到。

如何做…

请按照以下步骤开始：

1. 我们将使用与之前的配方完全相同的函数来计算统计信息——这是一个高度依赖 NumPy 的函数：

   import numpy as np
   import zarr
   def calc_stats(my_chunk):
       num_miss = np.sum(np.equal(my_chunk[:,:,0], -1), axis=1)
       num_anc_hom = np.sum(
           np.all([
               np.equal(my_chunk[:,:,0], 0),
               np.equal(my_chunk[:,:,0], my_chunk[:,:,1])], axis=0), axis=1)
       num_het = np.sum(
           np.not_equal(
               my_chunk[:,:,0],
               my_chunk[:,:,1]), axis=1)
       return num_miss, num_anc_hom, num_het
   

2. 让我们访问我们的蚊子数据：

   mosquito = zarr.open('data/AG1000G-AO')
   gt_2l = mosquito['/2L/calldata/GT']
   

3. 尽管我们使用相同的函数来计算统计信息，但我们对整个数据集的处理方式将有所不同。首先，我们计算所有将调用calc_stats的区间。这些区间将被设计成与变异体的块划分完美匹配：

   chunk_pos_size = gt_2l.chunks[0]
   max_pos = gt_2l.shape[0]
   intervals = []
   for chunk_pos in range(ceil(max_pos / chunk_pos_size)):
       start_pos = chunk_pos * chunk_pos_size
       end_pos = min(max_pos + 1, (chunk_pos + 1) * chunk_pos_size)
       intervals.append((start_pos, end_pos))
   

我们的区间列表必须与磁盘上的块划分相关。这项计算会很高效，只要这个映射尽可能接近。

1. 现在，我们将把计算每个区间的代码分离到一个函数中。这一点很重要，因为multiprocessing模块将在它创建的每个进程中多次执行这个函数：

   def compute_interval(interval):
       start_pos, end_pos = interval
       my_chunk = gt_2l[start_pos:end_pos, :, :]
       num_samples = my_chunk.shape[1]
       num_miss, num_anc_hom, num_het = calc_stats(my_chunk)
       chunk_complete_data = np.sum(np.equal(num_miss, 0))
       chunk_more_anc_hom = np.sum(num_anc_hom > num_het)
       return chunk_complete_data, chunk_more_anc_hom
   

2. 我们现在终于将让代码在多个核心上执行：

   with Pool() as p:
       print(p)
       chunk_returns = p.map(compute_interval, intervals)
       complete_data = sum(map(lambda x: x[0], chunk_returns))
       more_anc_hom = sum(map(lambda x: x[1], chunk_returns))
       print(complete_data, more_anc_hom)
   

第一行使用multiprocessing.Pool对象创建一个上下文管理器。Pool对象默认会创建多个编号为os.cpu_count()的进程。池提供了一个map函数，能够在所有创建的进程中调用我们的compute_interval函数。每次调用将处理一个区间。

还有更多…

本配方简要介绍了如何在 Python 中进行并行处理，而无需使用外部库。话虽如此，它展示了并发并行处理的最重要构建块。

由于 Python 中的线程管理方式，线程并不是实现真正并行处理的可行替代方案。纯 Python 代码无法通过多线程并行执行。

一些你可能使用的库——通常 NumPy 就是这样的——能够在执行顺序代码时利用所有底层处理器。确保在使用外部库时，不要过度占用处理器资源：当你有多个进程时，底层库也会使用多个核心。

另见

• 有关 multiprocessing 模块的讨论还有很多。你可以从标准文档开始了解：docs.python.org/3/library/multiprocessing.xhtml

• 要理解为什么基于 Python 的多线程不能充分利用所有 CPU 资源，请阅读有关 全局解释器锁 (GIL) 的内容：realpython.com/python-gil/.

使用 Dask 处理基于 NumPy 数组的基因组数据

Dask 是一个提供高级并行处理的库，可以从单个计算机扩展到非常大的集群或云操作。它还提供了处理比内存更大的数据集的能力。它能够提供与常见 Python 库如 NumPy、Pandas 或 scikit-learn 相似的接口。

我们将重复之前配方中的一个子集——即计算数据集中 SNP 的缺失情况。我们将使用 Dask 提供的类似于 NumPy 的接口。

在我们开始之前，请注意 Dask 的语义与 NumPy 或 Pandas 等库有很大不同：它是一个懒加载库。例如，当你指定一个等效于 np.sum 的调用时，你实际上并没有计算和求和，而是在未来会计算它的任务。让我们进入配方来进一步澄清这一点。

准备就绪

我们将以一种完全不同的方式重新分块 Zarr 数据。我们这么做的原因是为了在准备算法时能够可视化任务图。包含五个操作的任务图比包含数百个节点的任务图更容易可视化。为了实际目的，你不应该像我们这里做的那样将数据重新分块为如此小的块。实际上，如果你根本不重新分块这个数据集，也是完全可以的。我们这么做只是为了可视化的目的：

import zarr
mosquito = zarr.open('data/AG1000G-AO/2L/calldata/GT')
zarr.array(
    mosquito,
    chunks=(1 + 48525747 // 4, 81, 2),
    store='data/rechunk')

我们最终会得到非常大的块，虽然这对我们的可视化目的很有用，但它们可能太大而无法放入内存中。

此配方的代码可以在 Chapter11/Dask_Intro.py 中找到。

如何实现…

1. 让我们首先加载数据并检查 DataFrame 的大小：

   import numpy as np
   import dask.array as da
   
   mosquito = da.from_zarr('data/rechunk')
   mosquito
   

如果你在 Jupyter 中执行，这将是输出结果：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_11_001.jpg

图 11.1 - Dask 数组的 Jupyter 输出，汇总我们的数据

完整的数组占用 7.32 GB。最重要的数字是块的大小：1.83 GB。每个工作节点需要有足够的内存来处理一个块。记住，我们这里只是使用了较少的块数，以便能够在这里绘制任务。

由于大块数据的大小，我们最终只得到了四个块。

我们尚未将任何内容加载到内存中：我们只是指定了最终想要执行的操作。我们正在创建一个任务图来执行，而不是立即执行——至少目前如此。

1. 让我们来看看我们需要执行哪些任务来加载数据：

   mosquito.visualize()
   

这是输出结果：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_11_002.jpg

图 11.2 - 加载我们的 Zarr 数组所需执行的任务

因此，我们有四个任务要执行，每个块对应一个任务。

1. 现在，让我们看看计算每个块缺失值的函数：

   def calc_stats(variant):
       variant = variant.reshape(variant.shape[0] // 2, 2)
       misses = np.equal(variant, -1)
       return misses
   

每个块的函数将在 NumPy 数组上操作。请注意区别：我们在主循环中使用的代码是针对 Dask 数组的，但在块级别，数据以 NumPy 数组的形式呈现。因此，这些块必须适配内存，因为 NumPy 需要如此。

1. 后面，当我们实际使用这个函数时，我们需要一个二维（2D）数组。由于数组是三维（3D）的，我们需要对数组进行重塑：

   mosquito_2d = mosquito.reshape(
       mosquito.shape[0],
       mosquito.shape[1] * mosquito.shape[2])
   mosquito_2d.visualize()
   

这是当前的任务图：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_11_003.png

图 11.3 - 加载基因组数据并重塑的任务图

reshape 操作发生在 dask.array 层，而不是 NumPy 层，因此它仅向任务图中添加了节点。仍然没有执行。

1. 现在，让我们准备执行这个函数——意味着在整个数据集上向我们的任务图中添加任务。有很多种执行方式；在这里，我们将使用 dask.array 提供的 apply_along_axis 函数，它基于 NumPy 中同名的函数：

   max_pos = 10000000
   stats = da.apply_along_axis(
       calc_stats, 1, mosquito_2d[:max_pos,:],
       shape=(max_pos,), dtype=np.int64)
   stats.visualize()
   

目前，我们只打算研究前百万个位置。正如你在任务图中看到的，Dask 足够智能，只会对参与计算的块添加操作：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_11_004.jpg

图 11.4 - 包括统计计算在内的完整任务图

1. 记住，在此之前我们还没有进行任何计算。现在是时候真正执行任务图了：

   stats = stats.compute() 
   

这将启动计算。计算的具体方式是我们将在下一个配方中讨论的内容。

警告

由于块大小的问题，这段代码可能会导致你的计算机崩溃。至少需要 16 GB 内存才能保证安全。记住，你可以使用更小的块大小——而且你应该使用更小的块大小。我们之所以使用这样的块大小，是为了能够生成前面展示的任务图（否则，它们可能会有数百个节点，无法打印出来）。

还有更多内容…

我们没有在这里讨论如何优化 Dask 代码的策略——那将是另一本书的内容。对于非常复杂的算法，你需要进一步研究最佳实践。

Dask 提供的接口类似于其他常见的 Python 库，如 Pandas 或 scikit-learn，可以用于并行处理。你也可以将它用于不依赖现有库的通用算法。

参见

• 对于 Dask 的最佳实践，最好的起点是 Dask 文档本身，尤其是 docs.dask.org/en/latest/best-practices.xhtml。

使用 dask.distributed 调度任务

Dask 在执行方面非常灵活：我们可以在本地执行、在科学集群上执行，或者在云上执行。这种灵活性是有代价的：它需要被参数化。有多种配置 Dask 调度和执行的方式，但最通用的是 dask.distributed，因为它能够管理不同种类的基础设施。因为我不能假设你能够访问像 dask.distributed 这样存在于不同平台上的集群或云服务。

在这里，我们将再次计算 Anopheles 1000 基因组项目的不同变体的简单统计数据。

准备工作

在开始使用dask.distributed之前，我们需要注意，Dask 有一个默认的调度器，这个调度器实际上会根据你所使用的库而有所变化。例如，以下是我们 NumPy 示例的调度器：

import dask
from dask.base import get_scheduler
import dask.array as da
mosquito = da.from_zarr('data/AG1000G-AO/2L/calldata/GT')
print(get_scheduler(collections=[mosquito]).__module__)

输出将如下所示：

dask.threaded

Dask 在这里使用了一个线程调度器。对于 NumPy 数组来说，这样做是有道理的：NumPy 实现本身是多线程的（真正的多线程，带有并行性）。当底层库并行运行时，我们不希望有大量进程在后台运行。如果你使用的是 Pandas DataFrame，Dask 可能会选择一个多进程调度器。因为 Pandas 本身不支持并行，所以让 Dask 自己并行运行是有意义的。

好的——既然我们已经解决了这个重要细节，现在让我们回到环境准备工作。

dask.distributed 有一个集中式调度器和一组工作节点，我们需要启动它们。可以在命令行中运行以下代码来启动调度器：

dask-scheduler --port 8786 --dashboard-address 8787

我们可以在与调度器相同的机器上启动工作节点，方法如下：

dask-worker --nprocs 2 –nthreads 1 127.0.0.1:8786

我指定了每个进程使用一个线程。对于 NumPy 代码来说，这个配置是合理的，但实际配置将取决于你的工作负载（如果你在集群或云上，配置可能完全不同）。

小贴士

你实际上不需要像我这里所做的那样手动启动整个进程。dask.distributed 会为你启动一些东西——虽然它不会完全优化你的工作负载——如果你没有自己准备好系统（详情请见下一部分）。但我想给你一个概念，因为在很多情况下，你必须自己设置基础设施。

同样，我们将使用第一部分食谱中的数据。请确保按照准备工作部分的说明下载并准备好数据。我们不会使用重新分块的部分——我们将在下一部分的 Dask 代码中进行处理。我们的代码可以在Chapter11/Dask_distributed.py中找到。

如何做到这一点…

按照以下步骤开始：

1. 让我们从连接到之前创建的调度器开始：

   import numpy as np
   import zarr
   import dask.array as da
   from dask.distributed import Client
   
   client = Client('127.0.0.1:8786')
   client
   

如果你使用的是 Jupyter，你将看到一个很好的输出，汇总了你在此食谱的准备工作部分所创建的配置：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_11_005.jpg

图 11.5 - 使用 dask.distributed 时的执行环境摘要

你会注意到这里提到了一个仪表板。dask.distributed 提供了一个实时仪表板，允许你跟踪计算的状态。你可以在浏览器中输入 http://127.0.0.1:8787/ 来访问它，或者直接点击 图 11.5 中提供的链接。

由于我们还没有进行任何计算，仪表板大部分是空的。一定要探索顶部的许多标签：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_11_006.jpg

图 11.6 - dask.distributed 仪表板的初始状态

1. 让我们加载数据。更严格地说，让我们准备任务图以加载数据：

   mosquito = da.from_zarr('data/AG1000G-AO/2L/calldata/GT')
   mosquito
   

以下是在 Jupyter 上的输出：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_11_007.jpg

图 11.7 - 原始 Zarr 数组（2L 染色体）的汇总

1. 为了方便可视化，让我们再次进行分块。我们还将为第二个维度——样本——创建一个单一的块。这是因为我们缺失值的计算需要所有样本，而在我们的特定情况下，为每个样本创建两个块是没有意义的：

   mosquito = mosquito.rechunk((mosquito.shape[0]//8, 81, 2))
   

提醒一下，我们有非常大的块，你可能会遇到内存问题。如果是这样，你可以使用原始的块进行运行。只是可视化效果将无法读取。

1. 在继续之前，让我们要求 Dask 不仅执行重新分块操作，还要确保结果已经准备好并存储在工作节点中：

   mosquito = mosquito.persist()
   

persist 调用确保数据在工作节点中可用。在以下截图中，你可以看到计算过程中的仪表板。你可以查看每个节点上正在执行的任务、已完成和待完成的任务摘要，以及每个工作节点上存储的字节数。需要注意的是 溢写到磁盘 的概念（见屏幕左上角）。如果内存不足以容纳所有块，它们会暂时写入磁盘：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_11_008.jpg

图 11.8 - 执行持久化函数以重新分块数组时的仪表板状态

1. 现在，让我们计算统计信息。对于最后一个配方，我们将使用不同的方法——我们将请求 Dask 对每个块应用一个函数：

   def calc_stats(my_chunk):
       num_miss = np.sum(
           np.equal(my_chunk[0][0][:,:,0], -1),
           axis=1)
       return num_miss
   stats = da.blockwise(
       calc_stats, 'i', mosquito, 'ijk',
       dtype=np.uint8)
   stats.visualize()
   

请记住，每个块不是 dask.array 实例，而是一个 NumPy 数组，因此代码是在 NumPy 数组上运行的。以下是当前的任务图。没有加载数据的操作，因为之前执行的函数已经完成了所有这些操作：

https://github.com/OpenDocCN/freelearn-ds-zh/raw/master/docs/bioinfo-py-cb/img/B17942_11_009.jpg

图 11.9 - 从持久化数据开始的每个块对 calc_stats 函数的调用

1. 最后，我们可以得到我们的结果：

   stat_results = stats.compute()
   

还有更多…

关于 dask.distributed 接口，还有很多内容可以进一步讲解。在这里，我们介绍了其架构的基本概念和仪表盘。

dask.distributed 提供了基于 Python 标准 async 模块的异步接口。由于本章内容的介绍性性质，我们不会详细讨论它，但建议你查看相关内容。

另见

• 你可以从 dask.distributed 的文档开始，访问 distributed.dask.org/en/stable/。

• 在许多情况下，你需要将代码部署到集群或云端。请查看部署文档，了解不同平台的资源：docs.dask.org/en/latest/deploying.xhtml。

• 在掌握本章内容后，下一步可以学习 Python 中的异步计算。请查看 docs.python.org/3/library/asyncio-task.xhtml。