Python在生物信息学领域的应用秘籍（第一篇）

原文：Bioinformatics with Python Cookbook

协议：CC BY-NC-SA 4.0

零、前言

生物信息学是一个活跃的研究领域，它使用一系列简单到高级的计算从生物数据中提取有价值的信息，这本书将向您展示如何使用 Python 管理这些任务。

这本更新版的《Python 生物信息学指南》从 Python 生态系统中各种工具和库的快速概述开始，这些工具和库将帮助您转换、分析和可视化生物数据集。随着章节的深入，您将借助真实世界的例子，涵盖下一代测序、单细胞分析、基因组学、宏基因组学、群体遗传学、系统发育学和蛋白质组学的关键技术。您将学习如何使用重要的管道系统，如 Galaxy 服务器和 Snakemake，并了解 Python 中用于函数式和异步编程的各种模块。这本书还将帮助您探索一些主题，如在高性能计算框架下使用统计方法发现 SNP，包括 Dask 和 Spark，以及机器学习算法在生物信息学中的应用。

在这本生物信息学 Python 书籍结束时，你将具备实现最新编程技术和框架的知识，使你能够处理各种规模的生物信息学数据。

这本书是给谁的

这本书是为想要解决中高级生物学和生物信息学问题的生物信息学分析师、数据科学家、计算生物学家、研究人员和 Python 开发者而写的。Python 编程语言的工作知识是预期的。生物学基础知识会有帮助。

这本书涵盖了什么

第一章 ，Python 及周边软件生态，告诉你如何用 Python 搭建现代生物信息学环境。本章讨论了如何使用 Docker 部署软件、与 R 接口以及与 Jupyter 笔记本交互。

第二章 ，了解 NumPy、pandas、Arrow、Matplotlib ，介绍数据科学的基础 Python 库:数组和矩阵处理的 NumPy；Pandas 用于基于表格的数据操作；优化 Pandas 处理和 Matplotlib 图表的箭头。

第三章 ，下一代测序，提供处理下一代测序数据的具体解决方案。本章教你如何处理大的 FASTQ、BAM 和 VCF 文件。它还讨论了数据过滤。

第四章 ，高级 NGS 处理，涵盖了过滤 NGS 数据的高级编程技术。这包括使用孟德尔数据集，然后通过标准统计进行分析。我们还介绍了宏基因组分析

第五章 ，与基因组一起工作，不仅处理高质量的参照——如人类基因组——还讨论了如何分析非模式物种中典型的其他低质量参照。它介绍 GFF 处理，教你分析基因组特征信息，并讨论如何使用基因本体论。

第六章 ，群体遗传学，描述如何对经验数据集进行群体遗传学分析。例如，在 Python 中，我们可以执行主成分分析、计算机 FST 或结构/混合物图。

第七章 ，种系发生学，使用最近测序的埃博拉病毒的完整序列进行真正的种系发生分析，包括树重建和序列比较。本章讨论处理树状结构的递归算法。

第八章 ，利用蛋白质数据库，着重于处理 PDB 文件，例如，执行蛋白质的几何分析。本章着眼于蛋白质可视化。

第九章 ，生物信息管道，介绍两种管道。第一种管道是基于 Python 的 Galaxy，这是一个广泛使用的系统，其 web 界面主要面向非编程用户，尽管生物信息学家可能仍然需要通过编程与它进行交互。第二种将基于 snakemake 和 nextflow，这是一种面向程序员的管道。

第十章 ，生物信息学的机器学习，介绍了机器学习使用直观的方法来处理计算生物学问题。本章包括主成分分析、聚类、决策树和随机森林。

第十一章 ，用 Dask 和 Zarr 进行并行处理，介绍处理超大型数据集和计算密集型算法的技术。本章将解释如何在许多计算机(集群或云)上使用并行计算。我们还将讨论生物数据的有效存储。

第十二章 ，生物信息学的函数式编程，介绍了允许开发更复杂的 Python 程序的函数式编程，通过懒惰编程和不变性，这些程序更容易部署在具有复杂算法的并行环境中

为了充分利用这本书

| **书中涵盖的软件/硬件** | **操作系统要求** | | Python 3.9 | Windows、Mac OS X 和 Linux(首选) | | Numpy 熊猫 Matplolib | | | 生物 python | | | Dask，zarr，scikit-learn | |

如果你使用的是这本书的数字版本，我们建议你自己输入代码或者通过 GitHub 库获取代码(链接见下一节)。这样做将帮助您避免任何与复制和粘贴代码相关的潜在错误。

下载示例代码文件

你可以从 GitHub 的 https://GitHub . com/packt publishing/Bioinformatics-with-Python-Cookbook-third-edition 下载本书的示例代码文件。如果代码有更新，它将在现有的 GitHub 库中更新。

我们在github.com/PacktPublishing/也有丰富的书籍和视频目录中的其他代码包。看看他们！

下载彩色图片

我们还提供了一个 PDF 文件，其中有本书中使用的截图/图表的彩色图像。你可以在这里下载:packt.link/3KQQO。

使用的惯例

本书通篇使用了许多文本约定。

Code in text:表示文本中的码字、数据库表名、文件夹名、文件名、文件扩展名、路径名、伪 URL、用户输入和 Twitter 句柄。下面是一个例子:“call_genotype的形状为 56，241×1，1198，2，也就是说它是标注尺寸的变型、样本、ploidy。”

代码块设置如下:

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)

当我们希望将您的注意力吸引到代码块的特定部分时，相关的行或项目以粗体显示:

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

Bold :表示一个新术语、一个重要单词或您在屏幕上看到的单词。例如，菜单或对话框中的单词出现在文本中，如下所示。这里有一个例子:“关于块列，参见 第十一章——但是你现在可以安全地忽略它。”

提示或重要注意事项

像这样出现。

章节

在这本书里，你会发现几个经常出现的标题(做好准备，怎么做…、工作原理…，还有更多…和参见。

要给出如何完成配方的明确说明，请使用以下章节:

准备就绪

本节将告诉您制作方法的内容，并介绍如何设置制作方法所需的任何软件或任何初步设置。

怎么做……

本节包含遵循配方所需的步骤。

它是如何工作的……

这一部分通常包括对前一部分发生的事情的详细解释。

还有更多……

这一部分包含了关于配方的附加信息，以使你对配方有更多的了解。

参见

这个部分提供了一些有用的链接，可以链接到食谱的其他有用信息。

取得联系

我们随时欢迎读者的反馈。

总体反馈:如果您对本书的任何方面有疑问，请在邮件主题中提及书名，并发邮件至 customercare@packtpub.com 联系我们。

勘误表:虽然我们已经尽力确保内容的准确性，但错误还是会发生。如果你在这本书里发现了一个错误，请告诉我们，我们将不胜感激。请访问 www.packtpub.com/support/errata，选择您的图书，点击勘误表提交表格链接，并输入详细信息。

盗版:如果您在互联网上遇到我们作品的任何形式的非法拷贝，如果您能提供我们的地址或网站名称，我们将不胜感激。请通过 copyright@packt.com 的联系我们，并提供材料链接。

如果你有兴趣成为一名作家:如果有你擅长的主题，并且你有兴趣写书或投稿，请访问 authors.packtpub.com。

评论

请留下评论。一旦你阅读并使用了这本书，为什么不在你购买它的网站上留下评论呢？潜在的读者可以看到并使用您不带偏见的意见来做出购买决定，我们 Packt 可以了解您对我们产品的看法，我们的作者可以看到您对他们的书的反馈。谢谢大家！

更多关于 Packt 的信息，请访问packt.com。

分享你的想法

一旦你阅读了 Python 食谱中的生物信息学，我们很想听听你的想法！请点击这里直接进入亚马逊对这本书的评论页面，并分享你的反馈。

您的评论对我们和技术社区非常重要，将有助于我们确保提供高质量的内容。

一、Python 及周边软件生态

我们将从安装本书大部分内容所需的基本软件开始。这将包括 Python 发行版、一些基本的 Python 库和外部生物信息学软件。在这里，我们也将看看 Python 之外的世界。在生物信息学和大数据领域， R 也是主要玩家；因此，您将学习如何通过 rpy2 与它交互，这是一个 Python/R 桥。此外，我们将探索 IPython 框架(通过 Jupyter Lab)可以给我们带来的优势，以便有效地与 r 接口。这一章将为我们将在本书剩余部分中执行的所有计算生物学奠定基础。

由于不同的用户有不同的需求，我们将介绍两种不同的软件安装方法。一种方法是使用 Anaconda Python(docs.continuum.io/anaconda/)发行版，另一种安装软件的方法是通过 Docker(这是一种基于共享相同操作系统内核的容器的服务器虚拟化方法；请参考 https://www.docker.com/)。这仍然会为您安装 Anaconda，但是是在一个容器中。如果您使用的是基于 Windows 的操作系统，强烈建议您考虑更换操作系统或通过 Windows 上的一些现有选项使用 Docker。在 macOS 上，虽然 Docker 也是可用的，但你也许可以在本地安装大部分软件。学习使用本地发行版(Anaconda 或其他)比 Docker 容易，但是考虑到 Python 中的包管理可能很复杂，Docker 映像提供了一定程度的稳定性。

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

用 Anaconda 安装所需的软件

用 Docker 安装所需的软件

通过rpy2与 R 接口

和 Jupyter 一起表演 R 魔术

使用 Anaconda 安装所需的基础软件

在开始之前，我们需要安装一些基本的必备软件。接下来的部分将带您了解软件以及安装它们所需的步骤。每一章和部分可能都有额外的要求——我们会在本书的过程中明确这些要求。另一种开始的方法是使用 Docker 食谱，之后一切都会通过 Docker 容器为您处理。

如果您已经在使用不同的 Python 发行版，强烈建议您考虑 Anaconda，因为它已经成为数据科学和生物信息学的事实上的标准。同样，允许你安装来自 https://bioconda.github.io/ 的发行版软件。

准备就绪

Python 可以在不同的环境上运行。例如，你可以在 Java 虚拟机 ( JVM )中使用 Python(通过 Jython 或)。NET via IronPython )。然而，在这里，我们不仅关注 Python，还关注围绕它的完整的软件生态。因此，我们将使用标准( CPython )实现，因为 JVM 和。NET 版本的存在主要是为了与这些平台的本地库进行交互。

对于我们的代码，我们将使用 Python 3.10。如果你从 Python 和生物信息学开始，任何操作系统都可以。但是在这里，我们主要关心的是中级到高级的用法。因此，虽然你可能会使用 Windows 和 macOS，但大多数的重型分析将在 Linux 上完成(可能是在 ?? 的高性能计算集群或 ?? 的高性能计算集群上)。下一代测序 ( NGS )数据分析和复杂的机器学习大多在 Linux 集群上进行。

如果你在 Windows 上，你应该考虑升级到 Linux 来进行你的生物信息学工作，因为大多数现代生物信息学软件不能在 Windows 上运行。请注意，macOS 将适用于几乎所有的分析，除非您计划使用计算机集群，这可能是基于 Linux 的。

如果您使用的是 Windows 或 macOS，并且不容易访问 Linux，请不要担心。现代虚拟化软件(如 VirtualBox 和 Docker )将会拯救你，它将允许你在你的操作系统上安装一个虚拟 Linux。如果您正在使用 Windows，并且决定进行本地化而不使用 Anaconda，那么请谨慎选择库；如果您为所有东西(包括 Python 本身)安装 32 位版本，您可能会更安全。

注意

如果您使用的是 Windows，许多工具将不可用。

生物信息学和数据科学正以极快的速度发展；这不仅仅是炒作，这是现实。安装软件库时，选择版本可能会很棘手。根据您拥有的代码，它可能无法与一些旧版本一起工作，或者甚至无法与新版本一起工作。希望您使用的任何代码都能指出正确的依赖关系——尽管这并不能保证。在本书中，我们将修正所有软件包的精确版本，并且我们将确保代码将与它们一起工作。很自然，代码可能需要与其他包版本一起调整。

为这本书开发的软件可以从 https://github . com/packt publishing/Bioinformatics-with-Python-Cookbook-third-edition 获得。要访问它，您需要安装 Git。习惯 Git 可能是个好主意，因为许多科学计算软件都是用它开发的。

在正确安装 Python 堆栈之前，您需要安装您将与之交互操作的所有外部非 Python 软件。每一章的列表都会有所不同，所有特定章节的包都会在它们各自的章节中解释。幸运的是，从这本书的前几版开始，大多数生物信息学软件已经可以通过 Bioconda 项目获得；因此，安装通常很容易。

你将需要安装一些开发编译器和库，这些都是免费的。在 Ubuntu 上，考虑安装 build-essential 包(apt-get install build-essential)，在 macOS 上，考虑Xcode(developer.apple.com/xcode/)。

在下表中，您会发现使用 Python 开发生物信息学最重要的软件列表:

| **名称** | **应用** | **网址** | **目的** | | 木星计划 | 所有章节 | https://jupyter.org/ | 交互式计算 | | 熊猫 | 所有章节 | https://pandas.pydata.org/ | 数据处理 | | NumPy | 所有章节 | http://www.numpy.org/ | 阵列/矩阵处理 | | 我的天啊 | 所有章节 | https://www.scipy.org/ | 科学计算 | | 生物 python | 所有章节 | https://biopython.org/ | 生物信息学图书馆 | | 希伯恩 | 所有章节 | http://seaborn.pydata.org/ | 统计图表库 | | 稀有 | 生物信息学和统计学 | https://www.r-project.org/ | 统计计算语言 | | rpy2 | r 连通性 | https://rpy2.readthedocs.io | r 接口 | | PyVCF | 鼻胃管吸出 | https://pyvcf.readthedocs.io | VCF 加工 | | 我留下来 | 鼻胃管吸出 | https://github.com/pysam-developers/pysam | SAM/BAM 处理 | | HTSeq | NGS/基因组 | https://htseq.readthedocs.io | NGS 加工 | | 树木 | 系统发育学 | https://dendropy.org/ | 系统发育学 | | PyMol | 蛋白基因组学 | https://pymol.org | 分子可视化 | | scikit-learn | 机器学习 | http://scikit-learn.org | 机器学习库 | | Cython | 大数据 | http://cython.org/ | 高性能 | | Numba | 大数据 | https://numba.pydata.org/ | 高性能 | | Dask | 大数据 | http://dask.pydata.org | 并行处理 |

图 1.1-显示在生物信息学中有用的各种软件包的表格

我们将使用pandas来处理大部分表格数据。另一种选择是只使用标准 Python。pandas在数据科学中已经变得如此普遍，以至于用它来处理所有的表格数据可能是有意义的(如果它适合内存的话)。

我们所有的工作都将在 Jupyter 项目中进行，即 Jupyter 实验室。Jupyter 已经成为编写交互式数据分析脚本的事实上的 ?? 标准。不幸的是，Jupyter 笔记本的默认格式是基于 JSON 的。这种格式难以阅读，难以比较，并且需要导出以输入到普通的 Python 解释器中。为了避免这个问题，我们将使用jupytext(jupytext.readthedocs.io/)扩展 Jupyter，这允许我们将 Jupyter 笔记本保存为普通的 Python 程序。

怎么做…

为了让开始，请看下面的步骤:

1. 从从 https://www.anaconda.com/products/individual.下载 Anaconda 发行版开始，我们将使用 21.05 版本，尽管您可能对最新版本没有问题。您可以接受所有安装的默认设置，但是您可能希望确保conda二进制文件在您的路径中(不要忘记打开一个新窗口以便可以更新路径)。如果您有另一个 Python 发行版，请小心使用您的PYTHONPATH和现有的 Python 库。最好还是不要设置你的PYTHONPATH。尽可能卸载所有其他 Python 版本和已安装的 Python 库。

2. 我们去图书馆吧。我们现在将使用biopython=1.70创建一个名为bioinformatics_base的新conda环境，如下面的命令所示:

   conda create -n bioinformatics_base python=3.10
   

3. 让我们激活环境，如下:

   conda activate bioinformatics_base
   

4. 让我们将bioconda和conda-forge频道添加到我们的信号源列表:

   conda config --add channels bioconda
   conda config --add channels conda-forge
   

5. 另外，安装基本软件包:

   conda install \
   biopython==1.79 \
   jupyterlab==3.2.1 \
   jupytext==1.13 \
   matplotlib==3.4.3 \
   numpy==1.21.3 \
   pandas==1.3.4 \
   scipy==1.7.1
   

6. 现在，让我们保存我们的环境，以便我们稍后可以重用它来在其他机器上创建新的环境，或者如果您需要清理基础环境:

   conda list –explicit > bioinformatics_base.txt
   

7. 我们甚至可以从conda :

   conda install rpy2 r-essentials r-gridextra
   

   安装 R

注意r-essentials安装了很多 R 包，包括我们后面会用到的 ggplot2。此外，我们安装了r-gridextra，因为我们将在笔记本中使用它。

还有更多…

如果您不想使用 Anaconda，您可以使用您选择的任何发行版通过pip安装许多 Python 库。你可能需要相当多的编译器和构建工具——不仅是 C 编译器，还有 C++和 Fortran。

我们将不使用我们在前面步骤中创建的环境。相反，我们将把它作为一个基础来克隆工作环境。这是因为使用 Python 进行环境管理——即使有了conda包系统的帮助——仍然会非常痛苦。因此，我们将创造一个干净的环境，如果我们的开发环境变得难以管理，我们永远不会破坏这个环境。

比如，假设你想用scikit-learn创造一个机器学习的环境。您可以执行以下操作:

1. 使用以下内容创建原始环境的克隆:

   conda create -n scikit-learn --clone bioinformatics_base
   

2. 添加scikit-learn :

   conda activate scikit-learn
   conda install scikit-learn
   

在 JupyterLab 中，我们应该用笔记本打开 jupytext 文件，而不是文本编辑器。由于 jupytext 文件与 Python 文件具有相同的扩展名——这是一个特性，而不是一个错误——默认情况下，JupyterLab 将使用普通的文本编辑器。当我们打开一个 jupytext 文件时，我们需要覆盖缺省值。打开后右键选择笔记本，如下截图所示:

图 1.2–在笔记本中打开一个 jupytext 文件

我们的 jupytext 文件不会保存图形输出，这对于本书来说已经足够了。如果你想要一个带有图像的版本，可以使用配对笔记本。更多详情，请查看 Jupytext 页面(【https://github.com/mwouts/jupytext】??)。

警告

由于我们的代码是要在 Jupyter 中运行的，所以在本书中，我不会多次使用print来输出内容，因为单元格的最后一行会自动呈现。如果你不用笔记本，记得做一个print。

用 Docker 安装所需软件

Docker 是实现操作系统级虚拟化最广泛使用的框架。这项技术允许你拥有一个独立的容器:一个比虚拟机更轻的层，但仍然允许你划分软件。这基本上隔离了所有的进程，让人感觉每个容器都是一个虚拟机。

Docker 在开发领域的两个极端都能很好地工作:它是为学习目的设置本书内容的一种便利方式，并且可以成为您在复杂环境中部署应用程序的首选平台。这个配方是前一个配方的替代方案。

然而，对于长期的开发环境来说，遵循前面的方法可能是您的最佳途径，尽管它可能需要更费力的初始设置。

准备就绪

如果你在 Linux 上，你要做的第一件事就是安装 Docker。最安全的解决方案是从 https://www.docker.com/获得最新版本。虽然您的 Linux 发行版可能有一个 Docker 包，但它可能太旧了，而且有很多错误。

如果你是 Windows 或者 macOS 上的，不要绝望；看看 Docker 网站。有各种各样的选择可以拯救你，但没有明确的公式，因为 Docker 在这些平台上进步很快。运行我们的 64 位虚拟机需要一台相当新的计算机。如果您有任何问题，请重新启动您的机器，并确保 BIOS、VT-X 或 AMD-V 已启用。至少，你需要 6 GB 的内存，最好更多。

注意

这将需要从互联网上下载大量内容，因此请确保您有足够的带宽。还有，做好长时间等待的准备。

怎么做…

要开始，请按照下列步骤操作:

1. 在 Docker shell 上使用以下命令:

   docker build -t bio https://raw.githubusercontent.com/PacktPublishing/Bioinformatics-with-Python-Cookbook-third-edition/main/docker/main/Dockerfile
   

在 Linux 上，您需要拥有 root 权限或者被添加到 Docker Unix 组。

1. 现在您已经准备好运行容器了，如下所示:

   docker run -ti -p 9875:9875 -v YOUR_DIRECTORY:/data bio
   

2. 用操作系统上的目录替换YOUR_DIRECTORY。这将在您的主机操作系统和 Docker 容器之间共享。YOUR_DIRECTORY将在/data的容器中看到，反之亦然。

-p 9875:9875将容器的 TCP 端口9875暴露在主机端口9875上。

尤其是在 Windows 上(也许在 macOS 上)，确保你的目录在 Docker shell 环境中是可见的。如果没有，请查看 Docker 官方文档，了解如何公开目录。

1. 现在，您可以使用该系统了。将你的浏览器指向http://localhost:9875，你应该会得到 Jupyter 环境。

如果这个在 Windows 上不工作，检查官方 Docker 文档(docs.docker.com/)关于如何暴露端口。

参见

以下也值得了解:

Docker 是使用最广泛的容器化软件，最近使用量有了巨大的增长。你可以在 https://www.docker.com/了解更多。

出于安全考虑，Docker 的替代品是 rkt ，在coreos.com/rkt/可以找到。

如果你不能使用 Docker，例如，如果你没有必要的权限，就像大多数计算机集群的情况一样，那么看看 https://www.sylabs.io/singularity/ ?? 的奇点。

通过 rpy2 与 R 接口

如果有一些你需要的功能，而你在 Python 库中找不到，你的第一个调用端口是检查它的是否在 R 中实现了，对于统计方法，R 仍然是最完整的框架；此外，一些生物信息学功能在 R 中只有才有，并且可能作为属于 Bioconductor 项目的软件包提供。

rpy2 提供了一个从 Python 到 r 的声明性接口。正如你将看到的，你将能够编写非常优雅的 Python 代码来执行接口过程。为了显示接口(并尝试最常见的 R 数据结构之一 DataFrame 和最流行的 R 库之一ggplot2)，我们将从人类 1000 基因组计划(【http://www.1000genomes.org/】)下载它的元数据。这不是一本关于 R 的书，但是我们想提供一些有趣且实用的例子。

准备就绪

你需要从 1，000 个基因组序列索引中获取元数据文件。请查看github . com/packt publishing/Bioinformatics-with-Python-Cookbook-third-edition/blob/main/datasets . py，下载sequence.index文件。如果你用的是 Jupyter 笔记本，打开Chapter01/Interfacing_R.py文件，简单地执行上面的wget命令。

这个文件包含项目中所有 FASTQ 文件的信息(在接下来的章节中，我们将使用人类 1，000 基因组项目中的数据)。这包括 FASTQ 文件、样品 ID、原始群体和每个泳道的重要统计信息，如读取次数和读取的 DNA 碱基数。

要设置 Anaconda，可以运行以下命令:

conda create -n bioinformatics_r --clone bioinformatics_base
conda activate bioinformatics_r
conda install r-ggplot2=3.3.5 r-lazyeval r-gridextra rpy2

使用 Docker，您可以运行以下内容:

docker run -ti -p 9875:9875 -v YOUR_DIRECTORY:/data tiagoantao/bioinformatics_r

现在我们可以开始了。

怎么做…

要开始使用,请遵循以下步骤:

1. 我们先来做一些导入:

   import os
   from IPython.display import Image
   import rpy2.robjects as robjects
   import rpy2.robjects.lib.ggplot2 as ggplot2
   from rpy2.robjects.functions import SignatureTranslatedFunction
   import pandas as pd
   import rpy2.robjects as ro
   from rpy2.robjects import pandas2ri
   from rpy2.robjects import local_converter
   

我们将在 Python 端使用pandas。r 数据帧很好地映射到pandas。

1. 我们将使用 R 的read.delim函数:

   read_delim = robjects.r('read.delim')
   seq_data = read_delim('sequence.index', header=True, stringsAsFactors=False)
   #In R:
   # seq.data <- read.delim('sequence.index', header=TRUE, stringsAsFactors=FALSE)
   

   从我们的文件中读取数据

导入后我们做的第一件事是访问read.delim R 函数，它允许您读取文件。R 语言规范允许你在对象名中加点。所以我们要把一个函数名转换成read_delim。然后，我们称函数名为 proper 请注意以下高度声明性的特性。首先，大多数原子对象，比如字符串，可以不经过转换就被传递。其次，参数名可以无缝转换(除了点号问题)。最后，对象在 Python 命名空间中可用(然而，对象实际上在 R 命名空间中不可用；我们稍后将进一步讨论这一点)。

作为参考，我已经包含了相应的 R 代码。我希望这是一个简单的转换。seq_data对象是一个数据帧。如果你知道 basic R 或pandas，你可能知道这种类型的数据结构。如果不是，那么这本质上就是一个表，也就是一系列行，其中每一列都具有相同的类型。

1. 让我们对该数据帧进行基本检查，如下:

   print('This dataframe has %d columns and %d rows' %
   (seq_data.ncol, seq_data.nrow))
   print(seq_data.colnames)
   #In R:
   # print(colnames(seq.data))
   # print(nrow(seq.data))
   # print(ncol(seq.data))
   

再次注意代码的相似性。

1. 你甚至可以使用下面的代码混合风格:

   my_cols = robjects.r.ncol(seq_data)
   print(my_cols)
   

可以直接调用 R 函数；在这种情况下，如果他们的名字中没有点，我们将调用ncol;但是，要小心。这将显示一个输出，不是 26(列数)，而是[26]，这是一个由26元素组成的向量。这是因为，默认情况下，R 中的大多数操作都返回向量。如果你想要的列数，你必须执行my_cols[0]。此外，谈到陷阱，注意 R 数组索引是从 1 开始的，而 Python 是从 0 开始的。

1. 现在，我们需要执行一些数据清理。例如，一些列应该被解释为数字，但是它们却被读为字符串:

   as_integer = robjects.r('as.integer')
   match = robjects.r.match
   my_col = match('READ_COUNT', seq_data.colnames)[0] # vector returned
   print('Type of read count before as.integer: %s' % seq_data[my_col - 1].rclass[0])
   seq_data[my_col - 1] = as_integer(seq_data[my_col - 1])
   print('Type of read count after as.integer: %s' % seq_data[my_col - 1].rclass[0])
   

match函数有点类似于 Python 列表中的index方法。正如所料，它返回一个向量，这样我们就可以提取0元素。它也是 1 索引的，所以我们在 Python 上工作时减去 1。as_integer函数将把一列转换成整数。第一次打印将显示字符串(即被"包围的值)，而第二次打印将显示数字。

1. 我们需要多按摩一下这张桌子；这方面的细节可以在笔记本里找到。这里，我们将最终把数据帧转换成 R(记住，虽然它是一个 R 对象，但它实际上在 Python 名称空间上是可见的):

   robjects.r.assign('seq.data', seq_data)
   

这个将在 R 名称空间中创建一个名为seq.data的变量，其中包含来自 Python 名称空间的数据帧的内容。请注意，在此操作之后，两个对象将是独立的(如果您更改了其中一个，它将不会反映在另一个中)。

注意

虽然您可以在 Python 上执行绘图，但 R 具有默认的内置绘图功能(这里我们将忽略)。它还有一个名为 ?? 的库，实现了图形 ?? 的 ?? 语法(一种指定统计图表的声明性语言)。

1. 关于我们基于人类 1，000 个基因组项目的具体例子，首先，我们将绘制具有中心名称分布的直方图，其中产生了所有测序泳道。为此，我们将使用ggplot :

   from rpy2.robjects.functions import SignatureTranslatedFunction
   ggplot2.theme = SignatureTranslatedFunction(ggplot2.theme, init_prm_translate = {'axis_text_x': 'axis.text.x'})
   bar = ggplot2.ggplot(seq_data) + ggplot2.geom_bar() + ggplot2.aes_string(x='CENTER_NAME') + ggplot2.theme(axis_text_x=ggplot2.element_text(angle=90, hjust=1))
   robjects.r.png('out.png', type='cairo-png')
   bar.plot()
   dev_off = robjects.r('dev.off')
   dev_off()
   

第二行有点无趣，但却是一段重要的样板代码。我们将要调用的 R 函数中的一个有一个名字中带点的参数。因为 Python 函数调用不能这样，所以我们必须将函数主题中的axis.text.x R 参数名映射到axis_text_r Python 名。我们用猴子修补它(也就是说，我们用自己的修补版本替换ggplot2.theme)。

然后，我们画出图表本身。当我们向图表添加特性时，请注意ggplot2的声明性。首先，我们指定seq_data数据帧，然后我们使用一个叫做geom_bar的柱状图。接下来，我们注释x变量(CENTER_NAME)。最后，我们通过改变主题来旋转 x 轴的文本。我们通过关闭 R 打印设备来完成此操作。

1. 现在，我们可以在 Jupyter 笔记本中打印图像:

   Image(filename='out.png')
   

生成了以下图表:

图 1.3–gg plot 2 生成的中心名称直方图，负责对来自 1，000 基因组项目的人类基因组数据的泳道进行测序

1. 作为最后一个示例，我们现在将使用人类 1，000 个基因组项目(我们将彻底使用该项目的数据摘要，可以在 第三章 、下一代测序的使用现代序列格式*配方中看到)为约鲁班人(YRI)和具有北欧和西欧血统的犹他州居民(CEU)的所有测序泳道绘制一个读计数和碱基计数散点图此外，我们对不同类型测序之间的差异感兴趣(例如，外显子组覆盖、高覆盖和低覆盖)。首先，我们生成一个只有YRI和CEU通道的数据帧，并限制最大基址和读取计数:

   robjects.r('yri_ceu <- seq.data[seq.data$POPULATION %in% c("YRI", "CEU") & seq.data$BASE_COUNT < 2E9 & seq.data$READ_COUNT < 3E7, ]')
   yri_ceu = robjects.r('yri_ceu')
   ```* 
   

2. 现在我们准备绘图:

   scatter = ggplot2.ggplot(yri_ceu) + ggplot2.aes_string(x='BASE_COUNT', y='READ_COUNT', shape='factor(POPULATION)', col='factor(ANALYSIS_GROUP)') + ggplot2.geom_point()
   robjects.r.png('out.png')
   scatter.plot()
   

希望这个例子(请参考下面的截图)能让我们清楚地了解图形方法的语法的力量。我们将从声明数据框架和使用的图表类型开始(即由geom_point实现的散点图)。

请注意，表达每个点的形状取决于POPULATION变量，颜色取决于ANALYSIS_GROUP变量是多么容易:

图 1.4–gg plot 2 生成的散点图，带有所有测序泳道读数的碱基和读数计数；每个点的颜色和形状反映了分类数据(群体和排序数据的类型)

1. 因为 R 数据帧与和pandas非常接近，所以在两者之间进行转换是有意义的，因为这是由rpy2 :

   import rpy2.robjects as ro
   from rpy2.robjects import pandas2ri
   from rpy2.robjects.conversion import localconverter 
   with localconverter(ro.default_converter + pandas2ri.converter):
     pd_yri_ceu = ro.conversion.rpy2py(yri_ceu)
   del pd_yri_ceu['PAIRED_FASTQ']
   with localconverter(ro.default_converter + pandas2ri.converter):
     no_paired = ro.conversion.py2rpy(pd_yri_ceu)
   robjects.r.assign('no.paired', no_paired)
   robjects.r("print(colnames(no.paired))")
   

   支持的

我们从导入必要的转换模块开始—rpy2提供了许多将数据从 R 转换成 Python 的策略。这里，我们关注的是数据帧转换。然后我们转换 R 数据帧(注意，我们转换的是 R 名称空间中的yri_ceu，而不是 Python 名称空间中的那个)。我们在pandas数据帧上删除指示配对 FASTQ 文件名称的列，并将其复制回 R 名称空间。如果您打印新的 R 数据帧的列名，您将会看到PAIRED_FASTQ不见了。

还有更多…

值得重复的是，Python 软件生态的进步正在以极快的速度发生。这意味着，如果某个功能目前还不可用，它可能会在不久的将来发布。因此，如果您正在开发一个新项目，在使用 R 包中的功能之前，一定要检查 Python 前沿的最新发展。

生物导体项目中有大量的生物信息学 R 包(www.bioconductor.org/)。这可能是您在生物信息学功能的 R 世界中的第一个停靠站。不过注意很多 R 生物信息学包并不在 Bioconductor 上，所以一定要在综合 R 档案网 ( CRAN )(参考 CRAN 在cran.rproject.org/)上搜索更广泛的 R 包。

Python 有很多绘图库。Matplotlib 是最常见的库，但是您还有很多其他选择。在 R 的上下文中，值得注意的是有一个类似于 ggplot2 的 Python 实现，它基于图表的图形描述语言的语法，而且——令人惊讶的是——这叫做ggplot!(yhat.github.io/ggpy/)。

参见

要了解有关这些主题的更多信息，请参考以下资源:

有大量关于 R 的教程和书籍；查看 R 网页(www.r-project.org/)获取文档。

对于生物导体，查看位于manuals.bioinformatics.ucr.edu/home/R_BioCondManual的文件。

如果您与 NGS 合作，您可能还想看看在 http://manuals.bioinformatics.ucr.edu/home/ht-seq使用 Bioconductor 进行高通量序列分析。

rpy库文档是您通向 R 的 Python 门户，可以在 https://rpy2.bitbucket.io/的找到。

斯普林格的利兰·威尔金森在一本名为《?? 图形语法》的书中描述了图形语法的方法。

在数据结构方面，可以在pandas库中找到与 R 类似的功能。你可以在 http://pandas.pydata.org/pandas-docs/dev/tutorials.xhtml 找到一些教程。奥莱利媒体公司的韦斯·麦金尼所著的《用于数据分析的 Python》一书也是一个可以考虑的选择。在下一章，我们将讨论熊猫，并在整本书中使用它。

与朱庇特一起表演魔术

与标准 Python 相比，Jupyter 提供了相当多的额外特性。在这些特性中，它提供了一个名为 magics 的可扩展命令框架(实际上，这只适用于 Jupyter 的 IPython 内核，因为它实际上是一个 IPython 特性，但才是我们所关心的)。魔法允许你以许多有用的方式扩展语言。您可以使用一些神奇的函数来处理 R。正如您将在我们的示例中看到的，它使 R 接口变得更加容易和更具声明性。这份食谱不会引入任何新的 R 功能，但希望它能清楚地表明 IPython 如何在这方面成为科学计算的一个重要生产力提升。

准备就绪

您需要遵循先前准备好的步骤通过 rpy2 配方与 R 接口。笔记本是Chapter01/R_magic.py。笔记本比这里介绍的食谱更完整，包括更多的图表示例。为了简洁起见，我们将只关注使用魔法与 R 交互的基本构造。如果您正在使用 Docker，您可以使用以下内容:

docker run -ti -p 9875:9875 -v YOUR_DIRECTORY:/data tiagoantao/bioinformatics_r

怎么做…

这个配方是对前一个配方的积极简化，因为它展示了 R magics 的简洁和优雅:

1. 你需要做的第一件事是加载 R magics 和ggplot2 :

   import rpy2.robjects as robjects
   import rpy2.robjects.lib.ggplot2 as ggplot2
   %load_ext rpy2.ipython
   

注意，%启动了一个特定于 IPython 的指令。举个简单的例子，你可以在 Jupyter 单元格上写%R print(c(1, 2))。

看看不使用robjects包执行 R 代码有多容易。其实，rpy2是被用来看引擎盖下的。

1. 让来读一下在上一个菜谱中下载的sequence.index文件:

   %%R
   seq.data <- read.delim('sequence.index', header=TRUE, stringsAsFactors=FALSE)
   seq.data$READ_COUNT <- as.integer(seq.data$READ_COUNT)
   seq.data$BASE_COUNT <- as.integer(seq.data$BASE_COUNT)
   

然后，您可以使用%%R(注意双%%)指定整个单元格应该被解释为 R 代码。

1. 我们现在可以将变量转移到 Python 名称空间:

   seq_data = %R seq.data
   print(type(seq_data))  # pandas dataframe!
   

数据帧的类型不是标准的 Python 对象，而是一个pandas数据帧。这与之前版本的 R magic 界面有所不同。

1. 因为我们有一个pandas数据框架，我们可以使用pandas接口:

   my_col = list(seq_data.columns).index("CENTER_NAME")
   seq_data['CENTER_NAME'] = seq_data['CENTER_NAME'].apply(lambda` x: x.upper())
   

   很容易地对它进行操作

2. 让我们把这个数据帧放回 R 名称空间，如下:

   %R -i seq_data
   %R print(colnames(seq_data))
   

-i参数通知 magic 系统，Python 空间后面的变量将被复制到 R 名称空间中。第二行只是显示数据帧确实在 r 中可用。我们使用的名称不同于原来的——它是seq_data,而不是seq.data。

1. 让我们做一些最后的清理(更多细节，见之前的配方)并打印和之前一样的条形图:

   %%R
   bar <- ggplot(seq_data) +  aes(factor(CENTER_NAME)) + geom_bar() + theme(axis.text.x = element_text(angle = 90, hjust = 1))
   print(bar)
   

此外，R magic 系统允许您减少代码，因为它改变了 R 与 IPython 的交互行为。例如，在之前配方的ggplot2代码中，你不需要使用.png和dev.off R 函数，因为魔法系统会为你处理这些。当你告诉 R 打印一个图表时，它会神奇地出现在你的笔记本或图形控制台上。

还有更多…

随着时间的推移，R magics 的界面似乎已经发生了很大的变化。比如这本书第一版的 R 代码我就更新过几次。DataFrame 赋值的当前版本返回pandas对象，这是一个主要的变化。

参见

有关更多信息，请查看以下链接:

有关 IPython magics 的基本说明，请参见IPython . readthedocs . io/en/stable/interactive/magics . XHTML。

IPython 的第三方扩展列表，包括神奇的扩展，可以在github.com/ipython/ipython/wiki/Extensions-Index找到。

二、了解 NumPy、pandas、Arrow 和 Matplotlib

Python 最大的优势之一是它丰富的高质量科学和数据处理库。所有这些的核心是 NumPy ，它提供了高效的数组和矩阵支持。在 NumPy 上面，我们可以找到几乎所有的科学图书馆。例如，在我们的领域，有生物制药。但是其他通用数据分析库也可以用于我们的领域。例如，熊猫是处理表格数据的事实上的标准。最近， Apache Arrow 提供了一些 pandas 功能的高效实现，以及语言互操作性。最后， Matplotlib 是 Python 空间中最常见的绘图库，适合科学计算。虽然这些是具有广泛适用性的通用库，但它们是生物信息学处理的基础，所以我们将在本章研究它们。

我们将从熊猫开始，因为它提供了一个具有广泛实用性的高级库。然后，我们将介绍 Arrow，我们将只在支持熊猫的范围内使用它。之后，我们将讨论 NumPy，它是我们所做的几乎所有事情背后的主力。最后，我们将介绍 Matplotlib。

我们的食谱是非常入门的——这些图书馆中的每一个都很容易占据一整本书，但是食谱应该足以帮助你阅读这本书。如果您正在使用 Docker，并且因为所有这些库都是数据分析的基础，它们可以在来自 第一章 的tiagoantao/bioinformatics_base Docker 图像中找到。

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

用熊猫来处理疫苗不良事件

处理加入熊猫数据框架的陷阱

减少熊猫数据帧的内存使用

用 Apache Arrow 加速 pandas 处理

理解 NumPy 是 Python 数据科学和生物信息学背后的引擎

介绍用于图表生成的 Matplotlib

利用熊猫处理疫苗不良事件

我们将用一个具体的生物信息学数据分析例子来介绍熊猫:我们将研究来自疫苗不良事件报告系统 ( VAERS 、vaers.hhs.gov/)的数据。由美国卫生与公众服务部维护的 VAERS 包括一个追溯到 1990 年的疫苗不良事件数据库。

VAERS 让数据以逗号分隔值 ( CSV )格式可用。 CSV 格式非常简单，甚至可以用简单的文本编辑器(小心非常大的文件，因为它们可能会使您的编辑器崩溃)或 Excel 等电子表格打开。熊猫可以很容易地用这种格式工作。

准备就绪

首先，我们需要下载数据。在 https://vaers.hhs.gov/data/datasets.xhtml 有售。请下载 ZIP 文件:我们将使用 2021 文件；不要只下载一个 CSV 文件。下载完文件后，解压，然后用gzip –9 *csv单独重新压缩所有文件，节省磁盘空间。

您可以使用文本编辑器随意查看这些文件，或者最好使用诸如less(对于压缩文件使用zless)之类的工具。你可以在vaers . hhs . gov/docs/VAERSDataUseGuide _ en _ September 2021 . pdf找到文件内容的文档。

如果您正在使用笔记本，它们的开头会提供代码，以便您可以负责必要的处理。如果你用的是 Docker，基本镜像就够了。

代码可以在Chapter02/Pandas_Basic.py中找到。

怎么做…

请遵循以下步骤:

1. 让我们从加载主数据文件和收集基本统计数据开始:

   vdata = pd.read_csv(
       "2021VAERSDATA.csv.gz", encoding="iso-8859-1")
   vdata.columns
   vdata.dtypes
   vdata.shape
   

我们从加载数据开始。在大多数情况下，不需要担心默认的文本编码，UTF-8，将工作，但在这种情况下，文本编码是legacy iso-8859-1。然后，我们打印列名，以VAERS_ID、RECVDATE、STATE、AGE_YRS等开始。它们包括对应于每一列的 35 个条目。然后，我们打印每一列的类型。以下是最初的几个条目:

VAERS_ID          int64
RECVDATE         object
STATE            object
AGE_YRS         float64
CAGE_YR         float64
CAGE_MO         float64
SEX              object

通过这样做，我们得到数据的形状:(654986, 35)。这意味着 654，986 行和 35 列。您可以使用前面的任何策略来获得您需要的关于表元数据的信息。

1. 现在，我们来探究一下数据:

   vdata.iloc[0]
   vdata = vdata.set_index("VAERS_ID")
   vdata.loc[916600]
   vdata.head(3)
   vdata.iloc[:3]
   vdata.iloc[:5, 2:4]
   

我们可以通过许多方式来查看数据。我们将根据位置从检查第一行开始。以下是节略版:

VAERS_ID                                       916600
RECVDATE                                       01/01/2021
STATE                                          TX
AGE_YRS                                        33.0
CAGE_YR                                        33.0
CAGE_MO                                        NaN
SEX                                            F
…
TODAYS_DATE                                          01/01/2021
BIRTH_DEFECT                                  NaN
OFC_VISIT                                     Y
ER_ED_VISIT                                       NaN
ALLERGIES                                       Pcn and bee venom

在我们通过VAERS_ID索引之后，我们可以使用一个 ID 来获得一行。我们可以使用 916600(这是前面记录中的 ID)并获得相同的结果。

然后，我们检索前三个行。请注意我们可以用两种不同的方式来做到这一点:

使用head方法

使用更通用的数组规范；也就是iloc[:3]

最后，我们检索前五行，但只检索第二和第三列–iloc[:5, 2:4]。以下是输出:

 AGE_YRS  CAGE_YR
VAERS_ID 
916600       33.0     33.0
916601       73.0     73.0
916602       23.0     23.0
916603       58.0     58.0
916604       47.0     47.0

1. 现在让我们做一些基本的计算，即计算数据集中的最大年龄:

   vdata["AGE_YRS"].max()
   vdata.AGE_YRS.max()
   

最大值为 119 年。比结果更重要的是，注意访问访问列的两种方言AGE_YRS(作为字典键和作为对象字段)。

1. 现在，让我们画出相关的年龄:

   vdata["AGE_YRS"].sort_values().plot(use_index=False)
   vdata["AGE_YRS"].plot.hist(bins=20) 
   

这将生成两个图(在下面的步骤中显示了一个压缩版本)。我们在这里使用 pandas 绘图机器，它在下面使用 Matplotib。

1. 虽然我们已经有了使用 Matplotlib 制作图表的完整方法(介绍用于图表生成的 Matplotlib)，但是让我们通过直接使用它来先睹为快:

   import matplotlib.pylot as plt
   fig, ax = plt.subplots(1, 2, sharey=True)
   fig.suptitle("Age of adverse events")
   vdata["AGE_YRS"].sort_values().plot(
       use_index=False, ax=ax[0],
       xlabel="Obervation", ylabel="Age")
   vdata["AGE_YRS"].plot.hist(bins=20, orientation="horizontal")
   

这包括前面步骤中的两个数字。以下是输出:

图 2.1–左侧–每次观察不良反应的年龄；右图——显示年龄分布的直方图

1. 我们也可以采取非图形化的、更的分析方法，比如统计每年的事件:

   vdata["AGE_YRS"].dropna().apply(lambda x: int(x)).value_counts()
   

输出如下所示:

50     11006
65     10948
60     10616
51     10513
58     10362
 ...

1. 现在，让我们看看死了多少人:

   vdata.DIED.value_counts(dropna=False)
   vdata["is_dead"] = (vdata.DIED == "Y")
   

计数的输出如下:

NaN    646450
Y        8536
Name: DIED, dtype: int64

请注意，DIED的类型是而不是布尔值。有一个布尔特征的布尔表示更具有说明性，所以我们为它创建了is_dead。

小费

这里，我们假设 NaN 被解释为False。总的来说，对南的解读要慎重。这可能意味着False或者仅仅意味着——在大多数情况下——缺乏数据。如果是那样的话，就不应该改成False。

1. 现在，让我们将死亡的个体数据与所涉及的疫苗类型联系起来:

   dead = vdata[vdata.is_dead]
   vax = pd.read_csv("2021VAERSVAX.csv.gz", encoding="iso-8859-1").set_index("VAERS_ID")
   vax.groupby("VAX_TYPE").size().sort_values()
   vax19 = vax[vax.VAX_TYPE == "COVID19"]
   vax19_dead = dead.join(vax19)
   

在我们获得仅包含死亡的数据帧后，我们必须读取包含疫苗信息的数据。首先要对疫苗的种类及其不良事件做一些探索性分析。以下是节略输出:

 …
HPV9         1506
FLU4         3342
UNK          7941
VARZOS      11034
COVID19    648723

之后，我们必须只选择与 COVID 相关的疫苗，并将它们与个人数据结合起来。

1. 最后，让我们来看看在死亡方面被过度代表的前 10 个 COVID 疫苗批次，以及美国有多少个州受到每个批次的影响:

   baddies = vax19_dead.groupby("VAX_LOT").size().sort_values(ascending=False)
   for I, (lot, cnt) in enumerate(baddies.items()):
       print(lot, cnt, len(vax19_dead[vax19_dead.VAX_LOT == lot].groupby""STAT"")))
       if i == 10:
           break
   

输出如下所示:

Unknown 254 34
EN6201 120 30
EN5318 102 26
EN6200 101 22
EN6198 90 23
039K20A 89 13
EL3248 87 17
EL9261 86 21
EM9810 84 21
EL9269 76 18
EN6202 75 18

这就结束了这个食谱！

还有更多…

前面关于疫苗和批次的数据不完全正确；我们将在下一个菜谱中介绍一些数据分析陷阱。

在介绍用于图表生成的 Matplotlib配方中，我们将介绍 Matplotlib，一个为 pandas 绘图提供后端的图表库。它是 Python 数据分析生态系统的基础组件。

参见

以下是一些可能有用的额外信息:

虽然这一章的前三个食谱足以支持你看完这本书，但网上还有很多内容可以帮助你了解熊猫。您可以从主用户指南开始，该指南可在pandas.pydata.org/docs/user_guide/index.xhtml获得。

如果您需要绘制数据，不要忘记查看指南的可视化部分，因为它特别有用:pandas . pydata . org/docs/user _ guide/visualization . XHTML。

应对加入熊猫数据框架的陷阱

前一个食谱是一个旋风之旅，介绍了熊猫，并展示了我们将在本书中使用的大多数功能。关于熊猫的详尽讨论需要一本完整的书，在这本食谱中，以及在下一本食谱中，我们将讨论影响数据分析的主题，这些主题在文献中很少讨论，但非常重要。

在这份食谱中，我们将讨论一些通过连接处理相关数据帧的陷阱:事实证明，许多数据分析错误是由于不小心连接数据而引入的。我们将在这里介绍减少此类问题的技术。

准备就绪

我们将使用与前一个配方中相同的数据，但是我们将稍微混淆一下，以便我们可以讨论典型的数据分析缺陷。我们将再次把主要不良事件表与疫苗接种表结合起来，但我们将从每个表中随机抽取 90%的数据。例如，这模拟了你只有不完整信息的情况。这是表之间的连接没有直观明显结果的许多例子之一。

通过随机抽取 90%的数据，使用以下代码准备我们的文件:

vdata = pd.read_csv("2021VAERSDATA.csv.gz", encoding="iso-8859-1")
vdata.sample(frac=0.9).to_csv("vdata_sample.csv.gz", index=False)
vax = pd.read_csv("2021VAERSVAX.csv.gz", encoding="iso-8859-1")
vax.sample(frac=0.9).to_csv("vax_sample.csv.gz", index=False)

因为这段代码涉及随机抽样，所以您将得到的结果将与这里报告的结果不同。如果你想得到同样的结果，我已经提供了我在Chapter02目录中使用的文件。该配方的代码可在Chapter02/Pandas_Join.py中找到。

怎么做…

请遵循以下步骤:

1. 让我们从做个体和疫苗表的内部连接开始:

   vdata = pd.read_csv("vdata_sample.csv.gz")
   vax = pd.read_csv("vax_sample.csv.gz")
   vdata_with_vax = vdata.join(
       vax.set_index("VAERS_ID"),
       on="VAERS_ID",
       how="inner")
   len(vdata), len(vax), len(vdata_with_vax)
   

该代码的len输出为个人数据 589，487，疫苗接种数据 620，361，连接数据 558，220。这表明一些个体和疫苗数据没有被捕获。

1. 让我们找到没有被以下连接捕获的数据:

   lost_vdata = vdata.loc[~vdata.index.isin(vdata_with_vax.index)]
   lost_vdata
   lost_vax = vax[~vax["VAERS_ID"].isin(vdata.index)]
   lost_vax
   

您将看到 56，524 行个人数据未连接，而有 62，141 行接种数据。

1. 还有其他连接数据的方法。默认方式是通过执行左外连接:

   vdata_with_vax_left = vdata.join(
       vax.set_index("VAERS_ID"),
       on="VAERS_ID")
   vdata_with_vax_left.groupby("VAERS_ID").size().sort_values()
   

左外连接确保左表上的所有行总是被表示出来。如果右边没有行，那么所有右边的列都将用None值填充。

警告

有一个警告，你应该小心。请记住，左边的表格vdata——每个VAERS_ID有一个条目。当您离开 join 时，可能会出现左侧重复多次的情况。例如，我们之前做的groupby操作显示 962303 的VAERS_ID有 11 个条目。这是正确的，但是不正确的预期并不少见，即在左侧的每行输出中仍然有一行。这是因为左连接返回 1 个或多个左条目，而上面的内连接返回 0 个或 1 个条目，有时我们希望正好有 1 个条目。确保总是根据条目的数量来测试您想要的输出。

1. 还有一个右连接。让我们将 COVID 疫苗-左表-与死亡事件-右表连接起来:

   dead = vdata[vdata.DIED == "Y"]
   vax19 = vax[vax.VAX_TYPE == "COVID19"]
   vax19_dead = vax19.join(dead.set_index("VAERS_ID"), on="VAERS_ID", how="right")
   len(vax19), len(dead), len(vax19_dead)
   len(vax19_dead[vax19_dead.VAERS_ID.duplicated()])
   len(vax19_dead) - len(dead)
   

如您所料，右连接将确保右表上的所有行都被表示出来。因此，我们最终得到 583，817 个 COVID 条目，7，670 个死条目，以及 8，624 个右连接条目。

我们还检查连接表上重复条目的数量，我们得到 954。如果我们从连接的表中减去死表的长度，我们也会得到 954。确保在进行连接时有这样的检查。

1. 最后，我们将重新讨论有问题的 COVID 批次计算，因为我们现在知道我们可能会过度计算批次:

   vax19_dead["STATE"] = vax19_dead["STATE"].str.upper()
   dead_lot = vax19_dead[["VAERS_ID", "VAX_LOT", "STATE"]].set_index(["VAERS_ID", "VAX_LOT"])
   dead_lot_clean = dead_lot[~dead_lot.index.duplicated()]
   dead_lot_clean = dead_lot_clean.reset_index()
   dead_lot_clean[dead_lot_clean.VAERS_ID.isna()]
   baddies = dead_lot_clean.groupby("VAX_LOT").size().sort_values(ascending=False)
   for i, (lot, cnt) in enumerate(baddies.items()):
       print(lot, cnt, len(dead_lot_clean[dead_lot_clean.VAX_LOT == lot].groupby("STATE")))
       if i == 10:
           break
   

请注意，我们在这里使用的策略确保了不会出现重复:首先，我们将列的数量限制在我们将要使用的数量，然后我们删除重复的索引并清空VAERS_ID。这确保了VAERS_ID、VAX_LOT对不会重复，并且没有批次与 id 相关联。

还有更多…

除了左连接、内连接和右连接，还有其他类型的连接。最值得注意的是外部连接，它确保两个表中的所有条目都有表示。

确保您的连接有测试和断言:一个非常常见的错误是对连接的行为有错误的预期。您还应该确保要连接的列上没有空值，因为它们会产生大量多余的元组。

减少熊猫数据帧的内存使用

当您处理大量信息时——例如，当分析全基因组测序数据时——内存使用可能会成为您分析的一个限制。事实证明，从记忆的角度来看，天真的熊猫并不是很有效率，我们可以大大减少它的消耗。

在这个食谱中，我们将重新访问我们的 VAERS 数据，并研究几种减少熊猫内存使用的方法。这些变化的影响可能是巨大的:在许多情况下，减少内存消耗可能意味着能够使用 pandas 或需要更替代和复杂的方法，如 Dask 或 Spark。

准备就绪

我们将使用第一个配方的数据。如果你已经运行过了，你就万事俱备了；如果没有，请遵循那里讨论的步骤。你可以在Chapter02/Pandas_Memory.py中找到这段代码。

怎么做……

请遵循以下步骤:

1. 首先，让我们加载数据并检查数据帧的大小:

   import numpy as np
   import pandas as pd
   vdata = pd.read_csv("2021VAERSDATA.csv.gz", encoding="iso-8859-1")
   vdata.info(memory_usage="deep")
   

下面是输出的一个删节版本:

RangeIndex: 654986 entries, 0 to 654985
Data columns (total 35 columns):
#   Column        Non-Null Count   Dtype 
---  ------        --------------   ----- 
0   VAERS_ID      654986 non-null  int64 
2   STATE         572236 non-null  object 
3   AGE_YRS       583424 non-null  float64
6   SEX           654986 non-null  object 
8   SYMPTOM_TEXT  654828 non-null  object 
9   DIED          8536 non-null    object 
31  BIRTH_DEFECT  383 non-null     object 
34  ALLERGIES     330630 non-null  object 
dtypes: float64(5), int64(2), object(28)
memory usage: 1.3 GB

这里，我们有关于行数和每行的类型和非空值的信息。最后，我们可以看到数据帧需要 1.3 GB 的巨大空间。

1. 我们还可以检查每一列的大小:

   for name in vdata.columns:
       col_bytes = vdata[name].memory_usage(index=False, deep=True)
       col_type = vdata[name].dtype
       print(
           name,
           col_type, col_bytes // (1024 ** 2))
   

下面是输出的一个删节版本:

VAERS_ID int64 4
STATE object 34
AGE_YRS float64 4
SEX object 36
RPT_DATE object 20
SYMPTOM_TEXT object 442
DIED object 20
ALLERGIES object 34

SYMPTOM_TEXT占用 442 MB，所以我们整个表的 1/3。

1. 现在，让我们看看DIED列。我们能找到更有效的表达方式吗？

   vdata.DIED.memory_usage(index=False, deep=True)
   vdata.DIED.fillna(False).astype(bool).memory_usage(index=False, deep=True)
   

原始列占用 21，181，488 字节，而我们的压缩表示占用 656，986 字节。那就少了 32 倍！

1. 那STATE栏呢？我们能做得更好吗？

   vdata["STATE"] = vdata.STATE.str.upper()
   states = list(vdata["STATE"].unique())
   vdata["encoded_state"] = vdata.STATE.apply(lambda state: states.index(state))
   vdata["encoded_state"] = vdata["encoded_state"].astype(np.uint8)
   vdata["STATE"].memory_usage(index=False, deep=True)
   vdata["encoded_state"].memory_usage(index=False, deep=True)
   

这里，我们将文本列STATE转换为数字列encoded_state。这个数字是州名在列表状态中的位置。我们用这个号码来查找州列表。原始列占用大约 36 MB，而编码列占用 0.6 MB。

作为这种方法的替代，你可以看看熊猫的分类变量。我更喜欢使用它们，因为它们有更广泛的应用。

1. 当我们加载数据时，我们可以应用这些优化，所以让我们为此做好准备。但是现在，我们有一个先有鸡还是先有蛋的问题:为了能够知道状态表的内容，我们必须先通过一次来获得状态列表，就像这样:

   states = list(pd.read_csv(
       "vdata_sample.csv.gz",
       converters={
          "STATE": lambda state: state.upper()
       },
       usecols=["STATE"]
   )["STATE"].unique())
   

我们有一个简单地返回大写状态的转换器。我们只返回STATE列以节省内存和处理时间。最后，我们从 DataFrame 的中得到STATE列(它只有一列)。

1. 最终的优化是而不是加载数据。想象一下，我们不需要SYMPTOM_TEXT——那是大约 1/3 的数据。那样的话，我们可以跳过它。下面是最终版本:

   vdata = pd.read_csv(
       "vdata_sample.csv.gz",
       index_col="VAERS_ID",
       converters={
          "DIED": lambda died: died == "Y",
          "STATE": lambda state: states.index(state.upper())
       },
       usecols=lambda name: name != "SYMPTOM_TEXT"
   )
   vdata["STATE"] = vdata["STATE"].astype(np.uint8)
   vdata.info(memory_usage="deep") 
   

我们现在是 714 MB，比原来的一半多一点。通过将我们用于STATE和DIED的方法应用于所有其他列，这仍然可以大大减少。

参见

以下是一些可能有用的额外信息:

如果您愿意使用一个支持库来帮助 Python 处理，请查看 Apache Arrow 上的下一个方法，它将允许您节省额外的内存以提高内存效率。

如果您最终得到的数据帧占用的内存超过了您在单台机器上可用的内存，那么您必须加快游戏的速度，使用分块技术——我们不会在 Pandas 的上下文中介绍这一点——或者可以自动处理大量数据的技术。Dask，我们将在第十一章中讨论，Dask 和 Zarr 的并行处理，允许你使用一个类似熊猫的接口来处理大于内存的数据集。

用 Apache Arrow 加速熊猫加工

当处理大量数据时，比如在全基因组测序中，pandas 既慢又耗内存。Apache Arrow 为几个 pandas 操作提供了更快、更节省内存的实现，并且可以与之互操作。

Apache Arrow 是由 pandas 的创始人 Wes McKinney 共同创建的一个项目，它有几个目标，包括以语言无关的方式处理表格数据，这允许语言互操作性，同时提供内存和计算高效的实现。这里，我们只关注第二部分:提高大数据处理的效率。我们将用一种综合的方式来对待熊猫。

这里，我们将再次使用 VAERS 数据，并展示如何使用 Apache Arrow 来加速 pandas 数据加载并减少内存消耗。

准备就绪

同样，我们将使用第一个配方的数据。确保你下载并准备了它，正如在使用熊猫处理疫苗不良事件食谱的准备部分所解释的。代码可在Chapter02/Arrow.py中找到。

怎么做…

请遵循以下步骤:

1. 让我们从使用 pandas 和 Arrow:

   import gzip
   import pandas as pd
   from pyarrow import csv
   import pyarrow.compute as pc 
   vdata_pd = pd.read_csv("2021VAERSDATA.csv.gz", encoding="iso-8859-1")
   columns = list(vdata_pd.columns)
   vdata_pd.info(memory_usage="deep") 
   vdata_arrow = csv.read_csv("2021VAERSDATA.csv.gz")
   tot_bytes = sum([
       vdata_arrow[name].nbytes
       for name in vdata_arrow.column_names])
   print(f"Total {tot_bytes // (1024 ** 2)} MB")
   

   加载数据开始

pandas 需要 1.3 GB，而 Arrow 需要 614 MB:不到内存的一半。对于像这样的大文件，这可能意味着能够在内存中处理数据或者需要找到另一种解决方案，比如 Dask。虽然 Arrow 中的一些函数与熊猫有相似的名字(例如，read_csv)，但这并不是最常见的。例如，请注意我们计算数据帧总大小的方式:通过获得每一列的大小并执行求和，这是一种不同于 pandas 的方法。

1. 让我们对推断出的类型做一个并排的比较:

   for name in vdata_arrow.column_names:
       arr_bytes = vdata_arrow[name].nbytes
       arr_type = vdata_arrow[name].type
       pd_bytes = vdata_pd[name].memory_usage(index=False, deep=True)
       pd_type = vdata_pd[name].dtype
       print(
           name,
           arr_type, arr_bytes // (1024 ** 2),
           pd_type, pd_bytes // (1024 ** 2),)
   

下面是输出的一个删节版本:

VAERS_ID int64 4 int64 4
RECVDATE string 8 object 41
STATE string 3 object 34
CAGE_YR int64 5 float64 4
SEX string 3 object 36
RPT_DATE string 2 object 20
DIED string 2 object 20
L_THREAT string 2 object 20
ER_VISIT string 2 object 19
HOSPITAL string 2 object 20
HOSPDAYS int64 5 float64 4

正如您所看到的，Arrow 通常在类型推断方面更加具体，这也是内存使用率显著降低的主要原因之一。

1. 现在，让我们做一个时间性能比较:

   %timeit pd.read_csv("2021VAERSDATA.csv.gz", encoding="iso-8859-1")
   %timeit csv.read_csv("2021VAERSDATA.csv.gz")
   

在我的电脑上，结果如下:

7.36 s ± 201 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
2.28 s ± 70.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Arrow 的实现速度快了三倍。您计算机上的结果会有所不同，因为这取决于硬件。

1. 让我们在不加载SYMPTOM_TEXT列的情况下重复内存占用比较。这是一个更公平的比较，因为大多数数值数据集往往没有很大的文本列:

   vdata_pd = pd.read_csv("2021VAERSDATA.csv.gz", encoding="iso-8859-1", usecols=lambda x: x != "SYMPTOM_TEXT")
   vdata_pd.info(memory_usage="deep")
   columns.remove("SYMPTOM_TEXT")
   vdata_arrow = csv.read_csv(
       "2021VAERSDATA.csv.gz",
        convert_options=csv.ConvertOptions(include_columns=columns))
   vdata_arrow.nbytes
   

pandas 需要 847 MB，而 Arrow 需要 205 MB:少了四倍。

1. 我们的目标是使用 Arrow 将数据加载到 pandas 中。为此，我们需要转换数据结构:

   vdata = vdata_arrow.to_pandas()
   vdata.info(memory_usage="deep")
   

这里有两点非常重要:Arrow 创建的 pandas 表示只使用了 1 GB，而 pandas 表示来自其原生的read_csv，是 1.3 GB。这意味着，即使您使用 pandas 来处理数据，Arrow 也可以创建一个更紧凑的表示。

前面的代码有一个关于内存消耗的问题:当转换器运行时，它将需要内存来保存熊猫和箭头表示，因此违背了使用更少内存的目的。Arrow 可以在创建熊猫版本时自毁其表示，从而解决问题。这条线是vdata = vdata_arrow.to_pandas(self_destruct=True)。

还有更多…

如果你有一个熊猫不能处理的非常大的数据帧，甚至在它被 Arrow 加载之后，那么 Arrow 可能会做所有的处理，因为它也有一个计算引擎。也就是说，在撰写本文时，Arrow 的引擎在功能上远不如 pandas 完善。记住 Arrow 还有许多其他的特性，比如语言互操作性，但是我们不会在本书中用到它们。

了解 NumPy 作为 Python 数据科学和生物信息学背后的引擎

您的大部分分析都会使用 NumPy，即使您没有明确地使用它。NumPy 是一个数组操作库，位于 pandas、Matplotlib、Biopython 和 scikit-learn 等库之后。虽然您的许多生物信息学工作可能不需要直接使用 NumPy，但是您应该知道它的存在，因为它支持您所做的几乎所有事情，即使只是通过其他库间接支持。

在这个菜谱中，我们将使用 VAERS 数据来展示 NumPy 是如何支持我们使用的许多核心库的。这是一个非常简单的关于这个库的介绍，所以你会意识到它的存在，并且它是几乎所有事情的幕后推手。我们的例子将从美国五个有更多不利影响的州提取病例数，将它们分成年龄组:0 到 19 岁，20 到 39 岁，直到 100 到 119 岁。

准备就绪

我们将再次使用第一个食谱中的数据，所以要确保它是可用的。它的代码可以在Chapter02/NumPy.py中找到。

怎么做……

请遵循以下步骤:

1. 让我们从装载熊猫的数据开始，减少数据，使它只与美国前五个州相关:

   import numpy as np
   import pandas as pd
   import matplotlib.pyplot as plt
   vdata = pd.read_csv(
       "2021VAERSDATA.csv.gz", encoding="iso-8859-1")
   vdata["STATE"] = vdata["STATE"].str.upper()
   top_states = pd.DataFrame({
       "size": vdata.groupby("STATE").size().sort_values(ascending=False).head(5)}).reset_index()
   top_states["rank"] = top_states.index
   top_states = top_states.set_index("STATE")
   top_vdata = vdata[vdata["STATE"].isin(top_states.index)]
   top_vdata["state_code"] = top_vdata["STATE"].apply(
       lambda state: top_states["rank"].at[state]
   ).astype(np.uint8)
   top_vdata = top_vdata[top_vdata["AGE_YRS"].notna()]
   top_vdata.loc[:,"AGE_YRS"] = top_vdata["AGE_YRS"].astype(int)
   top_states
   

顶级状态如下。这个等级将在以后用于构建一个 NumPy 矩阵:

图 2.2–美国出现负面影响最多的州

1. 现在，让我们提取两个 NumPy 数组，其中包含年龄和状态数据:

   age_state = top_vdata[["state_code", "AGE_YRS"]]
   age_state["state_code"]
   state_code_arr = age_state["state_code"].values
   type(state_code_arr), state_code_arr.shape, state_code_arr.dtype
   age_arr = age_state["AGE_YRS"].values
   type(age_arr), age_arr.shape, age_arr.dtype
   

请注意，作为熊猫基础的数据是 NumPy 数据(对两个系列的values调用返回 NumPy 类型)。此外，你可能还记得熊猫有像.shape或.dtype这样的属性:这些属性的灵感来自于 NumPy，它们的行为是一样的。

1. 现在，让我们从头开始创建一个 NumPy 矩阵(一个 2D 数组)，其中每行代表一个州，每列代表一个年龄组:

   age_state_mat = np.zeros((5,6), dtype=np.uint64)
   for row in age_state.itertuples():
       age_state_mat[row.state_code, row.AGE_YRS//20] += 1
   age_state_mat
   

该数组有五行(每个州一行)和六列(每个年龄组一列)。数组中的所有单元格必须具有相同的类型。

我们用零初始化数组。有许多方法可以初始化数组，但是如果你有一个非常大的数组，初始化它可能要花很多时间。有时，根据您的任务，数组在开始时为空可能没问题(意味着它是用随机垃圾初始化的)。那样的话，用np.empty会快很多。我们在这里使用 pandas 迭代:从 pandas 的角度来看，这不是最好的方法，但是我们希望使 NumPy 部分非常明确。

1. 我们可以很容易地提取一行数据——在我们的例子中，是一个州的数据。这同样适用于列。让我们看看加州的数据，然后是 0-19 岁年龄组:

   cal = age_state_mat[0,:]
   kids = age_state_mat[:,0]
   

注意提取行或列的语法。鉴于 pandas 从 NumPy 复制了语法，并且我们在以前的食谱中也遇到过，所以你应该对它很熟悉。

1. 现在，让我们计算一个新的矩阵，其中有每个年龄组的病例比例:

   def compute_frac(arr_1d):
       return arr_1d / arr_1d.sum()
   frac_age_stat_mat = np.apply_along_axis(compute_frac, 1, age_state_mat)
   

最后一行将compute_frac函数应用于所有行。compute_frac获取单个行并返回一个新行，其中所有元素除以总和。

1. 现在，让我们创建一个新的矩阵，作为一个百分比而不是分数——仅仅是因为它读起来更好:

   perc_age_stat_mat = frac_age_stat_mat * 100
   perc_age_stat_mat = perc_age_stat_mat.astype(np.uint8)
   perc_age_stat_mat
   

第一行只是将 2D 数组的所有元素乘以 100。Matplotlib 足够智能，可以遍历不同的数组结构。如果用任意维数的数组表示，这一行将会起作用，并且会完全按照预期的那样工作。

结果如下:

图 2.3-代表美国五个病例最多的州的疫苗不良反应分布的矩阵

1. 最后，让我们使用 Matplotlib:

   fig = plt.figure()
   ax = fig.add_subplot()
   ax.matshow(perc_age_stat_mat, cmap=plt.get_cmap("Greys"))
   ax.set_yticks(range(5))
   ax.set_yticklabels(top_states.index)
   ax.set_xticks(range(6))
   ax.set_xticklabels(["0-19", "20-39", "40-59", "60-79", "80-99", "100-119"])
   fig.savefig("matrix.png")
   

   创建矩阵的图形表示

不要在 Matplotlib 代码上花太多时间——我们将在下一个菜谱中讨论它。这里的基本要点是，您可以将 NumPy 数据结构传递给 Matplotlib。Matplotlib 和熊猫一样，都是基于 NumPy 的。

参见

以下是一些可能有用的额外信息:

NumPy 有比我们在这里讨论的更多的特性。有很多关于它们的书籍和教程。官方文档是一个很好的起点:numpy.org/doc/stable/。

NumPy 有许多重要的问题需要发现，但可能最重要的问题之一是广播:NumPy 接受不同结构的数组并正确操作的能力。详情请见numpy.org/doc/stable/user/theory.broadcasting.xhtml。

引入 Matplotlib 用于图表生成

Matplotlib 是用于生成图表的最常见的 Python 库。还有更现代的选择，比如以网络为中心的散景，但是 Matplotlib 的优势不仅在于它是最广泛可用的和广泛记录的图表库，而且在计算生物学领域，我们想要一个既以网络为中心又以纸张为中心的图表库。这是因为我们的许多图表将提交给科学杂志，这些杂志同样关注这两种格式。Matplotlib 可以为我们处理这些。

这个食谱中的许多例子也可以直接用 pandas 来完成(因此间接用 Matplotlib)，但是这里的要点是练习 Matplotlib。

我们将再次使用 VAERS 数据绘制一些关于数据帧元数据的信息，并总结流行病学数据。

准备就绪

同样，我们将使用第一个配方中的数据。代码可以在Chapter02/Matplotlib.py中找到。

怎么做…

请遵循以下步骤:

1. 我们要做的第一件事是绘制每列的空值部分:

   import numpy as np
   import pandas as pd
   import matplotlib as mpl
   import matplotlib.pyplot as plt
   vdata = pd.read_csv(
       "2021VAERSDATA.csv.gz", encoding="iso-8859-1",
       usecols=lambda name: name != "SYMPTOM_TEXT")
   num_rows = len(vdata)
   perc_nan = {}
   for col_name in vdata.columns:
       num_nans = len(vdata[col_name][vdata[col_name].isna()])
       perc_nan[col_name] = 100 * num_nans / num_rows
   labels = perc_nan.keys()
   bar_values = list(perc_nan.values())
   x_positions = np.arange(len(labels))
   

labels是我们正在分析的列名,bar_values是空值的分数，x_positions是我们接下来要绘制的条形图上的条形的位置。

1. 下面是条形图第一个版本的代码:

   fig = plt.figure()
   fig.suptitle("Fraction of empty values per column")
   ax = fig.add_subplot()
   ax.bar(x_positions, bar_values)
   ax.set_ylabel("Percent of empty values")
   ax.set_ylabel("Column")
   ax.set_xticks(x_positions)
   ax.set_xticklabels(labels)
   ax.legend()
   fig.savefig("naive_chart.png")
   

我们首先创建一个带有标题的人物对象。该图将有一个包含条形图的子图。我们还设置了几个标签，并且只使用默认值。这是令人悲伤的结果:

图 2.4–我们的第一次图表尝试，只是使用了默认值

1. 当然，我们可以做得更好。让我们格式化图表实质上更多:

   fig = plt.figure(figsize=(16, 9), tight_layout=True, dpi=600)
   fig.suptitle("Fraction of empty values per column", fontsize="48")
   ax = fig.add_subplot()
   b1 = ax.bar(x_positions, bar_values)
   ax.set_ylabel("Percent of empty values", fontsize="xx-large")
   ax.set_xticks(x_positions)
   ax.set_xticklabels(labels, rotation=45, ha="right")
   ax.set_ylim(0, 100)
   ax.set_xlim(-0.5, len(labels))
   for i, x in enumerate(x_positions):
       ax.text(
           x, 2, "%.1f" % bar_values[i], rotation=90,
           va="bottom", ha="center",
           backgroundcolor="white")
   fig.text(0.2, 0.01, "Column", fontsize="xx-large")
   fig.savefig("cleaner_chart.png")
   

我们做的第一件事是为 Matplotlib 设置一个更大的图形，以提供更紧凑的布局。我们将的 x 轴刻度标签旋转 45 度，这样它们会更合适。我们也把数值标在横条上。最后，我们没有一个标准 x 轴标签，因为它会在刻度标签的顶部。相反，我们显式地编写文本。注意，图形的坐标系可以和子图的坐标系完全不同——比如比较ax.text和fig.text的坐标。结果如下:

图 2.5–我们的第二次图表尝试，同时注意布局

1. 现在，我们将根据单个图形上的四个图对我们的数据进行一些汇总分析。我们将绘制涉及死亡的疫苗、施用和死亡之间的天数、随时间推移的死亡人数以及死亡人数最多的 10 个州的性别:

   dead = vdata[vdata.DIED == "Y"]
   vax = pd.read_csv("2021VAERSVAX.csv.gz", encoding="iso-8859-1").set_index("VAERS_ID")
   vax_dead = dead.join(vax, on="VAERS_ID", how="inner")
   dead_counts = vax_dead["VAX_TYPE"].value_counts()
   large_values = dead_counts[dead_counts >= 10]
   other_sum = dead_counts[dead_counts < 10].sum()
   large_values = large_values.append(pd.Series({"OTHER": other_sum}))
   distance_df = vax_dead[vax_dead.DATEDIED.notna() & vax_dead.VAX_DATE.notna()]
   distance_df["DATEDIED"] = pd.to_datetime(distance_df["DATEDIED"])
   distance_df["VAX_DATE"] = pd.to_datetime(distance_df["VAX_DATE"])
   distance_df = distance_df[distance_df.DATEDIED >= "2021"]
   distance_df = distance_df[distance_df.VAX_DATE >= "2021"]
   distance_df = distance_df[distance_df.DATEDIED >= distance_df.VAX_DATE]
   time_distances = distance_df["DATEDIED"] - distance_df["VAX_DATE"]
   time_distances_d = time_distances.astype(int) / (10**9 * 60 * 60 * 24)
   date_died = pd.to_datetime(vax_dead[vax_dead.DATEDIED.notna()]["DATEDIED"])
   date_died = date_died[date_died >= "2021"]
   date_died_counts = date_died.value_counts().sort_index()
   cum_deaths = date_died_counts.cumsum()
   state_dead = vax_dead[vax_dead["STATE"].notna()][["STATE", "SEX"]]
   top_states = sorted(state_dead["STATE"].value_counts().head(10).index)
   top_state_dead = state_dead[state_dead["STATE"].isin(top_states)].groupby(["STATE", "SEX"]).size()#.reset_index()
   top_state_dead.loc["MN", "U"] = 0  # XXXX
   top_state_dead = top_state_dead.sort_index().reset_index()
   top_state_females = top_state_dead[top_state_dead.SEX == "F"][0]
   top_state_males = top_state_dead[top_state_dead.SEX == "M"][0]
   top_state_unk = top_state_dead[top_state_dead.SEX == "U"][0]
   

前面的代码是严格基于 pandas 的，是为绘图活动准备的。

1. 以下代码同时绘制所有信息。我们将会有四个 2 乘 2 格式的支线剧情:

   fig, ((vax_cnt, time_dist), (death_time, state_reps)) = plt.subplots(
       2, 2,
       figsize=(16, 9), tight_layout=True)
   vax_cnt.set_title("Vaccines involved in deaths")
   wedges, texts = vax_cnt.pie(large_values)
   vax_cnt.legend(wedges, large_values.index, loc="lower left")
   time_dist.hist(time_distances_d, bins=50)
   time_dist.set_title("Days between vaccine administration and death")
   time_dist.set_xlabel("Days")
   time_dist.set_ylabel("Observations")
   death_time.plot(date_died_counts.index, date_died_counts, ".")
   death_time.set_title("Deaths over time")
   death_time.set_ylabel("Daily deaths")
   death_time.set_xlabel("Date")
   tw = death_time.twinx()
   tw.plot(cum_deaths.index, cum_deaths)
   tw.set_ylabel("Cummulative deaths")
   state_reps.set_title("Deaths per state stratified by sex") state_reps.bar(top_states, top_state_females, label="Females")
   state_reps.bar(top_states, top_state_males, label="Males", bottom=top_state_females)
   state_reps.bar(top_states, top_state_unk, label="Unknown",
                  bottom=top_state_females.values + top_state_males.values)
   state_reps.legend()
   state_reps.set_xlabel("State")
   state_reps.set_ylabel("Deaths")
   fig.savefig("summary.png")
   

我们从创建一个有 2×2 支线剧情的人物开始。subplots函数返回 figure 对象和四个轴对象，我们可以用它们来创建图表。请注意，图例位于饼图中，我们在时间距离图上使用了双轴，并且我们有一种方法来计算每个州的死亡人数图表上的堆叠条形图。结果如下:

图 2.6-总结疫苗数据的四个组合图表

还有更多…

Matplotlib 有两个界面可以使用——一个较老的界面，设计类似于 MATLAB，和一个更强大的面向对象的 ( OO )界面。尽量避免两者混淆。使用面向对象的接口可能更经得起未来的考验。类似 MATLAB 的界面在matplotlib.pyplot模块下面。让事情变得混乱的是，面向对象接口的入口点在那个模块中——也就是说，matplotlib.pyplot.figure和matplotlib.pyplot.subplots。

参见

以下是一些可能有用的额外信息:

Matplolib 的文档非常非常好。例如，有一个可视化示例库，其中包含生成每个示例的代码的链接。这可以在matplotlib.org/stable/gallery/index.xhtml找到。API 文档通常非常完整。

改善 Matplotlib 图表外观的另一种方法是使用 Seaborn 库。Seaborn 的主要目的是添加统计可视化工件，但作为副作用，当导入时，它会将 Matplotlib 的默认值更改为更容易接受的值。我们将在本书中使用 Seaborn 查看下一章提供的情节。

三、下一代测序

下一代测序 ( NGS )是本世纪生命科学的基础技术发展之一。全基因组测序 ( WGS )，限制性位点相关 DNA 测序 ( RAD-Seq )，核糖核酸测序 ( RNA-Seq )，染色质免疫沉淀测序 ( ChIP-Seq )，以及其他几种技术被常规用于研究重要的生物学问题。这些也被称为高通量测序技术，这是有充分理由的:它们产生了大量需要处理的数据。NGS 是计算生物学成为大数据学科的主要原因。最重要的是，这是一个需要强大的生物信息学技术的领域。

在这里，我们不会讨论每一个单独的 NGS 技术本身(这将需要一整本书)。我们将使用现有的 WGS 数据集——1000 个基因组项目——来说明分析基因组数据所需的最常见步骤。这里介绍的配方将很容易适用于其他基因组测序方法。其中一些也可用于转录组分析(例如 RNA-Seq)。这些配方也是独立于物种的，所以你可以将它们应用于任何其他有测序数据的物种。处理不同物种数据的最大差异与基因组大小、多样性和参考基因组的质量有关(如果你的物种存在参考基因组的话)。这些不会对 NGS 处理的自动化 Python 部分产生太大影响。无论如何，我们将在 第五章 、与基因组一起工作中讨论不同的基因组。

由于这不是一本入门书，所以希望你至少知道什么是 FASTA ( FASTA )、FASTQ、二进制比对图 ( BAM )、以及变体调用格式 ( VCF )文件。我还将使用基本的基因组术语，但不介绍它(如外显子组、非同义突变等)。您需要熟悉基本的 Python。我们将利用这些知识介绍 Python 中的基本库来执行 NGS 分析。在这里，我们将遵循标准生物信息学管道的流程。

然而，在我们深入研究来自真实项目的真实数据之前，让我们熟悉一下访问现有的基因组数据库和基本的序列处理——这是暴风雨前的一个简单开端。

如果您通过 Docker 运行内容，您可以使用tiagoantao/bioinformatics_ngs图像。如果您使用的是 Anaconda Python，本章所需的软件将在每一个菜谱中介绍。

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

访问 GenBank 并在国家生物技术信息中心的数据库中移动

执行基本序列分析

使用现代序列格式

使用路线数据

从 VCF 文件中提取数据

研究基因组可达性，过滤单核苷酸多态性 ( SNP )数据

用 HTSeq 处理 NGS 数据

访问 GenBank 并在 NCBI 数据库中移动

虽然你可能有自己的数据要分析，但你可能需要现有的基因组数据集。在这里，我们将研究如何从 NCBI 访问这样的数据库。我们不会只讨论 GenBank，还会讨论 NCBI 的其他数据库。许多人(错误地)将整套 NCBI 数据库称为 GenBank，但 NCBI 包括核苷酸数据库和许多其他数据库，例如 PubMed。

由于测序分析是一个很长的主题，并且这本书的目标是中级到高级用户，我们不会对一个核心上不太复杂的主题进行详尽的讨论。

尽管如此，这是我们将在本章末尾看到的更复杂的食谱的一个很好的热身。

准备就绪

我们会使用你在 第一章**Python 以及周边软件生态中安装的 Biopython。Biopython 提供了一个与 NCBI 提供的数据检索系统Entrez的接口。

该配方可在Chapter03/Accessing_Databases.py文件中获得。

小费

你将从 NCBI 访问一个现场应用编程接口 ( API )。请注意，系统的性能在一天中可能会有所不同。此外，在使用它的时候，你应该是一个“好公民”。你会在www.ncbi.nlm.nih.gov/books/NBK25497/#chapter2.找到一些推荐用法指南和要求。值得注意的是，您需要在查询中指定一个电子邮件地址。您应该尽量避免在高峰时段(美国东部时间周一至周五上午 9:00 至下午 5:00)出现大量请求(100 个或更多)，并且每秒发布的查询不要超过三个(Biopython 会为您处理这些问题)。这不仅是良好的公民身份，而且如果你过度使用 NCBI 的服务器，你还有被屏蔽的风险(这是给出真实电子邮件地址的好理由，因为 NCBI 可能会试图联系你)。

怎么做…

现在，让我们看看如何从 NCBI 数据库中搜索和获取数据:

1. 我们将从导入相关模块和配置电子邮件地址开始:

   from Bio import Entrez, SeqIO
   Entrez.email = 'put@your.email.here'
   

我们还将模块导入到流程序列中。不要忘记输入正确的电子邮件地址。

1. 我们将在nucleotide数据库

   handle = Entrez.esearch(db='nucleotide', term='CRT[Gene Name] AND "Plasmodium falciparum"[Organism]')
   rec_list = Entrez.read(handle)
   if int(rec_list['RetMax']) < int(rec_list['Count']):
       handle = Entrez.esearch(db='nucleotide', term='CRT[Gene Name] AND "Plasmodium falciparum"[Organism]', retmax=rec_list['Count'])
       rec_list = Entrez.read(handle)
   

   中Plasmodium falciparum(导致最致命形式疟疾的寄生虫)中寻找氯喹抗性转运蛋白 ( CRT )基因

我们将在nucleotide数据库中搜索我们的基因和生物体(关于搜索字符串的语法，请查看 NCBI 网站)。然后，我们将读取返回的结果。请注意，标准的搜索将记录引用的数量限制为 20 个，因此如果有更多，您可能希望使用增加的最大限制来重复查询。在我们的例子中，我们实际上将使用retmax覆盖默认限制。Entrez系统提供了相当多的复杂方法来检索大量结果(更多信息，请查看 Biopython 或 NCBI·恩特雷兹的文档)。尽管您现在拥有了所有记录的标识符(id)，但是您仍然需要正确地检索记录。

1. 现在，让我们尝试检索所有这些记录。下面的查询将从 GenBank 下载所有匹配的核苷酸序列，在写这本书的时候是 1374。你可能不会一直想这么做:

   id_list = rec_list['IdList']
   hdl = Entrez.efetch(db='nucleotide', id=id_list, rettype='gb')
   

好吧，既然这样，那就去做吧。但是，使用这种技术时要小心，因为您将检索大量完整的记录，其中一些记录内部会有相当大的序列。你冒着下载大量数据的风险(这对你和 NCBI 服务器来说都是一种压力)。

有几种方法可以解决这个问题。一种方法是进行更严格的查询和/或一次只下载几个，当你找到你需要的时候就停下来。精确的策略将取决于你想要达到的目标。无论如何，我们将检索 GenBank 格式的记录列表(包括序列，加上许多有趣的元数据)。

1. 让我们读取并解析结果:

   recs = list(SeqIO.parse(hdl, 'gb'))
   

注意我们已经将一个迭代器(SeqIO.parse的结果)转换成了一个列表。这样做的好处是，我们可以根据需要多次使用结果(例如，多次迭代)，而无需在服务器上重复查询。

如果您计划多次迭代，这将节省时间、带宽和服务器的使用。缺点是它将为所有记录分配内存。这不适用于非常大的数据集；你可能不想像第五章 、中的 那样在全基因组范围内进行这种转换。我们将在本书的最后一部分回到这个话题。如果你在做交互式计算，你可能更喜欢有一个列表(这样你可以多次分析和试验)，但是如果你在开发一个库，迭代器可能是最好的方法。

1. 我们现在只关注一张唱片。只有当您使用完全相同的前面的查询:

   for rec in recs:
       if rec.name == 'KM288867':
           break
   print(rec.name)
   print(rec.description)
   

   时，这才会起作用

rec变量现在有了我们感兴趣的记录。rec.description文件将包含它的可读描述。

1. 现在，让我们提取一些序列特征，这些特征包含序列上的gene产品和exon位置等信息:

   for feature in rec.features:
        if feature.type == 'gene':
            print(feature.qualifiers['gene'])
        elif feature.type == 'exon':
            loc = feature.location
            print(loc.start, loc.end, loc.strand)
        else:
            print('not processed:\n%s' % feature)
   

如果feature.type的值是gene，我们将打印它的名字，它将在qualifiers字典中。我们也将打印所有外显子的位置。外显子，作为具有所有特征的，在这个序列中有位置:一个起点，一个终点，以及它们被读取的链。虽然我们的外显子的所有起始和结束位置都是ExactPosition，但是请注意，Biopython 支持许多其他类型的位置。一类位置是BeforePosition，指定一个定位点在某个序列位置之前。另一类位置是BetweenPosition，给出某个位置开始/结束的间隔。职位类型比较多；这些只是一些例子。

坐标将以这样一种方式指定，即您将能够轻松地从带有范围的 Python 数组中检索序列，因此一般来说，开始将比记录上的值早一个，结束将是相等的。坐标系统的问题将在未来的食谱中重新讨论。

对于其他特征类型，我们简单地打印它们。请注意，当您打印时，Biopython 将提供该功能的人类可读版本。

1. 我们现在将查看记录上的注释，其中大部分是与序列位置不相关的元数据:

   for name, value in rec.annotations.items():
       print('%s=%s' % (name, value))
   

注意，有些值不是字符串；它们可以是数字，甚至是列表(例如，分类法注释就是一个列表)。

1. 最后但同样重要的是，您可以访问一条基本信息——序列:

   print(len(rec.seq))
   

序列对象将是我们下一个配方的主题。

还有更多…

在 NCBI 有更多的数据库。如果您正在处理 NGS 数据，您可能会想要检查序列读取档案 ( SRA )数据库(以前称为短读取档案)。SNP 数据库包含关于 SNP 的信息，而蛋白质数据库包含蛋白质序列，等等。Entrez 中数据库的完整列表链接在中，参见本菜谱的部分。

另一个你可能已经知道的关于 NCBI 的数据库是 PubMed，它包括一个科学和医学引文、摘要、甚至全文的列表。你也可以通过 Biopython 访问它。此外，GenBank 记录通常包含 PubMed 的链接。例如，我们可以对以前的记录执行此操作，如下所示:

from Bio import Medline
refs = rec.annotations['references']
for ref in refs:
    if ref.pubmed_id != '':
        print(ref.pubmed_id)
        handle = Entrez.efetch(db='pubmed', id=[ref.pubmed_id], rettype='medline', retmode='text')
        records = Medline.parse(handle)
        for med_rec in records:
            for k, v in med_rec.items():
                print('%s: %s' % (k, v))

这将获取所有引用注释，检查它们是否有 PubMed ID，然后访问 PubMed 数据库来检索记录，解析它们，然后打印它们。

每个记录的输出是一个 Python 字典。注意，在一个典型的 GenBank 记录中有许多对外部数据库的引用。

当然，还有很多 NCBI 之外的其他生物数据库，比如恩森布尔(【http://www.ensembl.org】)和加州大学圣克鲁斯分校 ( UCSC )基因组生物信息学(genome.ucsc.edu/)。Python 对这些数据库的支持会有很大不同。

如果没有对基本局部比对搜索工具 ( BLAST )的引用，一本关于生物数据库的入门食谱将是不完整的。BLAST 是一种评估序列相似性的算法。NCBI 提供一项服务，允许你将你感兴趣的序列与它自己的数据库进行比较。当然，你可以使用你当地的 BLAST 数据库，而不是使用 NCBI 的服务。Biopython 为此提供了广泛的支持，但由于这太过初级，我将只向您推荐 Biopython 教程。

参见

这些附加信息也很有用:

你可以在biopython.org/DIST/docs/tutorial/Tutorial.xhtml找到更多关于 Biopython 教程的例子。

在 http://www.ncbi.nlm.nih.gov/gquery/可以找到可访问的 NCBI 数据库列表。

一个很棒的问答 ( Q & A )网站是 Biostars(www.biostars.org)，在那里你可以找到关于数据库和序列分析问题的帮助；你可以将它用于本书的所有内容，而不仅仅是这个食谱。

执行基本序列分析

我们现在将做一些 DNA 序列的基本分析。我们将使用 FASTA 文件并进行一些操作，比如反向互补或转录。和前面的菜谱一样，我们将使用你在 第一章**Python 以及周边软件生态中安装的 Biopython。这两个配方为您提供了必要的介绍性构件，我们将使用它们执行所有现代 NGS 分析，然后在本章和第五章、中进行基因组处理。

*## 准备就绪

该配方的代码可在Chapter03/Basic_Sequence_Processing.py中找到。我们将以人类的乳糖酶 ( LCT )基因为例；您可以通过使用Entrez研究界面，利用您从之前的配方中获得的知识来获得此信息:

from Bio import Entrez, SeqIO, SeqRecord
Entrez.email = "your@email.here"
hdl = Entrez.efetch(db='nucleotide', id=['NM_002299'], rettype='gb') # Lactase gene
gb_rec = SeqIO.read(hdl, 'gb')

现在我们有了 GenBank 记录，让我们提取基因序列。记录还不止这些，但让我们先找到基因的精确位置:

for feature in gb_rec.features:
    if feature.type == 'CDS':
        location = feature.location  # Note translation existing
cds = SeqRecord.SeqRecord(gb_rec.seq[location.start:location.end], 'NM_002299', description='LCT CDS only')

我们的示例序列可在 Biopython 序列记录中找到。

怎么做…

让我们看一看下面的步骤:

1. 由于我们感兴趣的序列在 Biopython sequence 对象中是可用的，所以让我们首先将它保存到本地磁盘上的 FASTA 文件中:

   from Bio import SeqIO
   w_hdl = open('example.fasta', 'w')
   SeqIO.write([cds], w_hdl, 'fasta')
   w_hdl.close()
   

SeqIO.write函数接受要写入的序列列表(在我们的例子中，只有一个)。小心这个习语。如果你想写很多序列(你可以用 NGS 轻松地写几百万)，不要使用列表(如前面的代码片段所示)，因为这将分配大量的内存。在每次写操作中，对序列的子集使用迭代器或SeqIO.write函数几次。

1. 在大多数情况下，你实际上会有序列在磁盘上，所以你会有兴趣去读它:

   recs = SeqIO.parse('example.fasta', 'fasta')
   for rec in recs:
       seq = rec.seq
       print(rec.description)
       print(seq[:10])
   

这里，我们关心的是处理单个序列，但是 FASTA 文件可以包含多个记录。Python 习语来执行这个是相当容易的。要读取 FASTA 文件，只需使用标准的迭代技术，如下面的代码片段所示。对于我们的示例，前面的代码将打印以下输出:

NM_002299 LCT CDS only
 ATGGAGCTGT

注意，我们打印了seq[:10]。序列对象可以使用典型的数组切片来获取序列的一部分。

1. 因为我们现在有一个明确的 DNA，我们可以转录它如下:

   rna = seq.transcribe()
   print(rna)
   

2. 最后，我们可以把我们的基因翻译成蛋白质:

   prot = seq.translate()
   print(prot)
   

现在，我们有了基因的氨基酸序列。

还有更多…

关于 Biopython 中的序列管理，还可以说得更多，但这主要是介绍性材料，您可以在 Biopython 教程中找到。我认为让您体验一下序列管理是很重要的，主要是为了完成的目的。为了支持那些可能在生物信息学的其他领域有一些经验但刚刚开始序列分析的人，有几点你应该知道:

当你进行 RNA 翻译以获得蛋白质时，一定要使用正确的遗传密码。即使你在和“普通”生物(比如人类)一起工作，记住线粒体遗传密码是不同的。

Biopython 的Seq对象比这里显示的更加灵活。关于一些好的例子，请参考 Biopython 教程。然而，这个食谱对于我们需要用 FASTQ 文件做的工作来说已经足够了(见下一个食谱)。

为了处理与链相关的问题，正如预期的那样，有一些序列函数，如reverse_complement。

我们开始的 GenBank 记录包含了大量关于该序列的元数据信息，所以一定要探索它。

参见

Biopython 已知的遗传密码是由 NCBI 指定的，可在 http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi 获得。

和前面的食谱一样，Biopython 教程是你的主要停靠站，可以在biopython.org/DIST/docs/tutorial/Tutorial.xhtml找到。

请务必在biopython.org/wiki/SeqIO查看 Biopython SeqIO 页面。

使用现代序列格式

在这里，我们将使用 FASTQ 文件，这是现代测序仪使用的标准格式输出。你将学习如何处理每个碱基的质量分数，并考虑来自不同测序机器和数据库的输出差异。这是第一个将使用来自人类 1000 个基因组项目的真实数据(大数据)的配方。我们将从项目的简要描述开始。

准备就绪

人类 1，000 基因组计划旨在对世界范围内的人类基因变异进行编目，并利用现代测序技术进行 WGS。该项目公开了所有数据,包括测序仪的输出、序列比对和 SNP 调用，以及许多其他工件。“1000 个基因组”这个名称实际上是用词不当，因为它目前包括了 2500 多个样本。这些样本被分成数百个种群，跨越整个星球。我们将主要使用来自四个人群的数据:非洲约鲁巴人(YRI)拥有北欧和西欧血统的犹他州居民(CEU)在东京的日本人 ( JPT )和在北京的汉族人 ( CHB )。我们选择这些特定人群的原因是他们是第一批来自 HapMap 的人群，hap map 是一个有着相似目标的老项目。他们使用基因分型阵列来发现更多关于这个子集的质量。我们将在 第六章**群体遗传学中重温 1000 个基因组和单体型图项目。

小费

下一代数据集通常非常大。因为我们将使用真实的数据，你将下载的一些文件将会很大。虽然我已经尽可能地选择了最小的实例，但是您仍然需要良好的网络连接和相当大的磁盘空间。等待下载可能是这个食谱中最大的障碍，但是数据管理对 NGS 来说是一个严重的问题。在现实生活中，您将需要为数据传输安排时间、分配磁盘空间(这在财务上可能很昂贵)，并考虑备份策略。NGS 最常见的初始错误是认为这些问题微不足道，但事实并非如此。像将一组 BAM 文件复制到网络，甚至复制到您的计算机这样的操作将会变得令人头痛。做好准备。下载大文件后，最起码要检查一下大小是否正确。一些数据库提供消息摘要 5 ( MD5 )校验和。您可以使用 md5sum 等工具将这些校验和与您下载的文件中的校验和进行比较。

下载数据的说明在笔记本的顶部，如Chapter03/Working_with_FASTQ.py的第一个单元格所示。这是一个相当小的文件(27 兆字节 ( 兆字节))，代表一个约鲁巴人女性(NA18489)的部分测序数据。如果你参考 1000 个基因组项目，你会发现绝大多数 FASTQ 文件都要大得多(大了两个数量级)。

FASTQ 序列文件的处理将主要使用 Biopython 来执行。

怎么做…

在我们开始编码之前，让我们看一下 FASTQ 文件，其中有许多记录，如下面的代码片段所示:

@SRR003258.1 30443AAXX:1:1:1053:1999 length=51
 ACCCCCCCCCACCCCCCCCCCCCCCCCCCCCCCCCCCACACACACCAACAC
 +
 =IIIIIIIII5IIIIIII>IIII+GIIIIIIIIIIIIII(IIIII01&III

第 1 行以@开始，后面是序列 ID 和描述字符串。描述字符串会因序列器或数据库来源而异，但通常会服从自动解析。

第二行是测序的 DNA，就像 FASTA 文件一样。第三行是一个+符号，有时第一行是描述行。

第四行包含在第二行读取的每个碱基的质量值。每个字母给一个 Phred 质量分数(en.wikipedia.org/wiki/Phred_quality_score)编码，它给每个读数分配一个错误概率。这种编码在不同平台之间会有所不同。请务必在您的特定平台上检查这一点。

让我们来看看以下步骤:

1. 让我们打开文件:

   import gzip
   from Bio import SeqIO
   recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz'),'rt', encoding='utf-8'), 'fastq')
   rec = next(recs)
   print(rec.id, rec.description, rec.seq)
   print(rec.letter_annotations)
   

我们将打开一个 GNU ZIP ( GZIP )文件，这样我们就可以使用 Python gzip模块了。我们还将指定fastq格式。请注意，这种格式的一些变化会影响 Phred 质量分数的解释。您可能希望指定稍微不同的格式。所有格式参见biopython.org/wiki/SeqIO。

小费

你应该通常以压缩格式存储你的 FASTQ 文件。因为这些是文本文件，所以您不仅获得了大量的磁盘空间，还可能获得一些处理时间。虽然解压缩是一个缓慢的过程，但它仍然比从磁盘读取大得多的(未压缩的)文件要快。

我们将先前配方中的标准字段和质量分数打印到rec.letter_annotations中。只要我们选择了正确的解析器，Biopython 就会将所有 Phred 编码字母转换成对数分数，我们很快就会用到。

现在，不要这样做:

recs = list(recs) # do not do it!

虽然这可能适用于一些 FASTA 文件(以及这个非常小的 FASTQ 文件)，但是如果您执行类似这样的操作，您将分配内存，以便可以在内存中加载完整的文件。对于一个普通的 FASTQ 文件，这是使你的计算机崩溃的最好方法。作为一个规则，总是迭代你的文件。如果您必须对它执行几个操作，您有两个主要的选择。第一种选择是执行一次迭代或者一次执行所有操作。第二种选择是多次打开一个文件并重复迭代。

1. 现在，让我们来看看核苷酸的分布读数:

   from collections import defaultdict
   recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')
   cnt = defaultdict(int)
   for rec in recs:
       for letter in rec.seq:
           cnt[letter] += 1
   tot = sum(cnt.values())
   for letter, cnt in cnt.items():
       print('%s: %.2f %d' % (letter, 100\. * cnt / tot, cnt))
   

我们将重新打开该文件，并使用defaultdict来维护 FASTQ 文件中的核苷酸参考计数。如果您从未使用过这种 Python 标准字典类型，您可能会考虑使用它，因为它消除了初始化字典条目的需要，并为每种类型假定默认值。

注意

N通话有剩余号码。在这些呼叫中，序列发生器报告一个未知碱基。在我们的 FASTQ 文件示例中，我们有一点作弊，因为我们使用了一个过滤文件(调用N的比例将非常低)。预计从序列器中出来的文件中会有更多的N调用。事实上，你甚至可以期待更多关于N叫声的空间分布。

1. 让我们根据它们的读取位置来绘制N的分布:

   import seaborn as sns
   import matplotlib.pyplot as plt
   recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')
   n_cnt = defaultdict(int)
   for rec in recs:
       for i, letter in enumerate(rec.seq):
           pos = i + 1
           if letter == 'N':
               n_cnt[pos] += 1
   seq_len = max(n_cnt.keys())
   positions = range(1, seq_len + 1)
   fig, ax = plt.subplots(figsize=(16,9))
   ax.plot(positions, [n_cnt[x] for x in positions])
   fig.suptitle('Number of N calls as a function of the distance from the start of the sequencer read')
   ax.set_xlim(1, seq_len)
   ax.set_xlabel('Read distance')
   ax.set_ylabel('Number of N Calls')
   

我们导入了seaborn库。尽管在这一点上我们没有明确地使用它，但是这个库的优点是让matplotlib图看起来更好，因为它调整了默认的matplotlib样式。

然后，我们再次打开文件进行解析(记住，您不使用列表，而是再次迭代)。我们遍历文件并获得对N的任何引用的位置。然后，我们将N s 的分布绘制成距离序列起点的函数:

图 3.1–N 次调用的次数与从序列器读取开始的距离的函数关系

您将看到直到位置25，没有错误。这不是典型序列器输出的结果。我们的示例文件已经被过滤，1000 个基因组过滤规则强制要求在位置25之前不能出现N调用。

虽然我们无法研究位置25之前数据集中N s 的行为(您可以随意使用您自己的未过滤 FASTQ 文件来查看N s 如何在读取位置上分布)，但我们可以看到，在位置25之后，分布远非均匀。这里有一个重要的教训，即未调用碱基的数量取决于位置。那么，阅读的质量如何呢？

1. 我们来研究一下 Phred 分数的分布(也就是我们阅读的质量):

   recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')
   cnt_qual = defaultdict(int)
   for rec in recs:
       for i, qual in enumerate(rec.letter_annotations['phred_quality']):
           if i < 25:
               continue
           cnt_qual[qual] += 1
   tot = sum(cnt_qual.values())
   for qual, cnt in cnt_qual.items():
       print('%d: %.2f %d' % (qual, 100\. * cnt / tot, cnt))
   

我们将从(再次)重新打开文件并初始化默认字典开始。然后我们得到phred_quality字母注释，但是我们从一开始就忽略了长达 24 个碱基对 ( bp )的测序位置(由于我们的 FASTQ 文件的过滤，如果您有一个未过滤的文件，您可能想要放弃这个规则)。我们将质量分数添加到默认字典中，最后打印出来。

注意

简单提醒一下， Phred 质量分数是准确呼叫概率的对数表示。这个概率给定为。因此，10 的 Q 代表 90%的呼叫准确率，20 代表 99%的呼叫准确率，30 将是 99.9%。对于我们的文件，最大准确率将是 99.99%(40)。在某些情况下，值 60 是可能的(99.9999%的准确性)。

1. 更有趣的是，我们可以根据他们的阅读位置来绘制质量的分布:

   recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')
   qual_pos = defaultdict(list)
   for rec in recs:
       for i, qual in enumerate(rec.letter_annotations['phred_quality']):
           if i < 25 or qual == 40:
              continue
           pos = i + 1
           qual_pos[pos].append(qual)
   vps = []
   poses = list(qual_pos.keys())
   poses.sort()
   for pos in poses:
       vps.append(qual_pos[pos])
   fig, ax = plt.subplots(figsize=(16,9))
   sns.boxplot(data=vps, ax=ax)
   ax.set_xticklabels([str(x) for x in range(26, max(qual_pos.keys()) + 1)])
   ax.set_xlabel('Read distance')
   ax.set_ylabel('PHRED score')
   fig.suptitle('Distribution of PHRED scores as a function of read distance')
   

在这种情况下，我们将忽略从一开始就被排序为25 bp 的两个位置(同样，如果您有未过滤的测序数据，请删除此规则)和该文件的最高质量分数(40)。但是，在您的情况下，您可以考虑以最大值开始绘图分析。您可能需要检查音序器硬件的最大可能值。一般来说，由于大多数呼叫可以以最高质量执行，如果您试图了解质量问题所在，您可能希望删除它们。

注意我们使用的是seaborn的boxplot功能；我们使用这个只是因为输出看起来比标准的matplotlib的boxplot函数稍微好一点。如果你不想依赖seaborn，就使用股票matplotlib功能。这种情况下，你会调用ax.boxplot(vps)而不是sns.boxplot(data=vps, ax=ax)。

正如所料，分布并不均匀，如下图所示:

图 3.2–ph red 分数的分布，作为从测序仪读数开始的距离的函数

还有更多…

尽管不可能讨论序列器文件输出的所有变化，但成对端读取值得一提，因为它们很常见，需要不同的处理方法。使用成对末端测序，DNA 片段的两端都被测序，中间有一个缺口(称为插入片段)。这种情况下会产生两个文件:X_1.FASTQ和X_2.FASTQ。两个文件将具有相同的顺序和完全相同的序列数。第一个序列将与第一个序列X_2成对出现在X_1中，依此类推。关于编程技术，如果您想保留配对信息，您可以执行如下操作:

f1 = gzip.open('X_1.filt.fastq.gz', 'rt, enconding='utf-8')
f2 = gzip.open('X_2.filt.fastq.gz', 'rt, enconding='utf-8')
recs1 = SeqIO.parse(f1, 'fastq')
recs2 = SeqIO.parse(f2, 'fastq')
cnt = 0
for rec1, rec2 in zip(recs1, recs2):
    cnt +=1
print('Number of pairs: %d' % cnt)

前面的代码按顺序读取所有的线对，并且只计算线对的数量。您可能想做更多的事情，但是这公开了一种基于 Python zip函数的方言，允许您同时遍历两个文件。记得用你的FASTQ前缀替换X。

最后，如果你正在对人类基因组进行测序，你可能想要使用来自完整基因组学的测序数据。在这种情况下，请阅读下一个食谱中的还有更多… 部分，在那里我们将简要讨论完整的基因组数据。

参见

以下是一些包含更多信息的链接:

FASTQ 格式的维基百科页面信息量很大(en.wikipedia.org/wiki/FASTQ_format)。

你可以在 http://www.1000genomes.org/找到关于 1000 个基因组项目的更多信息。

有关 Phred 质量分数的信息可在en.wikipedia.org/wiki/Phred_quality_score找到。

Illumina 在www . illumina . com/science/technology/next-generation-sequencing/paired-end-vs-single-read-sequencing . XHTML提供了一个很好的关于配对末端阅读的介绍页面。

Medvedev 等人关于 nature methods 的论文利用下一代测序发现结构变异的计算方法(www . nature . com/nmeth/journal/V6/n11s/ABS/nmeth . 1374 . XHTML)；请注意，这不是开放访问。

使用校准数据

在你从测序仪收到你的数据后，你通常会使用工具，如布伦斯-惠勒比对器 ( bwa)将你的序列与参考基因组进行比对。大多数用户都有自己物种的参考基因组。你可以在 第五章 、使用基因组中阅读更多关于参考基因组的内容。

比对数据最常见的表示是序列比对图 ( 山姆)格式。由于这些文件大部分都很大，你可能会使用它的压缩版本(BAM)。压缩格式是可索引的，用于极快的随机访问(例如，快速找到染色体某一部分的比对)。注意，BAM 文件需要一个索引，通常由 SAMtools 的tabix实用程序创建。SAMtools 可能是使用最广泛的操作 SAM/BAM 文件的工具。

准备就绪

正如在前面的食谱中所讨论的，我们将使用来自 1000 个基因组项目的数据。我们将对女性NA18489的 20 号染色体使用外显子组比对。这才 312 MB。这个个体的全外显子组比对是 14.2 千兆字节 ( GB )，全基因组比对(在 4x 的低覆盖率下)是 40.1 GB。该数据是成对的末端，读数为 76 bp。这在今天很常见，但是处理起来稍微复杂一些。我们会考虑这一点。如果你的数据没有配对，就把下面的菜谱适当简化一下。

Chapter03/Working_with_BAM.py顶部的单元格会为你下载数据。你想要的文件是NA18490_20_exome.bam和NA18490_20_exome.bam.bai。

我们将使用pysam，一个 SAMtools C API 的 Python 包装器。您可以使用以下命令安装它:

conda install –c bioconda pysam

好的——我们开始吧。

怎么做…

在开始编码之前，注意你可以使用samtools view -h检查 BAM 文件(这是如果你安装了 SAMtools，我们推荐，甚至如果你使用基因组分析工具包 ( GATK )或其他用于变量调用的工具)。我们建议您查看一下头文件和前几条记录。SAM 格式太复杂，无法在此描述。互联网上有大量关于它的信息；尽管如此，有时，这些头文件中隐藏着一些非常有趣的信息。

小费

NGS 中最复杂的操作之一是从原始序列数据生成良好的比对文件。这不仅调用了对齐器，还清理了数据。现在，在高质量 BAM 文件的@PG头文件中，您会发现大多数——如果不是全部——用于生成该 BAM 文件的过程的实际命令行。在我们的示例 BAM 文件中，您将找到运行 bwa、SAMtools、GATK IndelRealigner 和 Picard 应用程序套件来清理数据所需的所有信息。请记住，虽然您可以轻松地生成 BAM 文件，但是它后面的程序在 BAM 输入的正确性方面会非常挑剔。例如，如果您使用 GATK 的变体调用程序来生成基因型调用，则必须对文件进行大范围的清理。因此，其他 BAM 文件的头文件可以为您提供生成自己的文件的最佳方式。最后一个建议是，如果你不处理人类数据，试着为你的物种找到好的 bam，因为给定程序的参数可能略有不同。此外，如果您使用 WGS 数据以外的数据，请检查相似类型的测序数据。

让我们看一看下面的步骤:

1. 让我们检查一下头文件:

   import pysam
   bam = pysam.AlignmentFile('NA18489.chrom20.ILLUMINA.bwa.YRI.exome.20121211.bam', 'rb')
   headers = bam.header
   for record_type, records in headers.items():
       print (record_type)
       for i, record in enumerate(records):
           if type(record) == dict:
               print('\t%d' % (i + 1))
               for field, value in record.items():
                   print('\t\t%s\t%s' % (field, value))
           else:
               print('\t\t%s' % record)
   

头表示为一个字典(其中键是record_type)。由于同一个record_type可能有多个实例，所以字典的值是一个列表(其中每个元素都是一个字典，或者有时是一个包含标签/值对的字符串)。

1. 我们现在将检查单个记录。每条记录的数据量相当复杂。这里，我们将关注成对末端阅读的一些基本领域。查看 SAM 文件规范和pysam API 文档了解更多详情:

   for rec in bam:
       if rec.cigarstring.find('M') > -1 and rec.cigarstring.find('S') > -1 and not rec.is_unmapped and not rec.mate_is_unmapped:
       break
   print(rec.query_name, rec.reference_id, bam.getrname(rec.reference_id), rec.reference_start, rec.reference_end)
   print(rec.cigarstring)
   print(rec.query_alignment_start, rec.query_alignment_end, rec.query_alignment_length)
   print(rec.next_reference_id, rec.next_reference_start,rec.template_length)
   print(rec.is_paired, rec.is_proper_pair, rec.is_unmapped, rec.mapping_quality)
   print(rec.query_qualities)
   print(rec.query_alignment_qualities)
   print(rec.query_sequence)
   

请注意，BAM 文件对象可以在其记录上迭代。我们将遍历它，直到我们找到一个记录，它的简洁的特质间隙对齐报告 ( 雪茄)字符串包含一个火柴和一个软夹子。

雪茄线给出了单个碱基排列的指示。序列中被剪切的部分是对齐器未能对齐的部分(但没有从序列中去除)。我们还需要读数、其配对 ID 和位置(配对的，因为我们有配对的末端读数),它们被映射到参考基因组。

首先，我们打印查询模板名称，后面跟着引用 ID。引用 ID 是一个指针，指向引用查找表中给定引用的序列名称。举个例子就能说明这一点。对于这个 BAM 文件上的所有记录，引用 ID 是19(一个无信息的数字)，但是如果应用bam.getrname(19)，就会得到20，这是染色体的名称。因此，不要混淆参考 ID(在本例中为19)和染色体名称(20)。随后是参考起点和参考终点。pysam是从 0 开始的，不是从 1 开始的，所以当你把坐标转换到其他库的时候要小心。您会注意到本例的起点和终点分别是 59，996 和 60，048，这意味着 52 个碱基的比对。当读取大小为 76 时，为什么只有 52 个碱基(还记得这个 BAM 文件中使用的读取大小吗)?答案可以在雪茄串上找到，在我们的例子中是52M24S，它是 52 个碱基的匹配，然后是 24 个碱基的软剪辑。

然后，我们打印对齐的起点和终点，并计算其长度。顺便说一下，你可以通过观察雪茄绳来计算。它从 0 开始(因为读取的第一部分被映射)并在 52 结束。长度又是 76。

现在，我们查询 mate(只有在有成对末端读取的情况下才会这样做)。我们得到它的引用 ID(如前面的代码片段所示)、它的起始位置以及两对之间的距离度量。只有当配偶双方都映射到同一条染色体上时，这种距离的度量才有意义。

然后，我们绘制该序列的 Phred 分数(参考之前的配方，使用现代序列格式，关于 Phred 分数)，然后仅绘制比对部分的 ph red 分数。最后，我们打印序列(不要忘记这样做！).这是完整的序列，不是剪辑过的序列(当然可以用前面的坐标来剪辑)。

1. 现在，让我们在 BAM 文件中的序列子集中绘制成功作图位置的分布:

   import seaborn as sns
   import matplotlib.pyplot as plt
   counts = [0] * 76
   for n, rec in enumerate(bam.fetch('20', 0, 10000000)):
       for i in range(rec.query_alignment_start, rec.query_alignment_end):
           counts[i] += 1
   freqs = [x / (n + 1.) for x in counts]
   fig, ax = plt.subplots(figsize=(16,9))
   ax.plot(range(1, 77), freqs)
   ax.set_xlabel('Read distance')
   ax.set_ylabel('PHRED score')
   fig.suptitle('Percentage of mapped calls as a function of the position from the start of the sequencer read')
   

我们将开始初始化一个数组来保存整个76位置的计数。注意，然后我们只获取 20 号染色体在 0 号和 10 号位置兆碱基对 ( Mbp )之间的记录。我们将只使用染色体的一小部分。对于这些类型的获取操作，拥有一个索引(由tabix生成)是非常重要的；执行的速度会完全不同。

我们遍历 10 Mbp 边界内的所有记录。对于每个边界，我们得到对齐的开始和结束，并增加对齐的位置之间的可映射性的计数器。最后，我们将其转换为频率，然后绘制出来，如下面的屏幕截图所示:

图 3.3–从序列器读取开始，作为位置函数的映射调用百分比

很明显，可映射性的分布远不是均匀的；极端情况更糟，中间有所下降。

1. 最后，让我们得到 Phred 分数在读数映射部分的分布。正如你可能会怀疑的，这可能不会是统一的:

   from collections import defaultdict
   import numpy as np
   phreds = defaultdict(list)
   for rec in bam.fetch('20', 0, None):
       for i in range(rec.query_alignment_start, rec.query_alignment_end):
           phreds[i].append(rec.query_qualities[i])
   maxs = [max(phreds[i]) for i in range(76)]
   tops = [np.percentile(phreds[i], 95) for i in range(76)]
   medians = [np.percentile(phreds[i], 50) for i in range(76)]
   bottoms = [np.percentile(phreds[i], 5) for i in range(76)]
   medians_fig = [x - y for x, y in zip(medians, bottoms)]
   tops_fig = [x - y for x, y in zip(tops, medians)]
   maxs_fig = [x - y for x, y in zip(maxs, tops)]
   fig, ax = plt.subplots(figsize=(16,9))
   ax.stackplot(range(1, 77), (bottoms, medians_fig,tops_fig))
   ax.plot(range(1, 77), maxs, 'k-')
   ax.set_xlabel('Read distance')
   ax.set_ylabel('PHRED score')
   fig.suptitle('Distribution of PHRED scores as a function of the position in the read')
   

这里，我们再次使用默认字典，允许您使用一些初始化代码。我们现在从头到尾取数据，并在字典中创建一个 Phred 分数列表，该列表的索引是序列读取中的相对位置。

然后我们使用 NumPy 计算第 95、第 50(中间值)和第 5 个百分点，以及每个位置的最大质量分数。对于大多数计算生物学分析，拥有数据的统计汇总视图是很常见的。因此，您可能不仅熟悉百分位数计算，还熟悉计算均值、标准差、最大值和最小值的其他 Pythonic 方法。

最后，我们将绘制每个职位的 Phred 分数分布的堆积图。由于matplotlib期望堆栈的方式，我们必须用stackplot调用之前的值减去下百分位的值。我们可以使用底部百分位数的列表，但是我们必须修正中间值和顶部，如下所示:

图 3.4–ph red 分数的分布，作为读数中位置的函数；底部的蓝色从 0 到第 5 个百分点；绿色到中间值，红色到第 95 百分位，紫色到最大值

还有更多…

虽然我们将在本章的研究基因组可及性和筛选 SNP 数据方法中讨论数据筛选，但我们的目的不是详细解释 SAM 格式或给出数据筛选的详细过程。这个任务需要一本自己的书，但是有了pysam的基础知识，您可以浏览 SAM/BAM 文件。然而，在本章的最后一个方法中，我们将从 BAM 文件中提取全基因组的指标(通过代表 BAM 文件指标的 VCF 文件上的注释),以了解我们数据集的整体质量。

您可能需要处理非常大的数据文件。某些 BAM 处理可能会花费太多时间。减少计算时间的第一种方法是二次采样。例如，如果以 10%进行二次抽样，就会忽略 10 条记录中的 9 条。对于许多任务，例如为 BAM 文件的质量评估所做的一些分析，百分之十(甚至百分之一)的二次抽样就足以获得文件质量的要点。

如果你使用人类数据，你可能会在 Complete Genomics 得到你的数据。在这种情况下，校准文件会有所不同。尽管 Complete Genomics 提供了转换成标准格式的工具，但如果你使用它自己的数据，你可能会得到更好的服务。

参见

可以在以下链接中找到更多信息:

SAM/BAM 格式在samtools.github.io/hts-specs/SAMv1.pdf进行了描述。

你可以在位于genome.sph.umich.edu/wiki/SAM的阿贝卡西集团维基页面上找到关于 SAM 格式的介绍性解释。

如果你真的需要从 BAM 文件中获取复杂的统计数据，Alistair Miles 的pysamstats库是你的停靠港，在 https://github.com/alimanfoo/pysamstats。

为了将你的原始序列数据转换成比对数据，你需要一个比对器；应用最广泛的是 BWA(bio-bwa.sourceforge.net/)。

Picard(肯定是指星际旅行:下一代)是清理 BAM 文件最常用的工具；参考 http://broadinstitute.github.io/picard/的。

用于序列分析的技术论坛被称为SEQanswers(seqanswers.com/)。

我想在这里重复一下关于 Biostars 的建议(在前一个配方中提到过，使用现代序列格式)；它是一个信息宝库，在 http://www.biostars.org/的有一个非常友好的社区。

如果你有完整的基因组学数据，看看 http://www.completegenomics.com/customer-support/faqs/的常见问题 ( 常见问题)。

从 VCF 文件中提取数据

在运行基因型调用者(例如 GATK 或 SAMtools)后，你将有一个 VCF 文件报告关于基因组变异的，比如 SNPs、插入/缺失 ( INDELs )、拷贝数变异 ( CNVs )等等。在这个配方中，我们将通过cyvcf2模块讨论 VCF 加工。

准备就绪

虽然《NGS》讲的都是大数据，但我可以要求你下载多少作为这本书的数据集是有限制的。我相信 2 到 20 GB 的数据对于一个教程来说要求太高了。虽然这个 OOM 中有 1000 个带有真实注释的 VCF 基因组文件，但我们希望在这里使用更少的数据。幸运的是，生物信息学社区已经开发了允许部分下载数据的工具。作为 SAMtools/ htslib包的一部分(【http://www.htslib.org/】)你可以下载tabix和bgzip，它们会负责数据管理。在命令行上，执行以下操作:

tabix -fh ftp://ftp-
trace.ncbi.nih.gov/1000genomes/ftp/release/20130502/supporting/vcf_with_sample_level_annotation/ALL.chr22.phase3_shapeit2_mvncall_integrated_v5_extra_anno.20130502.genotypes.vcf.gz 22:1-17000000 | bgzip -c > genotypes.vcf.gz
tabix -p vcf genotypes.vcf.gz

第一行将部分下载 1000 个基因组项目的 22 号染色体的 VCF 文件(高达 17 Mbp)。然后，bgzip会压缩它。

第二行将创建一个索引，我们将需要它来直接访问基因组的一部分。像往常一样，在一个笔记本中有这样做的代码(Chapter03/Working_with_VCF.py文件)。

您需要安装cyvcf2:

conda install –c bioconda cyvcf2

小费

如果您有冲突解决问题，您可以尝试使用pip来代替。这是你会发现自己用conda做的最后一招，因为它不能解决包依赖，这是经常发生的事情。你可以执行pip install cyvcf2。

怎么做…

看一看下面的步骤:

1. 让我们从检查每条记录可以获得的信息开始:

   from cyvcf2 import VCF
   v = VCF('genotypes.vcf.gz')
   rec = next(v)
   print('Variant Level information')
   info = rec.INFO
   for info in rec.INFO:
       print(info)
   print('Sample Level information')
   for fmt in rec.FORMAT:
       print(fmt)
   

我们首先检查每个记录可用的注释(记住每个记录编码一个变体，如 SNP、CNV、INDELs 等，以及每个样本中该变体的状态)。在变异(记录)水平，我们发现AC—被调用基因型中的ALT等位基因总数、AF—估计的等位基因频率、NS—有数据的样本数、AN—被调用基因型中的等位基因总数、DP—总读取深度。还有其他的，但它们大多是针对 1000 个基因组项目的(在这里，我们将尝试尽可能地通用)。您自己的数据集可能有更多注释(或者没有注释)。

在样本级别，这个文件中只有两个注释:GT—基因型，和DP—每个样本的读取深度。您有每个变量(总)的读取深度和每个样本的读取深度；一定不要混淆两者。

1. 现在我们知道了什么信息是可用的，让我们检查一个单独的 VCF 记录:

   v = VCF('genotypes.vcf.gz')
   samples = v.samples
   print(len(samples))
   variant = next(v)
   print(variant.CHROM, variant.POS, variant.ID, variant.REF, variant.ALT, variant.QUAL, variant.FILTER)
   print(variant.INFO)
   print(variant.FORMAT)
   print(variant.is_snp)
   str_alleles = variant.gt_bases[0]
   alleles = variant.genotypes[0][0:2]
   is_phased = variant.genotypes[0][2]
   print(str_alleles, alleles, is_phased)
   print(variant.format('DP')[0])
   

我们将从检索标准信息开始:染色体、位置、ID、参考碱基(通常只有一个)和备选碱基(可以有多个，但作为第一种过滤方法，通常只接受单个ALT，例如，只接受双等位基因 SNPs)、质量(如您所料，Phred-scaled)和过滤器状态。关于过滤器状态，记住无论 VCF 的文件说什么，你可能仍然想要应用额外的过滤器(在下一个食谱中，研究基因组可及性和过滤 SNP 数据)。

然后我们打印附加的变量级信息(AC、AS、AF、AN、DP等等)，接着是样本格式(在本例中是DP和GT)。最后，我们计算样本的数量，并检查单个样本，以检查它是否是针对此变体调用的。此外，还包括报道的等位基因、杂合性和分期状态(这个数据集恰好是分期的，这并不常见)。

1. 让我们检查一次变异的类型和非双碱基 SNP 的数量:

   from collections import defaultdict
   f = VCF('genotypes.vcf.gz')
   my_type = defaultdict(int)
   num_alts = defaultdict(int)
   for variant in f:
       my_type[variant.var_type, variant.var_subtype] += 1
       if variant.var_type == 'snp':
           num_alts[len(variant.ALT)] += 1
   print(my_type)
   

我们现在将使用现在通用的 Python 默认字典。我们发现这个数据集有 INDELs、CNVs，当然还有 SNPs(大约三分之二是转换，三分之一是颠换)。有剩余数量(79)的三等位基因 SNPs。

还有更多…

本食谱的目的是让你快速掌握cyvcf2模块。在这个阶段，您应该对这个 API 很熟悉了。我们不会在用法细节上花太多时间，因为这将是下一个诀窍的主要目的:使用 VCF 模块来研究你的变体调用的质量。

虽然cyvcf2非常快，但处理基于文本的 VCF 文件仍然需要很多时间。处理这个问题有两个主要策略。一种策略是并行处理，我们将在最后一章 第九章**生物信息管道中讨论。第二个策略是转换成更有效的格式；我们将在第六章 、群体遗传学中提供一个这样的例子。请注意，VCF 开发者正在开发一个二进制变体调用格式 ( BCF )版本来处理这些问题的一部分(www . 1000 genomes . org/wiki/analysis/Variant-Call-Format/BCF-Binary-vcf-version-2)。

参见

一些有用的链接如下:

在 http://samtools.github.io/hts-specs/VCFv4.2.pdf 可以得到 VCF 的规格。

GATK 是使用最广泛的变体调用者之一；详情请查看www.broadinstitute.org/gatk/。

SAMtools 和htslib用于变量调用和 SAM/BAM 管理；详情请看htslib.org。

研究基因组可及性和筛选 SNP 数据

虽然之前的方法侧重于给出 Python 库的概述，以处理比对和变体调用数据，但在本方法中，我们将专注于在头脑中清楚地使用它们。

如果您正在使用 NGS 数据，那么您要分析的最重要的文件很可能是 VCF 文件，它是由基因型调用者如 SAMtools、mpileup或 GATK 生成的。您的 VCF 通话质量可能需要评估和过滤。在这里，我们将建立一个框架来过滤 SNP 数据。我们将为您提供评估数据质量的程序，而不是为您提供过滤规则(一般情况下不可能完成的任务)。有了这个，你可以设计自己的过滤器。

准备就绪

在最好的情况下，您有一个应用了适当过滤器的 VCF 文件。如果是这种情况，您可以继续使用您的文件。注意，所有的 VCF 文件都有一个FILTER列，但是这并不意味着所有合适的过滤器都被应用了。你必须确保你的数据得到了适当的过滤。

在第二种情况下，这是最常见的情况之一，您的文件将有未过滤的数据，但您将有足够的注释，并可以应用硬过滤器(不需要编程过滤)。如果你有一个 GATK 注释文件，请参考gatkforums . broadinstitute . org/discussion/2806/how to-apply-hard-filters-to-a-call-set。

在第三种情况下，您有一个 VCF 文件，其中包含您需要的所有注释，但您可能希望应用更灵活的过滤器(例如，“如果读取深度> 20，如果贴图质量> 30，则接受；否则，如果贴图质量> 40，则接受”)。

在第四种情况下，您的 VCF 文件没有所有必需的注释，您必须重新访问您的 BAM 文件(或者甚至其他信息源)。在这种情况下，最好的解决方案是找到尽可能多的额外信息，并用所需的注释创建一个新的 VCF 文件。一些基因型调用者(如 GATK)允许您指定想要的注释；你可能还想使用额外的程序来提供更多的注释。例如，SNP eff(snpeff.sourceforge.net/)会用它们的效果预测来注释你的 SNP(例如，如果它们在外显子中，它们是编码的还是非编码的？).

不可能提供一个明确的方法，因为它会随着你的测序数据类型、你的研究种类、你对错误的容忍度以及其他变量而变化。我们能做的是提供一组典型的分析，以便进行高质量的过滤。

在这个食谱中，我们不会使用来自人类 1，000 基因组计划的数据。我们想要脏的，未过滤的数据，有许多可以用来过滤它的公共注释。我们将使用冈比亚按蚊 1，000 基因组项目(按蚊是一种蚊子媒介，参与传播导致疟疾的寄生虫)的数据，该项目提供过滤和未过滤的数据。你可以在 http://www.malariagen.net/projects/vector/ag1000g 的找到更多关于这个项目的信息。

我们将得到大约 100 只蚊子的染色体着丝粒的一部分(?? ),随后是该染色体的中间部分(?? )(并标记两者):

tabix -fh ftp://ngs.sanger.ac.uk/production/ag1000g/phase1/preview/ag1000g.AC.phase1.AR1.vcf.gz 3L:1-200000 |bgzip -c > centro.vcf.gz
tabix -fh ftp://ngs.sanger.ac.uk/production/ag1000g/phase1/preview/ag1000g.AC.phase1.AR1.vcf.gz 3L:21000001-21200000 |bgzip -c > standard.vcf.gz
tabix -p vcf centro.vcf.gz
tabix -p vcf standard.vcf.gz

如果链接不起作用，请务必查看github . com/packt publishing/Bioinformatics-with-Python-Cookbook-third-edition/blob/main/datasets . py获取更新。像往常一样，下载这些数据的代码可以在Chapter02/Filtering_SNPs.ipynb笔记本中找到。

最后，对这个食谱提出一个警告:这里的 Python 水平会比平常稍微复杂一些。我们编写的代码越通用，您就越容易在特定情况下重用它。我们将广泛使用函数式编程技术(lambda函数)和partial函数应用。

怎么做…

看看下面的步骤:

1. 让我们从开始，在两个文件中绘制跨基因组的变异分布:

   from collections import defaultdict
   import functools
   import numpy as np
   import seaborn as sns
   import matplotlib.pyplot as plt
   from cyvcf2 import VCF
   def do_window(recs, size, fun):
       start = None
       win_res = []
       for rec in recs:
           if not rec.is_snp or len(rec.ALT) > 1:
               continue
           if start is None:
               start = rec.POS
           my_win = 1 + (rec.POS - start) // size
           while len(win_res) < my_win:
               win_res.append([])
           win_res[my_win - 1].extend(fun(rec))
       return win_res
   wins = {}
   size = 2000
   names = ['centro.vcf.gz', 'standard.vcf.gz']
   for name in names:
    recs = VCF(name)
    wins[name] = do_window(recs, size, lambda x: [1])
   

我们将从执行所需的导入开始(和往常一样，如果您没有使用 IPython 笔记本，记得删除第一行)。在我解释这个函数之前，请注意我们在做什么。

对于这两个文件，我们将计算窗口统计。我们将把包含 200，000 bp 数据的文件分成大小为 2，000 的窗口(100 个窗口)。我们每发现一个双等位基因 SNP，就会在window函数中给这个窗口相关的列表加 1。

window函数将获取一个 VCF 记录(一个非双等位基因透镜(rec.ALT) == 1的rec.is_snp SNP)，确定该记录所属的窗口(通过按大小执行rec.POS的整数除法)，并通过作为fun参数传递给它的函数扩展该窗口的结果列表(在我们的示例中，该参数仅为 1)。

所以，现在我们有一个 100 个元素的列表(每个元素代表 2000 个 bp)。每个元素将是另一个列表，对于找到的每个双等位基因 SNP，该列表将具有 1。

因此，如果在第一个 2,000 bp 中有 200 个 SNP，列表的第一个元素将有 200 个。

1. 我们继续，如下:

   def apply_win_funs(wins, funs):
       fun_results = []
       for win in wins:
           my_funs = {}
           for name, fun in funs.items():
               try:
                   my_funs[name] = fun(win)
               except:
                   my_funs[name] = None
           fun_results.append(my_funs)
       return fun_results
   stats = {}
   fig, ax = plt.subplots(figsize=(16, 9))
   for name, nwins in wins.items():
       stats[name] = apply_win_funs(nwins, {'sum': sum})
       x_lim = [i * size for i in range(len(stats[name]))]
       ax.plot(x_lim, [x['sum'] for x in stats[name]], label=name)
   ax.legend()
   ax.set_xlabel('Genomic location in the downloaded segment')
   ax.set_ylabel('Number of variant sites (bi-allelic SNPs)')
   fig.suptitle('Number of bi-allelic SNPs along the genome', fontsize='xx-large')
   

在这里，我们为我们的 100 个窗口中的每一个执行包含统计信息的绘图。apply_win_funs将为每个窗口计算一组统计数据。在这种情况下，它将对窗口中的所有数字求和。记住，每当我们发现一个 SNP，我们就在窗口列表中加 1。这意味着如果我们有 200 个 SNP，我们将有 200 个；因此，对它们求和将得出 200。

因此，我们能够以一种明显复杂的方式计算每个窗口中 SNPs 的数量。为什么我们要用这种策略做事，很快就会变得很明显。然而，现在，让我们检查两个文件的计算结果,如下面的截图所示:

图 3.5-着丝粒附近(橙色)和染色体中间(蓝色)200 千碱基对(kbp)区域的 2,000 bp 大小的双等位基因 SNP 分布窗口的数量；这两个区域都来自 3L 按蚊 1000 基因组项目中的大约 100 只乌干达蚊子的染色体

小费

注意着丝粒中的 SNPs 数量比染色体中部的少。这是意料之中的，因为在染色体中调用变体比在中间调用更困难。此外，着丝粒中的基因组多样性可能更少。如果你习惯于人类或其他哺乳动物，你会发现变种的密度高得令人讨厌——对你来说就是蚊子！

1. 让我们看一下样本级注释。我们将检查映射质量零(参见www . broadinstitute . org/gatk/guide/tool docs/org _ broadinstitute _ gatk _ tools _ walkers _ annotator _ mappingqualityzerobysample . PHP)，这是一个衡量调用此变体所涉及的序列如何清楚地映射到此位置的指标。注意，在变量级别还有一个MQ0注释:

   mq0_wins = {}
   size = 5000
   def get_sample(rec, annot, my_type):
       return [v for v in rec.format(annot) if v > np.iinfo(my_type).min]
   for vcf_name in vcf_names:
       recs = vcf.Reader(filename=vcf_name)
       mq0_wins[vcf_name] = do_window(recs, size, functools.partial(get_sample, annot='MQ0', my_type=np.int32))
   

通过查看最后一个for开始检查这个；我们将通过从每个记录中读取MQ0注释来执行窗口分析。我们通过调用get_sample函数来实现这一点，该函数将返回我们首选的注释(在本例中为MQ0),该注释已被转换为某种类型(my_type=np.int32)。这里我们使用partial应用函数。Python 允许您指定一个函数的一些参数，并等待稍后指定其他参数。注意，这里最复杂的是函数式编程风格。另外，请注意，这使得计算其他样本级注释变得非常容易。把MQ0换成AB、AD、GQ就行了，以此类推。您会立即得到该注释的计算结果。如果注释不是整数类型，没有问题；改编my_type就好。如果你不习惯这种编程风格，那么它是一种很难的编程风格，但是你很快就会收获它的好处。

1. 现在，让我们打印每个窗口(在本例中，大小为 5000)的中间值和前 75 个百分位数:

   stats = {}
   colors = ['b', 'g']
   i = 0
   fig, ax = plt.subplots(figsize=(16, 9))
   for name, nwins in mq0_wins.items():
       stats[name] = apply_win_funs(nwins, {'median':np.median, '75': functools.partial(np.percentile, q=75)})
       x_lim = [j * size for j in range(len(stats[name]))]
       ax.plot(x_lim, [x['median'] for x in stats[name]], label=name, color=colors[i])
       ax.plot(x_lim, [x['75'] for x in stats[name]], '--', color=colors[i])
       i += 1
   ax.legend()
   ax.set_xlabel('Genomic location in the downloaded segment')
   ax.set_ylabel('MQ0')
   fig.suptitle('Distribution of MQ0 along the genome', fontsize='xx-large')
   

注意，我们现在有两个不同的关于apply_win_funs的统计数据(百分位数和中位数)。同样，我们将函数作为参数传递(np.median和np.percentile)，在np.percentile上完成partial函数应用。结果看起来像这样:

图 3.6–样本 SNP 的 MQ0 的中位数(实线)和第 75 百分位(虚线)，分布在着丝粒附近(蓝色)和染色体中间(绿色)200 kbp 区域的 5,000 bp 大小的窗口上；这两个区域都来自按蚊 1000 基因组项目中大约 100 只乌干达蚊子的 3L 染色体

对于standard.vcf.gz文件，中间值MQ0是0(它被绘制在最底部，几乎看不见)。这很好，因为它表明大多数与变体命名相关的序列都清晰地映射到基因组的这个区域。对于centro.vcf.gz文件，MQ0质量较差。此外，在一些区域，基因分型者根本找不到任何变异(因此图表不完整)。

1. 让我们用样本水平的注释 DP 来比较杂合性。这里，我们将对每个 SNP 的样本读取深度 ( DP )绘制杂合性调用的分数。首先，我们将解释结果，然后解释生成它的代码。

下面的屏幕截图显示了在某一深度杂合的呼叫比例:

图 3.7–实线代表在某一深度计算的杂合位点调用的分数；橙色的是着丝粒区；蓝色是“标准”区域；虚线表示每个深度的样本调用次数；这两个区域都来自按蚊 1000 基因组项目中大约 100 只乌干达蚊子的 3L 染色体

在前面的截图中，有两点需要考虑。在非常低的深度，杂合子的比例是有偏差的——在这种情况下，更低。这是有意义的，因为每个位置的读数不允许您对样品中两个等位基因的存在做出正确的估计。因此，您不应该信任深度非常低的调用。

不出所料，着丝粒内的呼叫次数远低于着丝粒外。着丝粒外 SNP 的分布遵循一种常见的模式，这种模式在许多数据集中是可以预期的。

这里显示了的代码:

def get_sample_relation(recs, f1, f2):
    rel = defaultdict(int)
    for rec in recs:
        if not rec.is_snp:
             continue
        for pos in range(len(rec.genotypes)):
            v1 = f1(rec, pos)
            v2 = f2(rec, pos)
            if v1 is None or v2 == np.iinfo(type(v2)).min:
                continue  # We ignore Nones
            rel[(v1, v2)] += 1
            # careful with the size, floats: round?
        #break
    return rel get_sample_relation(recs, f1, f2):
rels = {}
for vcf_name in vcf_names:
    recs = VCF(filename=vcf_name)
    rels[vcf_name] = get_sample_relation(
        recs,
        lambda rec, pos: 1 if rec.genotypes[pos][0] != rec.genotypes[pos][1] else 0,
        lambda rec, pos: rec.format('DP')[pos][0])

从寻找for循环开始。同样，我们使用函数式编程；get_sample_relation函数将遍历所有 SNP 记录，并应用两个函数参数。第一个参数决定杂合性，而第二个参数获得样本DP(记住还有一个DP的变体)。

现在，由于代码如此复杂，我选择了一个由get_sample_relation返回的简单数据结构:一个字典，其中的关键字是一对结果(在本例中是杂合性和DP)以及共享这两个值的 SNP 的总和。有更优雅的数据结构，有不同的权衡。为此，有 SciPy 稀疏矩阵、pandas 数据框架，或者您可以考虑 PyTables。这里的基本点是有一个足够通用的框架来计算几个样本注释之间的关系。

另外，要注意几个标注的维度空间。例如，如果您的注释是 float 类型的，您可能必须对它进行舍入(否则，您的数据结构的大小可能会变得太大)。

1. 现在，让我们看看绘图代码。让我们分两部分来执行。下面是第一部分:

   def plot_hz_rel(dps, ax, ax2, name, rel):
       frac_hz = []
       cnt_dp = []
       for dp in dps:
           hz = 0.0
           cnt = 0
           for khz, kdp in rel.keys():
               if kdp != dp:
                   continue
               cnt += rel[(khz, dp)]
               if khz == 1:
                   hz += rel[(khz, dp)]
           frac_hz.append(hz / cnt)
           cnt_dp.append(cnt)
       ax.plot(dps, frac_hz, label=name)
       ax2.plot(dps, cnt_dp, '--', label=name)
   

该函数将采用由get_sample_relation生成的数据结构，期望键元组的第一个参数是杂合性状态(0=纯合子，1=杂合子)，第二个参数是DP。这样，它将生成两条线:一条是样本分数(在某一深度是杂合子),另一条是 SNP 计数。

1. 现在，让我们调用这个函数:

   fig, ax = plt.subplots(figsize=(16, 9))
   ax2 = ax.twinx()
   for name, rel in rels.items():
       dps = list(set([x[1] for x in rel.keys()]))
   dps.sort()
   plot_hz_rel(dps, ax, ax2, name, rel)
   ax.set_xlim(0, 75)
   ax.set_ylim(0, 0.2)
   ax2.set_ylabel('Quantity of calls')
   ax.set_ylabel('Fraction of Heterozygote calls')
   ax.set_xlabel('Sample Read Depth (DP)')
   ax.legend()
   fig.suptitle('Number of calls per depth and fraction of calls which are Hz', fontsize='xx-large')
   

这里，我们将使用两个轴。在左手边，我们将有杂合 SNPs 的分数。在右边，我们会有 SNP 的数量。然后我们为两个数据文件调用plot_hz_rel。剩下的就是标准的matplotlib代码了。

1. 最后，让我们将DP变体与分类变体级注释(EFF)进行比较。EFF是由 SnpEff 提供的，它告诉我们(除了许多其他事情之外)SNP 的类型(例如，基因间、内含子、编码同义和编码非同义)。按蚊数据集提供了这一有用的注释。让我们从提取变量级注释和函数式编程风格开始:

   def get_variant_relation(recs, f1, f2):
       rel = defaultdict(int)
       for rec in recs:
           if not rec.is_snp:
               continue
       try:
           v1 = f1(rec)
           v2 = f2(rec)
           if v1 is None or v2 is None:
               continue # We ignore Nones
           rel[(v1, v2)] += 1
       except:
           pass
       return rel
   

这里的编程风格类似于get_sample_relation，但我们不会深究任何样本。现在，我们定义将使用的效果类型，并将其效果转换为整数(因为这将允许我们将其用作索引—例如，矩阵)。现在，考虑编码一个分类变量:

accepted_eff = ['INTERGENIC', 'INTRON', 'NON_SYNONYMOUS_CODING', 'SYNONYMOUS_CODING']
def eff_to_int(rec):
    try:
        annot = rec.INFO['EFF']
        master_type = annot.split('(')[0]
        return accepted_eff.index(master_type)
    except ValueError:
        return len(accepted_eff)

1. 我们现在将遍历文件；现在你应该清楚这种风格了:

   eff_mq0s = {}
   for vcf_name in vcf_names:
       recs = VCF(filename=vcf_name)
       eff_mq0s[vcf_name] = get_variant_relation(recs, lambda r: eff_to_int(r), lambda r: int(r.INFO['DP']))
   

2. 最后，我们使用 SNP 效应绘制了DP的分布:

   fig, ax = plt.subplots(figsize=(16,9))
   vcf_name = 'standard.vcf.gz'
   bp_vals = [[] for x in range(len(accepted_eff) + 1)]
   for k, cnt in eff_mq0s[vcf_name].items():
       my_eff, mq0 = k
       bp_vals[my_eff].extend([mq0] * cnt)
   sns.boxplot(data=bp_vals, sym='', ax=ax)
   ax.set_xticklabels(accepted_eff + ['OTHER'])
   ax.set_ylabel('DP (variant)')
   fig.suptitle('Distribution of variant DP per SNP type', fontsize='xx-large')
   

这里，我们只是为非着丝粒文件打印一个boxplot，如下图所示。结果与预期一致:编码区的 SNP 可能更有深度，因为它们位于比基因间 SNP 更容易调用的更复杂的区域；

图 3.8–不同 SNP 效应的变异阅读深度分布的箱线图

还有更多…

过滤单核苷酸多态性和其他基因组特征的整个问题将需要一本自己的书。这种方法将取决于您拥有的测序数据类型、样本数量和潜在的额外信息(例如，样本中的谱系)。

这份食谱实际上非常复杂，但它的某些部分非常幼稚(在一份简单的食谱中，我能强加给你的复杂性是有限度的)。例如，窗口代码不支持重叠窗口。此外，数据结构过于简单。然而，我希望它们能让你了解处理基因组高通量测序数据的一般策略。更多内容可以在 第四章**高级 NGS 处理中阅读。

参见

更多信息可通过以下链接找到:

有许多过滤规则，但我想提醒您注意合理的良好覆盖(明显高于 10 倍)的必要性。参见 Meynert et al. ，全基因组和外显子组测序中的变异体检测灵敏度和偏倚，在www.biomedcentral.com/1471-2105/15/247/。

bcbio-nextgen是一个基于 Python 的管道，用于高通量测序分析，值得一试(bcbio-nextgen.readthedocs.org)。

用 HTSeq 处理 NGS 数据

HTSeq(HTSeq . readthedocs . io)是一个用于处理 NGS 数据的备选库。HTSeq 提供的大多数功能实际上都可以在本书涵盖的其他库中获得，但是您应该知道它是处理 NGS 数据的另一种方式。HTSeq 支持 FASTA、FASTQ、SAM(通过pysam)、VCF、通用特征格式 ( GFF )和浏览器可扩展数据 ( 床)文件格式。它还包括一组用于处理(绘制的)基因组数据的抽象，包括基因组位置和间隔或比对等概念。对这个库的特性进行全面的检查超出了我们的范围，所以我们将集中讨论一小部分特性。我们将借此机会介绍 BED 文件格式。

BED 格式允许指定注释轨迹的特征。它有许多用途，但通常将 BED 文件加载到基因组浏览器中以可视化特征。每行至少包括位置信息(染色体、起点和终点)以及可选字段，如名称或链。关于这种形式的全部细节可以在genome.ucsc.edu/FAQ/FAQformat.xhtml#format1找到。

准备就绪

我们简单的示例将使用来自人类基因组中 LCT 基因所在区域的数据。LCT 基因编码乳糖酶，一种参与乳糖消化的酶。

我们将从 Ensembl 获取这些信息。去 http://uswest.ensembl.org/Homo_sapiens/Gene/Summary?的db =核心；g = ensg 000000115850并选择导出数据。输出格式应为床格式。基因信息要选(如果愿意可以多选)。为了方便起见，在Chapter03目录中有一个名为LCT.bed的下载文件。

这个代码的笔记本叫做Chapter03/Processing_BED_with_HTSeq.py。

在我们开始之前看一下文件。此处提供了该文件的几行示例:

track name=gene description="Gene information"
 2       135836529       135837180       ENSE00002202258 0       -
 2       135833110       135833190       ENSE00001660765 0       -
 2       135789570       135789798       NM_002299.2.16  0       -
 2       135787844       135788544       NM_002299.2.17  0       -
 2       135836529       135837169       CCDS2178.117    0       -
 2       135833110       135833190       CCDS2178.116    0       -

第四个列是特性名称。这将因文件而异，你将不得不每次都检查它。然而，在我们的例子中，似乎很明显我们有 Ensembl 外显子(ENSE…)，GenBank 记录(NM _…)，以及来自共有编码序列 ( CCDS )数据库(www.ncbi.nlm.nih.gov/CCDS/CcdsBrowse.cgi)的编码区信息(CCDS)。

您需要安装 HTSeq:

conda install –c bioconda htseq

现在，我们可以开始了。

怎么做…

看看下面的步骤:

1. 我们将首先为我们的文件设置一个阅读器。记住这个文件已经提供给你了，应该在你当前的工作目录:

   from collections import defaultdict
   import re
   import HTSeq
   lct_bed = HTSeq.BED_Reader('LCT.bed')
   

2. 我们现在将通过名称提取所有类型的特征:

   feature_types = defaultdict(int)
   for rec in lct_bed:
       last_rec = rec
       feature_types[re.search('([A-Z]+)', rec.name).group(0)] += 1
   print(feature_types)
   

记住这段代码是特定于我们的例子的。你必须使它适应你的情况。

小费

你会发现前面的代码使用了一个正则表达式 ( regex )。小心使用正则表达式，因为它们会生成难以维护的只读代码。你可能有更好的选择。无论如何，正则表达式是存在的，你会不时地发现它们。

我们案例的输出如下所示:

defaultdict(<class 'int'>, {'ENSE': 27, 'NM': 17, 'CCDS': 17})

1. 我们存储了最后一条记录，以便我们可以检查它:

   print(last_rec)
   print(last_rec.name)
   print(type(last_rec))
   interval = last_rec.iv
   print(interval)
   print(type(interval))
   

有许多字段可用，最著名的是name和interval。对于前面的代码，输出如下所示:

<GenomicFeature: BED line 'CCDS2178.11' at 2: 135788543 -> 135788322 (strand '-')>
 CCDS2178.11
 <class 'HTSeq.GenomicFeature'>
 2:[135788323,135788544)/-
 <class 'HTSeq._HTSeq.GenomicInterval'>

1. 让我们更深入地挖掘区间:

   print(interval.chrom, interval.start, interval.end)
   print(interval.strand)
   print(interval.length)
   print(interval.start_d)
   print(interval.start_as_pos)
   print(type(interval.start_as_pos))
   

输出如下所示:

2 135788323 135788544
 -
 221
 135788543
 2:135788323/-
 <class 'HTSeq._HTSeq.GenomicPosition'>

注意基因组的位置(染色体，开始和结束)。最复杂的问题是如何处理这个问题。如果特征编码在负链中，你必须小心处理。HTSeq 提供了start_d和end_d字段来帮助您完成这个任务(也就是说，如果链是负的，那么它们的起点和终点将会相反)。

最后，让我们从我们的编码区域(CCDS 记录)中提取一些统计数据。我们将使用 CCDS，因为它可能比这里的策划数据库更好:

exon_start = None
exon_end = None
sizes = []
for rec in lct_bed:
    if not rec.name.startswith('CCDS'):
        continue
    interval = rec.iv
    exon_start = min(interval.start, exon_start or interval.start)
    exon_end = max(interval.length, exon_end or interval.end)
    sizes.append(interval.length)
sizes.sort()
print("Num exons: %d / Begin: %d / End %d" % (len(sizes), exon_start, exon_end))
print("Smaller exon: %d / Larger exon: %d / Mean size: %.1f" % (sizes[0], sizes[-1], sum(sizes)/len(sizes)))

输出应该是不言自明的:

Num exons: 17 / Begin: 135788323 / End 135837169
 Smaller exon: 79 / Larger exon: 1551 / Mean size: 340.2

还有更多…

床的形式可以比这更复杂一点。此外，前面的代码基于关于我们文件内容的非常具体的前提。然而，这个例子应该足以让你开始。即使在最坏的情况下，床的格式也不是很复杂。

HTSeq 有比这多得多的功能，但是这个配方主要是作为整个包的起点。HTSeq 具有一些功能，可以作为我们到目前为止讨论过的大多数食谱的替代方案。*

作者：绝不原创的飞龙