Python-生物信息学秘籍第三版-一-Python 生物信息学秘籍第三版（一）第一章：《Python 生物信息学实用

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

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

译者：飞龙

协议：CC BY-NC-SA 4.0

第一章：《Python 生物信息学实用手册》第三版

使用现代 Python 库和应用程序解决现实世界的计算生物学问题

Tiago Antao

伯明翰—孟买

《Python 生物信息学实用手册》第三版

本书在编写过程中已尽最大努力确保所呈现信息的准确性。然而，本书所包含的信息是没有任何明示或暗示的担保的。无论是作者、Packt Publishing 还是其代理商与分销商，都不对因本书直接或间接造成的或声称已造成的任何损害负责。

Packt Publishing 力求通过适当使用大写字母提供本书中提到的所有公司和产品的商标信息，但无法保证这些信息的准确性。

出版产品经理：Devika Battike

高级编辑：David Sugarman

内容开发编辑：Joseph Sunil

技术编辑：Rahul Limbachiya

文案编辑：Safis Editing

项目协调员：Farheen Fathima

校对员：Safis Editing

索引员：Pratik Shirodkar

制作设计师：Shankar Kalbhor

市场协调员：Priyanka Mhatre

初版：2015 年 6 月

第二版：2018 年 11 月

第三版：2022 年 9 月

出版参考：1090922

由 Packt Publishing Ltd. 出版

Livery Place

35 Livery Street

伯明翰

B3 2PB，英国

ISBN 978-1-80323-642-1

www.packt.com

贡献者

关于作者

Tiago Antao 是一位生物信息学家，目前在基因组学领域工作。曾是计算机科学家，Tiago 拥有葡萄牙波尔图大学科学学院的生物信息学硕士学位，并在英国利物浦热带医学学院获得了关于抗药性疟疾传播的博士学位。博士后阶段，Tiago 在英国剑桥大学从事人类数据集研究，并在英国牛津大学处理蚊子全基因组测序数据，随后协助建立了美国蒙大拿大学的生物信息学基础设施。现在，他在马萨诸塞州波士顿的生物技术领域担任数据工程师。他是 Biopython 这一重要生物信息学 Python 包的共同作者之一。

关于审阅者

Urminder Singh 是一名生物信息学家、计算机科学家以及多个开源生物信息学工具的开发者。他的教育背景涵盖了物理学、计算机科学和计算生物学学位，包括在美国爱荷华州立大学获得的生物信息学博士学位。

他的研究兴趣广泛，包括新型基因进化、精准医学、社会基因组学、医学中的机器学习以及开发大规模异构数据的工具和算法。你可以在网上访问他：urmi-21.github.io。

Tiffany Ho 在 Embark Veterinary 担任生物信息学助理。她拥有加利福尼亚大学戴维斯分校的遗传学和基因组学学士学位，并持有康奈尔大学的植物育种与遗传学硕士学位。

前言

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

本书的更新版《Python 生物信息学实用手册》首先概述了 Python 生态系统中各种工具和库，它们将帮助你转换、分析和可视化生物数据集。在后续章节中，你将学习下一代测序、单细胞分析、基因组学、宏基因组学、群体遗传学、系统发育学和蛋白质组学等关键技术，并通过实际案例深入了解。你将学习如何使用重要的管道系统，如 Galaxy 服务器和 Snakemake，以及如何理解 Python 中用于函数式和异步编程的各种模块。本书还将帮助你探索一些课题，比如在高性能计算框架（包括 Dask 和 Spark）下使用统计方法进行 SNP 发现，以及将机器学习算法应用于生物信息学领域。

本书结束时，你将掌握最新的编程技术和框架，能够处理各种规模的生物信息学数据。

本书适用对象

本书适用于生物信息学分析师、数据科学家、计算生物学家、研究人员以及希望解决中级到高级生物学和生物信息学问题的 Python 开发人员。读者需要具备 Python 编程语言的工作知识，具备基础生物学知识会有所帮助。

本书内容

第一章*，《Python 与相关软件生态》*，告诉你如何使用 Python 搭建现代生物信息学环境。本章讨论了如何使用 Docker 部署软件，如何与 R 交互，以及如何使用 Jupyter Notebooks 进行互动。

第二章*，《了解 NumPy、pandas、Arrow 和 Matplotlib》*，介绍了数据科学中的基础 Python 库：用于数组和矩阵处理的 NumPy；用于表格数据处理的 Pandas；优化 Pandas 处理的 Arrow，以及用于绘图的 Matplotlib。

第三章*，《下一代测序技术（NGS）》*，提供了解决下一代测序数据的具体方案。本章教你如何处理大型的 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、Pandas、Matplolib
Biopython
Dask、zarr、scikit-learn

如果您正在使用本书的数字版，我们建议您自己输入代码或通过 GitHub 仓库访问代码（下一个章节中会提供链接）。这样可以帮助您避免与复制和粘贴代码相关的潜在错误。

下载示例代码文件

您可以从 GitHub 下载本书的示例代码文件，网址为github.com/PacktPublishing/Bioinformatics-with-Python-Cookbook-third-edition。如果代码有任何更新，它会在现有的 GitHub 仓库中更新。

我们还提供了来自丰富书籍和视频目录的其他代码包，您可以在github.com/PacktPublishing/找到。请查看！

下载彩色图像

我们还提供了一份 PDF 文件，其中包含本书中使用的截图/图表的彩色图像。您可以在此处下载：packt.link/3KQQO。

使用的约定

本书中使用了许多文本约定。

文本中的代码：表示文本中的代码词汇、数据库表名、文件夹名、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 用户名。以下是一个示例：“call_genotype 的形状是 56,241x1,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

粗体：表示新术语、重要词汇或屏幕上显示的词汇。例如，菜单或对话框中的文字通常以这种方式出现在文本中。以下是一个示例：“对于Chunk列，请参见 第十一章 —— 但目前您可以暂时忽略它。”

提示或重要说明

以这种方式显示。

部分

本书中，您会看到一些经常出现的标题（准备工作、如何操作...、工作原理...、更多内容... 和另见）。

为了清晰地指导如何完成食谱，请按照以下方式使用这些部分：

准备工作

本节告诉您在食谱中可以期待什么，并描述如何设置所需的软件或进行任何前期设置。

如何操作…

本节包含了遵循食谱所需的步骤。

工作原理…

本节通常包含对上一节所做内容的详细解释。

另见

本节提供了有助于理解食谱的其他有用信息链接。

联系我们

我们始终欢迎读者的反馈。

一般反馈：如果您对本书的任何部分有疑问，请在邮件主题中提及书名，并通过customercare@packtpub.com与我们联系。

勘误：尽管我们已尽力确保内容的准确性，但难免会有错误。如果你在本书中发现了错误，我们将非常感激你向我们报告。请访问 www.packtpub.com/support/err…，选择你的书籍，点击“勘误提交表单”链接，并填写详细信息。

盗版：如果你在互联网上发现任何我们作品的非法复制品，我们将感激不尽，如果你能提供该位置地址或网站名称。请通过 copyright@packt.com 联系我们，并附上材料的链接。

如果你有兴趣成为一名作者：如果你在某个领域有专业知识，并且有兴趣撰写或参与书籍的编写，请访问 authors.packtpub.com。

请留下评论。在阅读并使用本书后，为什么不在你购买书籍的网站上留下评论呢？潜在读者可以通过你的公正意见做出购买决定，我们在 Packt 可以了解你对我们产品的看法，而我们的作者也可以看到你对他们书籍的反馈。谢谢！

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

分享你的想法

一旦你阅读了*《Python 生物信息学实用指南》*，我们很希望听到你的想法！请点击这里直接前往本书的亚马逊评论页面，分享你的反馈。

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

第二章：Python 与周边软件生态

我们将从安装本书大部分内容所需的基本软件开始。包括Python发行版，一些基本的 Python 库，以及外部生物信息学软件。在这里，我们还将关注 Python 之外的世界。在生物信息学和大数据领域，R 也是一个重要的角色；因此，你将学习如何通过 rpy2 与其进行交互，rpy2 是一个 Python/R 桥接工具。此外，我们将探索 IPython 框架（通过 Jupyter Lab）所能带来的优势，以便高效地与 R 进行接口。由于 Git 和 GitHub 在源代码管理中的广泛应用，我们将确保我们的设置能够与它们兼容。本章将为本书剩余部分的所有计算生物学工作奠定基础。

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

本章将涵盖以下内容：

使用 Anaconda 安装所需的软件
使用 Docker 安装所需的软件
通过 rpy2 与 R 交互
在 Jupyter 中使用 R 魔法

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

在我们开始之前，我们需要安装一些基本的前置软件。接下来的章节将引导你完成所需软件及其安装步骤。每个章节和部分可能会有额外的要求，我们会在本书后续部分中明确说明。另一种启动方式是使用 Docker 配方，之后所有工作将通过 Docker 容器自动处理。

如果你已经在使用其他 Python 发行版，强烈建议考虑使用 Anaconda，因为它已经成为数据科学和生物信息学的事实标准。此外，它是可以让你从 Bioconda 安装软件的发行版（bioconda.github.io/）。

准备工作

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

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

如果你使用的是 Windows，考虑升级到 Linux 来进行生物信息学工作，因为大多数现代生物信息学软件无法在 Windows 上运行。请注意，除非你计划使用计算机集群，否则 macOS 对几乎所有分析都是可以的，计算机集群通常基于 Linux。

如果你使用的是 Windows 或 macOS，并且没有轻松访问 Linux 的方式，别担心。现代虚拟化软件（如 VirtualBox 和 Docker）将会帮你解决问题，它们允许你在操作系统上安装一个虚拟的 Linux。如果你正在使用 Windows 并决定使用原生方式而不是 Anaconda，要小心选择库；如果你为所有内容（包括 Python 本身）安装 32 位版本，可能会更安全。

注意

如果你使用的是 Windows，许多工具将无法使用。

生物信息学和数据科学正在以惊人的速度发展，这不仅仅是炒作，而是现实。在安装软件库时，选择版本可能会很棘手。根据你所使用的代码，某些旧版本可能无法工作，或者可能甚至无法与新版本兼容。希望你使用的任何代码会指出正确的依赖项——尽管这并不保证。在本书中，我们将固定所有软件包的精确版本，并确保代码可以与这些版本一起工作。代码可能需要根据其他软件包版本进行调整，这是非常自然的。

本书开发的软件可以在 github.com/PacktPublishing/Bioinformatics-with-Python-Cookbook-third-edition 找到。要访问它，你需要安装 Git。习惯使用 Git 可能是个好主意，因为许多科学计算软件都是在它的基础上开发的。

在正确安装 Python 环境之前，你需要先安装所有与 Python 交互的外部非 Python 软件。这个列表会根据章节有所不同，所有章节特定的软件包将在相应的章节中解释。幸运的是，自本书的前几个版本以来，大多数生物信息学软件已经可以通过 Bioconda 项目安装，因此安装通常很容易。

你需要安装一些开发编译器和库，所有这些都是免费的。在 Ubuntu 上，建议安装 build-essential 包（apt-get install build-essential），在 macOS 上，建议安装 Xcode（developer.apple.com/xcode/）。

在下表中，你将找到开发生物信息学所需的最重要的 Python 软件列表：

名称	应用	网址	目的
Project Jupyter	所有章节	jupyter.org/	交互式计算
pandas	所有章节	pandas.pydata.org/	数据处理
NumPy	所有章节	www.numpy.org/	数组/矩阵处理
SciPy	所有章节	www.scipy.org/	科学计算
Biopython	所有章节	biopython.org/	生物信息学库
seaborn	所有章节	seaborn.pydata.org/	统计图表库
R	生物信息学与统计	www.r-project.org/	统计计算语言
rpy2	R 连接	rpy2.readthedocs.io	R 接口
PyVCF	NGS	pyvcf.readthedocs.io	VCF 处理
Pysam	NGS	github.com/pysam-devel…	SAM/BAM 处理
HTSeq	NGS/基因组	htseq.readthedocs.io	NGS 处理
DendroPY	系统发育学	dendropy.org/	系统发育学
PyMol	蛋白质组学	pymol.org	分子可视化
scikit-learn	机器学习	scikit-learn.org	机器学习库
Cython	大数据	cython.org/	高性能
Numba	大数据	numba.pydata.org/	高性能
Dask	大数据	dask.pydata.org	并行处理

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

我们将使用 pandas 来处理大多数表格数据。另一种选择是仅使用标准 Python。pandas 在数据科学中已经变得非常普及，因此如果数据量适合内存，使用它处理所有表格数据可能是明智之举。

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

怎么做...

要开始，请查看以下步骤：

首先从 www.anaconda.com/products/in… 下载 Anaconda 发行版。我们将使用版本 21.05，尽管您可能可以使用最新版本。您可以接受安装的所有默认设置，但您可能希望确保conda二进制文件在您的路径中（不要忘记打开一个新窗口，以便路径可以更新）。如果您有另一个 Python 发行版，请注意您的PYTHONPATH和现有的 Python 库。最好卸载所有其他 Python 版本和已安装的 Python 库。尽可能地，删除所有其他 Python 版本和已安装的 Python 库。
让我们继续使用库。我们现在将使用以下命令创建一个名为bioinformatics_base的新conda环境，并安装biopython=1.70：
```
conda create -n bioinformatics_base python=3.10
```
让我们按以下步骤激活环境：
```
conda activate bioinformatics_base
```

让我们将bioconda和conda-forge通道添加到我们的源列表：

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

还要安装基本软件包：

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

现在，让我们保存我们的环境，以便以后在其他机器上创建新的环境，或者如果需要清理基础环境：
```
conda list –explicit > bioinformatics_base.txt
```

我们甚至可以从conda安装 R：

conda install rpy2 r-essentials r-gridextra

请注意，r-essentials会安装许多 R 包，包括后面要用到的 ggplot2。另外，我们还会安装r-gridextra，因为我们在 Notebook 中会用到它。

还有更多...

如果您不喜欢使用 Anaconda，您可以选择使用任何发行版通过pip安装许多 Python 库。您可能需要许多编译器和构建工具，不仅包括 C 编译器，还包括 C++和 Fortran。

我们将不再使用前面步骤中创建的环境。相反，我们将把它用作克隆工作环境的基础。这是因为使用 Python 进行环境管理，即使借助conda包系统的帮助，仍可能非常痛苦。因此，我们将创建一个干净的环境，以免造成损害，并可以从中派生出开发环境，如果开发环境变得难以管理。

例如，假设您想创建一个带有scikit-learn的机器学习环境。您可以执行以下操作：

使用以下命令创建原始环境的克隆：

conda create -n scikit-learn --clone bioinformatics_base

添加scikit-learn：

conda activate scikit-learn
conda install scikit-learn

在 JupyterLab 中，我们应该通过 notebook 打开我们的 jupytext 文件，而不是通过文本编辑器。由于 jupytext 文件与 Python 文件具有相同的扩展名——这是一个特性，而非 bug——JupyterLab 默认会使用普通的文本编辑器。当我们打开 jupytext 文件时，我们需要覆盖默认设置。右键点击并选择 Notebook，如下图所示：

图 1.2 – 在 Notebook 中打开 jupytext 文件

我们的 jupytext 文件将不会保存图形输出，这对于本书来说已经足够。如果你想要一个带有图片的版本，这是可以通过配对的 notebooks 实现的。更多详情，请查看 Jupytext 页面 (github.com/mwouts/jupytext)。

警告

由于我们的代码是为了在 Jupyter 中运行的，本书中很多地方，我不会使用 print 输出内容，因为单元格的最后一行会被自动渲染。如果你没有使用 notebook，请记得使用 print。

使用 Docker 安装所需的软件

Docker 是实现操作系统级虚拟化的最广泛使用的框架。这项技术使你能够拥有一个独立的容器：一个比虚拟机更轻量的层，但仍然允许你将软件进行隔离。这基本上隔离了所有进程，使得每个容器看起来像一个虚拟机。

Docker 在开发领域的两端表现得非常好：它是设置本书内容以用于学习目的的便捷方式，并且可以成为你在复杂环境中部署应用程序的首选平台。这个方案是前一个方案的替代。

然而，对于长期的开发环境，类似前一个方案的方法可能是你的最佳选择，尽管它可能需要更繁琐的初始设置。

准备工作

如果你使用的是 Linux，首先需要做的是安装 Docker。最安全的解决方案是从 www.docker.com/ 获取最新版本。虽然你的 Linux 发行版可能有 Docker 包，但它可能过时且有 bug。

如果你使用的是 Windows 或 macOS，不要灰心；查看 Docker 网站。有多种可供选择的方案，但没有明确的公式，因为 Docker 在这些平台上发展得很快。你需要一台相对较新的计算机来运行我们的 64 位虚拟机。如果遇到任何问题，重启计算机并确保 BIOS、VT-X 或 AMD-V 已启用。至少你需要 6 GB 内存，最好更多。

注意

这将需要从互联网进行非常大的下载，因此请确保你有足够的带宽。另外，准备好等待很长时间。

如何操作...

要开始，请按照以下步骤操作：

在 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 组中。

现在你准备好运行容器，具体如下：

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

将 YOUR_DIRECTORY 替换为你操作系统中的一个目录。该目录将在主机操作系统和 Docker 容器之间共享。YOUR_DIRECTORY 在容器内会显示为 /data，反之亦然。

-p 9875:9875 将暴露容器的 TCP 端口 9875 到主机计算机的端口 9875。

特别是在 Windows 上（也许在 macOS 上），请确保你的目录在 Docker shell 环境中确实可见。如果不可见，请查看官方 Docker 文档，了解如何暴露目录。

现在你准备好使用系统了。将浏览器指向 http://localhost:9875，你应该会看到 Jupyter 环境。

如果在 Windows 上无法运行，请查看官方 Docker 文档 (docs.docker.com/)，了解如何暴露端口。

另请参见

以下内容也值得了解：

Docker 是最广泛使用的容器化软件，最近在使用量上增长巨大。你可以在 www.docker.com/ 阅读更多关于它的信息。
作为 Docker 的一个面向安全的替代方案，rkt 可以在 coreos.com/rkt/ 找到。
如果你无法使用 Docker，例如，如果你没有必要的权限（这通常发生在大多数计算集群中），可以查看 Singularity，网址为 www.sylabs.io/singularity/。

通过 rpy2 与 R 进行交互

如果你需要某些功能，但在 Python 库中找不到，首先应该检查该功能是否已经在 R 中实现。对于统计方法，R 仍然是最完整的框架；此外，某些生物信息学功能仅在 R 中可用，并且可能作为 Bioconductor 项目的一部分提供。

ggplot2，我们将从人类 1,000 基因组计划 (www.1000genomes.org/) 下载它的元数据。这本书不是关于 R 的，但我们希望提供有趣且实用的示例。

准备工作

你需要从 1,000 基因组序列索引中获取元数据文件。请查看 github.com/PacktPublishing/Bioinformatics-with-Python-Cookbook-third-edition/blob/main/Datasets.py，并下载 sequence.index 文件。如果你使用 Jupyter Notebook，打开 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

现在我们可以开始了。

如何操作...

开始之前，请按照以下步骤操作：

让我们从导入一些库开始：

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 中的 DataFrame 与pandas非常匹配。

我们将使用 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)

导入后，我们做的第一件事是访问 R 中的read.delim函数，它允许你读取文件。R 语言规范允许在对象名称中使用点。因此，我们必须将函数名转换为read_delim。然后，我们调用正确的函数名；请注意以下几个显著的特性。首先，大多数原子对象（如字符串）可以直接传递而无需转换。其次，参数名也可以无缝转换（除了点问题）。最后，对象可在 Python 命名空间中使用（然而，对象实际上不能在 R 命名空间中使用；我们稍后会详细讨论）。

作为参考，我附上了相应的 R 代码。我希望你能看出它的转换是非常简单的。seq_data对象是一个 DataFrame。如果你了解基本的 R 或者pandas，你可能对这种数据结构有所了解。如果不了解，那么它基本上是一个表格，即一系列行，其中每一列都具有相同的数据类型。

让我们对这个 DataFrame 进行基本的检查，如下所示：

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

再次注意代码的相似性。

你甚至可以使用以下代码混合不同的风格：
```
my_cols = robjects.r.ncol(seq_data)
print(my_cols)
```

你可以直接调用 R 函数；在这种情况下，如果函数名称中没有点，我们将调用ncol；但是要小心。这样做会输出结果，而不是 26（列数），而是[26]，这是一个包含26元素的向量。这是因为默认情况下，R 中的大多数操作返回向量。如果你想要列数，必须执行my_cols[0]。另外，谈到陷阱，注意 R 数组的索引是从 1 开始的，而 Python 是从 0 开始的。

现在，我们需要进行一些数据清理。例如，有些列应该解释为数字，但它们却被当作字符串读取：

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函数会将一列数据转换为整数。第一次打印会显示字符串（即被"包围的值），而第二次打印会显示数字。

我们需要对这张表进行更多处理；有关详细信息可以在笔记本中找到。在这里，我们将完成将 DataFrame 获取到 R 中的操作（请记住，虽然它是一个 R 对象，但实际上在 Python 命名空间中可见）：
```
robjects.r.assign('seq.data', seq_data)
```

这将在 R 命名空间中创建一个名为seq.data的变量，其内容来自 Python 命名空间中的 DataFrame。请注意，在此操作之后，两个对象将是独立的（如果更改其中一个对象，将不会反映到另一个对象中）。

注意

尽管你可以在 Python 上执行绘图，R 默认内置了绘图功能（这里我们将忽略）。它还有一个叫做ggplot2的库，实现了图形语法（一种声明性语言，用于指定统计图表）。

关于我们基于人类 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 打印设备来完成这一操作。

现在，我们可以在 Jupyter Notebook 中打印图像：
```
Image(filename='out.png')
```

生成以下图表：

图 1.3 – 使用 ggplot2 生成的中心名称直方图，负责从 1,000 个基因组计划中测序人类基因组数据的通道

作为最后一个示例，我们现在将对 Yoruban（YRI）和来自北欧和西欧的犹他州居民（CEU）的所有测序通道的读取和碱基计数进行散点图绘制，使用人类 1,000 个基因组计划（该项目的数据总结，我们将充分使用，可以在第三章的使用现代序列格式食谱中看到）。此外，我们对不同类型的测序之间的差异感兴趣（例如，外显子覆盖率、高覆盖率和低覆盖率）。首先，我们生成一个仅包含YRI和CEU通道的 DataFrame，并限制最大碱基和读取计数：
```
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')
```

现在我们已经准备好绘图了：

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

希望这个例子（请参考下面的截图）能清楚地展示图形语法方法的强大。我们将首先声明 DataFrame 和正在使用的图表类型（即通过 geom_point 实现的散点图）。

注意到很容易表达每个点的形状取决于 POPULATION 变量，而颜色取决于 ANALYSIS_GROUP 变量：

图 1.4 – 由 ggplot2 生成的散点图，显示所有测序通道的基础和读取计数；每个点的颜色和形状反映了类别数据（种群和所测序的数据类型）

由于 R DataFrame 与 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 DataFrame（请注意，我们正在转换 R 命名空间中的 yri_ceu，而不是 Python 命名空间中的）。我们删除了 pandas DataFrame 中指示配对 FASTQ 文件名称的列，并将其复制回 R 命名空间。如果你打印新 R DataFrame 的列名，你会发现 PAIRED_FASTQ 列丢失了。

还有更多内容...

值得重复的是，Python 软件生态的进展速度非常快。这意味着如果今天某个功能不可用，它可能在不久的将来会发布。因此，如果你正在开发一个新项目，务必在使用 R 包中的功能之前，先检查 Python 领域的最新进展。

Bioconductor 项目中有许多适用于生物信息学的 R 包（www.bioconductor.org/）。这应该是你在 R 世界中进行生物信息学功能开发的首选。不过，需要注意的是，许多 R 生物信息学包并不在 Bioconductor 上，因此务必在 综合 R 存档网络 (CRAN) 中查找更广泛的 R 包（请参阅 CRAN：cran.rproject.org/）。

Python 中有许多绘图库。Matplotlib 是最常用的库，但你也可以选择其他众多的库。在 R 的背景下，值得注意的是，Python 中有一个类似于 ggplot2 的实现，它基于图形语法描述语言，用于图表的绘制，惊讶吗？这就是 ggplot！（yhat.github.io/ggpy/）。

另见

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

有许多关于 R 的教程和书籍；请访问 R 的官方网站 (www.r-project.org/) 查看文档。
对于 Bioconductor，请查看文档：manuals.bioinformatics.ucr.edu/home/R_BioCondManual。
如果你从事 NGS 工作，可能还想查看 Bioconductor 中的高通量序列分析：manuals.bioinformatics.ucr.edu/home/ht-seq。
rpy库的文档是你通往 R 的 Python 接口，文档可以在rpy2.bitbucket.io/找到。
图形语法方法在 Leland Wilkinson 所著的书《图形语法》中有详细描述，由 Springer 出版。
在数据结构方面，可以在pandas库中找到类似于 R 的功能。你可以在pandas.pydata.org/pandas-docs/dev/tutorials.xhtml找到一些教程。Wes McKinney 的《Python 数据分析》一书（由 O'Reilly Media 出版）也是一个可以考虑的替代方案。在下一章中，我们将讨论 pandas，并在整本书中使用它。

在 Jupyter 中执行 R 魔法

Jupyter 相比标准 Python 提供了许多额外的功能。在这些功能中，它提供了一种名为魔法的可扩展命令框架（实际上，这只适用于 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 魔法的简洁性和优雅性：

首先需要做的是加载 R 魔法和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被用来查看背后的实现。

让我们读取在之前配方中下载的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 代码（注意双 %%）。

现在，我们可以将变量传递到 Python 命名空间：

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

DataFrame 的类型不是标准的 Python 对象，而是一个 pandas DataFrame。这与之前版本的 R magic 接口有所不同。

由于我们有一个 pandas DataFrame，我们可以通过 pandas 接口轻松操作它：

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

让我们将这个 DataFrame 重新放入 R 命名空间，如下所示：
```
%R -i seq_data
%R print(colnames(seq_data))
```

-i 参数通知 magic 系统，将后面的变量从 Python 空间复制到 R 命名空间。第二行仅显示 DataFrame 确实在 R 中可用。我们使用的名称与原始名称不同——它是 seq_data，而不是 seq.data。

让我们进行最后的清理（有关更多详细信息，请参见之前的食谱），并打印出与之前相同的条形图：

%%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 函数，因为 magic 系统会为您处理这些。当您告诉 R 打印图表时，它会神奇地出现在您的笔记本或图形控制台中。

还有更多内容...

随着时间的推移，R magic 的接口似乎发生了很大的变化。例如，我已经多次更新了本书第一版中的 R 代码。当前版本的 DataFrame 赋值返回 pandas 对象，这是一个重大变化。

另请参见

欲了解更多信息，请查看以下链接：

有关 IPython 魔法命令的基本说明，请参见 ipython.readthedocs.io/en/stable/interactive/magics.xhtml。
包括魔法命令在内的 IPython 第三方扩展的列表可以在 github.com/ipython/ipython/wiki/Extensions-Index 找到。

第三章：了解 NumPy、pandas、Arrow 和 Matplotlib

Python 的最大优势之一是其丰富的高质量科学和数据处理库。所有这些库的核心是NumPy，它提供了高效的数组和矩阵支持。在 NumPy 之上，我们几乎可以找到所有的科学库。例如，在我们这一领域，有Biopython。但其他通用的数据分析库也可以在我们这一领域使用。例如，pandas 是处理表格数据的事实标准。最近，Apache Arrow 提供了一些 pandas 功能的高效实现，并且支持语言互操作性。最后，Matplotlib 是 Python 领域中最常见的绘图库，适用于科学计算。虽然这些库都是广泛应用的通用库，但它们对生物信息学处理至关重要，因此我们将在本章中学习它们。

我们将从 pandas 开始，因为它提供了一个高层次的库，具有非常广泛的实际应用性。然后，我们将介绍 Arrow，我们只在支持 pandas 的范围内使用它。接下来，我们将讨论 NumPy，这是几乎所有工作背后的驱动力。最后，我们将介绍 Matplotlib。

我们的教程非常基础——这些库中的每一个都可以轻松占据一本完整的书，但这些教程应该足够帮助你完成本书的内容。如果你使用 Docker，并且由于所有这些库对于数据分析至关重要，它们可以在来自第一章的 tiagoantao/bioinformatics_base Docker 镜像中找到。

在本章中，我们将涵盖以下教程：

使用 pandas 处理疫苗不良事件
处理 pandas DataFrame 合并的陷阱
降低 pandas DataFrame 的内存使用
使用 Apache Arrow 加速 pandas 处理
理解 NumPy 作为 Python 数据科学和生物信息学的引擎
介绍 Matplotlib 用于图表生成

使用 pandas 处理疫苗不良事件

我们将通过一个具体的生物信息学数据分析示例来介绍 pandas：我们将研究来自疫苗不良事件报告系统（VAERS，vaers.hhs.gov/）的数据。VAERS 由美国卫生与公共服务部维护，包含自 1990 年以来的疫苗不良事件数据库。

VAERS 提供的数据是逗号分隔值（CSV）格式。CSV 格式非常简单，甚至可以用简单的文本编辑器打开（请注意，文件过大会导致编辑器崩溃），或者使用类似 Excel 的电子表格程序打开。pandas 可以非常轻松地处理这种格式。

准备工作

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

随时可以使用文本编辑器查看文件，或者最好使用诸如less（压缩文件用zless）的工具。您可以在vaers.hhs.gov/docs/VAERSDataUseGuide_en_September2021.pdf找到文件内容的文档。

如果您使用的是笔记本，代码已在开头提供，您可以处理所需的处理步骤。如果您使用的是 Docker，基础镜像已足够。

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

如何做到这一点...

请按照以下步骤操作：

让我们从加载主要数据文件并收集基本统计信息开始：

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 列。您可以使用上述任何一种策略来获取有关表格元数据的信息。

现在，让我们探索数据：

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

现在，让我们做一些基本的计算，即计算数据集中最大年龄：
```
vdata["AGE_YRS"].max()
vdata.AGE_YRS.max()
```

最大值为 119 岁。比结果更重要的是，注意两种访问AGE_YRS的方式（作为字典键和作为对象字段）来访问列。

现在，让我们绘制涉及的年龄分布：

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

这将生成两个图表（以下步骤显示的是简化版本）。我们在这里使用的是 pandas 的绘图工具，它底层使用 Matplotlib。

虽然我们已经有完整的 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 – 左侧 – 每个不良反应观察的年龄；右侧 – 显示年龄分布的直方图

我们也可以采取一种非图形的、更分析性的方法，比如按年计数事件：
```
vdata["AGE_YRS"].dropna().apply(lambda x: int(x)).value_counts()
```

输出结果如下：

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

现在，让我们看看有多少人死亡：

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。一般来说，我们必须小心解读 NaN。它可能表示 False，或者像大多数情况一样，仅表示数据缺失。如果是这种情况，它不应该被转换为 False。

现在，让我们将死亡的个人数据与所涉及的疫苗类型进行关联：

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)

获取仅包含死亡数据的 DataFrame 后，我们需要读取包含疫苗信息的数据。首先，我们需要进行一些关于疫苗类型及其不良反应的探索性分析。以下是简化后的输出：

           …
HPV9         1506
FLU4         3342
UNK          7941
VARZOS      11034
COVID19    648723

之后，我们必须选择与 COVID 相关的疫苗，并将其与个人数据进行合并。

最后，让我们看看前 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。你可以从主要的用户指南开始，网址为 pandas.pydata.org/docs/user_guide/index.xhtml。
如果你需要绘制数据图表，不要忘记查看指南中的可视化部分，因为它特别有帮助：pandas.pydata.org/docs/user_guide/visualization.xhtml。

处理连接 pandas DataFrame 时的陷阱

上一个食谱是对 pandas 的快速介绍，涵盖了我们在本书中将使用的大部分功能。尽管关于 pandas 的详细讨论需要一本完整的书，但在本食谱（以及下一个食谱）中，我们将讨论一些对数据分析有影响的主题，这些主题在文献中很少讨论，但却非常重要。

在这个食谱中，我们将讨论通过连接（joins）关联 DataFrame 时的一些陷阱：事实证明，许多数据分析错误是由于不小心连接数据所引入的。我们将在这里介绍一些减少此类问题的技巧。

准备工作

我们将使用与上一食谱相同的数据，但会稍微打乱它，以便讨论典型的数据分析陷阱。我们将再次将主要的不良事件表与疫苗表连接，但会从每个表中随机采样 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 中找到。

如何操作...

请按照以下步骤操作：

让我们先对个体数据和疫苗数据表做一个内连接：

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。这表明一些个体数据和疫苗数据没有被捕获。

让我们通过以下连接查找未被捕获的数据：

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 行疫苗数据。

还有其他方式可以连接数据。默认的方法是执行左外连接：

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 都有一条记录。当你进行左连接时，可能会遇到左侧数据被重复多次的情况。例如，我们之前做的 groupby 操作显示，VAERS_ID 为 962303 的数据有 11 条记录。这是正确的，但也不罕见的是，很多人错误地期望左侧的每一行在输出中仍然是单独一行。这是因为左连接会返回一个或多个左侧条目，而上述的内连接返回的是 0 或 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。做连接时，确保进行这样的检查。

最后，我们将重新审视有问题的 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 关联的批次。

还有更多...

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

确保你对连接操作有测试和断言：一个非常常见的 bug 是对连接行为的预期错误。你还应该确保在连接的列上没有空值，因为空值可能会产生大量多余的元组。

减少 pandas DataFrame 的内存使用

当你处理大量信息时——例如在分析全基因组测序数据时——内存使用可能会成为分析的限制因素。事实证明，天真的 pandas 在内存方面并不是很高效，我们可以大幅减少它的内存消耗。

在这个方案中，我们将重新审视我们的 VAERS 数据，并探讨几种减少 pandas 内存使用的方法。这些变化的影响可能是巨大的：在许多情况下，减少内存消耗可能意味着能否使用 pandas，或者需要采用更复杂的替代方法，如 Dask 或 Spark。

准备就绪

我们将使用第一个方案中的数据。如果你已经运行过它，你就可以开始了；如果没有，请按照其中讨论的步骤操作。你可以在Chapter02/Pandas_Memory.py找到这段代码。

如何操作...

按照这些步骤操作：

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

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

在这里，我们有关于行数、每行类型和非空值的相关信息。最后，我们可以看到，DataFrame 需要高达 1.3 GB 的内存。

我们还可以检查每一列的大小：

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。

现在，让我们看看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 倍！

那么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。

作为这种方法的替代方案，您可以查看 pandas 中的分类变量。我更喜欢使用它们，因为它们有更广泛的应用。

我们可以在加载数据时应用大多数这些优化，因此让我们为此做好准备。但现在，我们遇到了一个先有鸡还是先有蛋的问题：为了能够了解州表的内容，我们必须进行第一次遍历，获取州列表，如下所示：
```
states = list(pd.read_csv(
    "vdata_sample.csv.gz",
    converters={
       "STATE": lambda state: state.upper()
    },
    usecols=["STATE"]
)["STATE"].unique())
```

我们有一个转换器，简单地返回州的全大写版本。我们只返回 STATE 列，以节省内存和处理时间。最后，我们从 DataFrame 中获取 STATE 列（该列只有一个字段）。

最终的优化是不加载数据。假设我们不需要 SYMPTOM_TEXT —— 这大约占数据的三分之一。在这种情况下，我们可以跳过它。以下是最终版本：

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 的食谱，它将帮助您在更多内存效率上节省额外的内存。
如果最终得到的 DataFrame 占用了比单台机器可用内存更多的内存，那么您必须提高处理能力并使用分块处理——我们在 Pandas 上下文中不会涉及——或者使用可以自动处理大数据的工具。Dask（我们将在 第十一章 “*使用 Dask 和 Zarr 进行并行处理” 中讨论）允许您使用类似 pandas 的接口处理超大内存数据集。

使用 Apache Arrow 加速 pandas 处理

当处理大量数据时，例如全基因组测序，pandas 的速度较慢且占用内存较大。Apache Arrow 提供了几种 pandas 操作的更快且更节省内存的实现，并且可以与 pandas 进行互操作。

Apache Arrow 是由 pandas 的创始人 Wes McKinney 共同创立的一个项目，它有多个目标，包括以与语言无关的方式处理表格数据，这样可以实现语言间的互操作性，同时提供高效的内存和计算实现。在这里，我们将只关注第二部分：提高大数据处理的效率。我们将与 pandas 一起以集成的方式实现这一点。

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

准备工作

再次，我们将使用第一个食谱中的数据。确保你已经按照 准备工作 部分中解释的方式下载并准备好它，代码可以在 Chapter02/Arrow.py 中找到。

如何做到...

按照以下步骤进行：

让我们开始使用 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 中某些函数与 pandas 的名称相似（例如，read_csv），但这并不是最常见的情况。例如，注意我们计算 DataFrame 总大小的方法：通过获取每一列的大小并求和，这与 pandas 的方法不同。

让我们并排比较推断出的类型：

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 在类型推断方面通常更为具体，这也是其内存使用显著更低的主要原因之一。

现在，让我们做一个时间性能比较：

%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 的实现速度是三倍。由于这取决于硬件，因此你电脑上的结果可能有所不同。

让我们在不加载 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：少了四倍。

我们的目标是使用 Arrow 将数据加载到 pandas 中。为此，我们需要转换数据结构：
```
vdata = vdata_arrow.to_pandas()
vdata.info(memory_usage="deep")
```

这里有两点非常重要：Arrow 创建的 pandas 表示只用了 1 GB，而 pandas 从其本地的 read_csv 创建的表示则需要 1.3 GB。这意味着即使你使用 pandas 处理数据，Arrow 也可以首先创建一个更紧凑的表示。

上述代码存在一个内存消耗问题：当转换器运行时，它将需要内存来存储 pandas 和 Arrow 的两个表示，从而违背了使用更少内存的初衷。Arrow 可以在创建 pandas 版本的同时自我销毁其表示，从而解决这个问题。相关代码行是 vdata = vdata_arrow.to_pandas(self_destruct=True)。

还有更多...

如果你有一个非常大的 DataFrame，甚至在 Arrow 加载后也无法由 pandas 处理，那么或许 Arrow 可以完成所有处理，因为它也有计算引擎。尽管如此，Arrow 的计算引擎在写作时，功能上远不如 pandas 完善。记住，Arrow 还有许多其他功能，例如语言互操作性，但我们在本书中不会使用到这些。

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

即使你没有显式使用 NumPy，大多数分析都将使用 NumPy。NumPy 是一个数组操作库，它是 pandas、Matplotlib、Biopython、scikit-learn 等许多库背后的基础。虽然你在生物信息学工作中可能不需要直接使用 NumPy，但你应该了解它的存在，因为它几乎支持了你所做的一切，即使是通过其他库间接使用。

在这个例子中，我们将使用 VAERS 数据演示 NumPy 如何在我们使用的许多核心库中发挥作用。这是一个非常简要的 NumPy 入门介绍，目的是让你了解它的存在，并知道它几乎在所有东西背后。我们的示例将从五个美国州中提取不良反应案例的数量，并按年龄分组：0 至 19 岁、20 至 39 岁，直到 100 至 119 岁。

准备工作

再次提醒，我们将使用第一个例子中的数据，因此请确保数据可用。相关的代码可以在 Chapter02/NumPy.py 找到。

如何实现…

按照以下步骤操作：

我们首先通过 pandas 加载数据，并将数据减少到仅与前五个美国州相关：

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 – 美国具有最大不良反应数量的州

现在，让我们提取包含年龄和州数据的两个 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

请注意，pandas 背后的数据是 NumPy 数据（对于 Series，values 调用返回的是 NumPy 类型）。此外，你可能还记得 pandas 有 .shape 或 .dtype 等属性：这些属性的设计灵感来源于 NumPy，且行为相同。

现在，让我们从头开始创建一个 NumPy 矩阵（一个二维数组），其中每一行代表一个州，每一列代表一个年龄组：

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 的部分非常明确。

我们可以非常轻松地提取单一的行——在我们这个例子中，是某个州的数据。对列同样适用。我们来看看加利福尼亚州的数据，然后是 0-19 岁的年龄组：
```
cal = age_state_mat[0,:]
kids = age_state_mat[:,0]
```

注意提取行或列的语法。鉴于 pandas 复制了 NumPy 的语法，且我们在之前的例子中也遇到过，它应该对你来说很熟悉。

现在，让我们计算一个新的矩阵，矩阵中的每个元素表示每个年龄组的案例比例：

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 接受一行数据并返回一个新行，其中所有元素都被总和除以。

现在，让我们创建一个新的矩阵，表示百分比而不是比例——这样看起来更清晰：

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 – 一个矩阵，表示美国五个病例最多的州中疫苗不良反应的分布

最后，让我们使用 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 就像 pandas 一样，是基于 NumPy 的。

另见

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

NumPy 有许多功能超出了我们在这里讨论的范围。市面上有许多书籍和教程可以帮助学习这些功能。官方文档是一个很好的起点：numpy.org/doc/stable/。
NumPy 有许多重要的功能值得探索，但其中最重要的可能是广播：NumPy 能够处理不同结构的数组，并正确地进行操作。详细信息请参见 numpy.org/doc/stable/user/theory.broadcasting.xhtml。

介绍用于图表生成的 Matplotlib

Matplotlib 是最常用的 Python 图表生成库。虽然也有一些更现代的替代库，如Bokeh，它是以 web 为中心的，但 Matplotlib 的优势不仅在于它是最广泛可用且文档最为丰富的图表库，还因为在计算生物学领域，我们需要一个既适用于 web 又适用于纸质文档的图表库。这是因为我们许多图表将会提交给科学期刊，这些期刊同样关注这两种格式。Matplotlib 能够为我们处理这一需求。

本示例中的许多例子也可以直接使用 pandas（间接使用 Matplotlib）完成，但这里的重点是练习使用 Matplotlib。

我们将再次使用 VAERS 数据来绘制有关 DataFrame 元数据的信息，并总结流行病学数据。

准备工作

再次，我们将使用第一个示例中的数据。代码可以在 Chapter02/Matplotlib.py 中找到。

如何做到…

按照以下步骤操作：

我们要做的第一件事是绘制每列空值的比例：

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 是接下来要绘制的条形图上条形的位置。

以下是第一个版本的条形图代码：

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 – 我们的第一次图表尝试，使用默认设置

当然，我们可以做得更好。让我们对图表进行更大幅度的格式化：

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 – 我们的第二次图表尝试，已考虑布局问题

现在，我们将根据一个图形上的四个子图来对数据进行一些汇总分析。我们将展示与死亡相关的疫苗、接种与死亡之间的天数、随时间变化的死亡情况以及十大州的死亡者性别：

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，并为绘图活动做了准备。

以下代码同时绘制所有信息。我们将有四个子图，按 2x2 格式排列：

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

我们首先创建一个 2x2 的子图图形。subplots 函数返回图形对象及四个坐标轴对象，我们可以使用这些坐标轴对象来创建我们的图表。请注意，图例位于饼图中，我们在时间距离图上使用了双坐标轴，并且在每个州的死亡率图上计算了堆积条形图。以下是结果：

图 2.6 – 汇总疫苗数据的四个合并图表

还有更多...

Matplotlib 有两个接口可供使用——一个较旧的接口，设计上类似于 MATLAB，另一个是功能更强大的 matplotlib.pyplot 模块。为了增加混淆，面向对象接口的入口点就在这个模块中——即 matplotlib.pyplot.figure 和 matplotlib.pyplot.subplots。

另见

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

Matplotlib 的文档真的非常出色。例如，它提供了一个包含每个示例代码链接的可视化样本画廊。你可以在matplotlib.org/stable/gallery/index.xhtml找到这个画廊。API 文档通常非常完整。
改善 Matplotlib 图表外观的另一种方法是使用 Seaborn 库。Seaborn 的主要目的是添加统计可视化工具，但作为副作用，它在导入时会将 Matplotlib 的默认设置更改为更易接受的样式。我们将在本书中始终使用 Seaborn；请查看下一章中提供的图表。

第四章：下一代测序

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

在这里，我们不会讨论每一种单独的 NGS 技术本身（这将需要一本完整的书）。我们将使用现有的 WGS 数据集——千人基因组计划，来说明分析基因组数据所需的最常见步骤。这里提供的步骤可以轻松应用于其他基因组测序方法。其中一些也可以用于转录组分析（例如，RNA-Seq）。这些步骤也与物种无关，因此你可以将它们应用于任何其他已测序物种的数据。来自不同物种的数据处理之间最大的差异与基因组大小、多样性以及参考基因组的质量（如果存在的话）有关。这些因素不会对 NGS 处理的自动化 Python 部分产生太大影响。无论如何，我们将在第五章中讨论不同的基因组，基因组分析。

由于这不是一本入门书籍，你需要至少了解FASTA（FASTA）、FASTQ、二进制比对映射（BAM）和变异调用格式（VCF）文件是什么。我还将使用一些基本的基因组术语而不做介绍（例如外显子、非同义突变等）。你需要具备基本的 Python 知识。我们将利用这些知识来介绍 Python 中进行 NGS 分析的基本库。在这里，我们将遵循标准的生物信息学流程。

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

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

在本章中，我们将涵盖以下步骤：

访问 GenBank 并浏览国家生物技术信息中心（NCBI）数据库
执行基本的序列分析
处理现代序列格式
处理比对数据
从 VCF 文件中提取数据
研究基因组可访问性并筛选单核苷酸多态性（SNP）数据
使用 HTSeq 处理 NGS 数据

访问 GenBank 并在 NCBI 数据库中浏览

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

由于测序分析是一个庞大的主题，而本书面向的是中高级用户，我们不会对这个本质上并不复杂的主题进行非常详尽的讲解。

尽管如此，这也是一个很好的热身练习，为我们在本章末尾看到的更复杂的教程做准备。

准备工作

我们将使用你在第一章中安装的 Biopython，Python 与周围的软件生态。Biopython 为Entrez提供了接口，Entrez是 NCBI 提供的数据检索系统。

这个教程可以在Chapter03/Accessing_Databases.py文件中找到。

小贴士

你将访问 NCBI 的实时应用程序编程接口（API）。请注意，系统的性能可能会在一天中有所波动。此外，使用时需要保持“良好的公民行为”。你可以在www.ncbi.nlm.nih.gov/books/NBK25497/#chapter2.Usage_Guidelines_and_Requiremen找到一些使用建议。特别需要注意的是，查询时你需要指定一个电子邮件地址。你应尽量避免在高峰时段（工作日美国东部时间上午 9:00 到下午 5:00 之间）发出大量请求（100 个或更多），并且每秒不超过三次查询（Biopython 会为你处理这个问题）。这不仅是为了遵守规定，而且如果你过度使用 NCBI 的服务器，你可能会被封锁（这也是一个好理由提供真实的电子邮件地址，因为 NCBI 可能会联系你）。

如何操作...

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

我们将从导入相关模块并配置电子邮件地址开始：

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

我们还将导入用于处理序列的模块。请不要忘记填写正确的电子邮件地址。

我们现在将尝试在nucleotide数据库中查找Plasmodium falciparum（导致最致命疟疾的寄生虫）：

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)

我们将搜索 nucleotide 数据库中的基因和生物体（关于搜索字符串的语法，参考 NCBI 网站）。然后，我们将读取返回的结果。请注意，标准搜索会将记录引用的数量限制为 20，因此如果有更多记录，可能需要重复查询并增加最大限制。在我们的例子中，我们将使用 retmax 来覆盖默认限制。Entrez 系统提供了多种复杂的方法来检索大量结果（更多信息请参考 Biopython 或 NCBI Entrez 文档）。尽管现在你已经拥有了所有记录的标识符（ID），但你仍然需要正确地检索记录。

现在，让我们尝试检索所有这些记录。以下查询将从 GenBank 下载所有匹配的核苷酸序列，截至本书编写时，共有 1,374 条。你可能不想每次都这样做：
```
id_list = rec_list['IdList']
hdl = Entrez.efetch(db='nucleotide', id=id_list, rettype='gb')
```

好吧，在这种情况下，继续进行吧。然而，使用这种技术时要小心，因为你将检索大量完整的记录，其中一些记录包含相当大的序列。你有可能下载大量数据（这会对你的机器和 NCBI 服务器都造成压力）。

有几种方法可以解决这个问题。一个方法是进行更严格的查询和/或每次只下载少量记录，当你找到需要的记录时就停止。具体策略取决于你的目标。无论如何，我们将检索 GenBank 格式的记录列表（该格式包含序列以及大量有趣的元数据）。

让我们读取并解析结果：
```
recs = list(SeqIO.parse(hdl, 'gb'))
```

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

如果你打算进行多次迭代，这可以节省时间、带宽和服务器使用。缺点是它会为所有记录分配内存。对于非常大的数据集，这种方法不可行；你可能不想像在 第五章 《处理基因组》那样进行全基因组的转换。如果你正在进行交互式计算，你可能更倾向于获取一个列表（以便你可以多次分析和实验），但如果你在开发一个库，迭代器可能是最佳的选择。

现在我们将集中精力处理单条记录。只有当你使用完全相同的前一个查询时，这种方法才有效：

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

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

现在，让我们提取一些包含 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 数组中检索序列，因此通常开始位置会比记录上的值小 1，结束位置将相等。坐标系统的问题将在未来的教程中再次讨论。

对于其他特征类型，我们只需打印它们。请注意，当你打印特征时，Biopython 会提供一个可读的人类版本。

现在我们来看一下记录上的注释，这些注释主要是与序列位置无关的元数据：
```
for name, value in rec.annotations.items():
    print('%s=%s' % (name, value))
```

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

最后但同样重要的是，你可以访问一个基础信息——序列：
```
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 之外还有许多其他生物数据库，如 Ensembl（www.ensembl.org）和加利福尼亚大学圣克鲁兹分校（UCSC）基因组生物信息学（genome.ucsc.edu/）。这些数据库在 Python 中的支持程度差异很大。

介绍生物数据库的食谱如果没有至少提到基础本地比对搜索工具（BLAST），那就不完整了。BLAST 是一种评估序列相似性的算法。NCBI 提供了一项服务，允许你将目标序列与其数据库进行比对。当然，你也可以使用本地 BLAST 数据库，而不是 NCBI 的服务。Biopython 为此提供了广泛的支持，但由于这是一个入门教程，我仅将你指引至 Biopython 的教程。

另见

这些附加信息也会很有用：

你可以在 Biopython 教程中找到更多示例，链接为biopython.org/DIST/docs/tutorial/Tutorial.xhtml。
可访问的 NCBI 数据库列表可以在www.ncbi.nlm.nih.gov/gquery/找到。
一个很棒的问答（Q&A）网站，你可以在上面找到关于数据库和序列分析问题的帮助，就是 Biostars (www.biostars.org)；你可以使用它来查找本书中的所有内容，而不仅仅是本章的这道示例。

执行基本的序列分析

现在我们将对 DNA 序列进行一些基本分析。我们将处理 FASTA 文件，并进行一些操作，例如反向互补或转录。与前面介绍的食谱一样，我们将使用你在第一章中安装的 Biopython，Python 与周边软件生态。这两道食谱为你提供了执行本章和第五章中所有现代 NGS 分析和基因组处理所需的入门基础。

准备工作

本食谱的代码可以在Chapter03/Basic_Sequence_Processing.py中找到。我们将使用人类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 序列记录中。

如何操作...

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

由于我们的目标序列已经存储在一个 Biopython 序列对象中，首先让我们把它保存为本地磁盘上的 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函数。

在大多数情况下，你实际上会把序列保存在磁盘上，因此你会对读取它感兴趣：

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]。序列对象可以使用典型的数组切片来获取序列的一部分。

既然我们现在有了明确无误的 DNA，我们可以按如下方式转录它：
```
rna = seq.transcribe()
print(rna)
```
最终，我们可以将我们的基因转化为蛋白质：
```
prot = seq.translate()
print(prot)
```

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

还有更多...

关于在 Biopython 中管理序列，还可以说很多，但这些主要是介绍性的内容，你可以在 Biopython 教程中找到。我认为给你们一个序列管理的概述是很重要的，主要是为了完整性。为了支持那些可能在生物信息学其他领域有一些经验，但刚刚开始做序列分析的人，尽管如此，仍然有几个点是你们应该注意的：

当你执行 RNA 翻译以获得你的蛋白质时，请确保使用正确的遗传密码。即使你正在处理“常见”的生物体（如人类），也要记住线粒体的遗传密码是不同的。
Biopython 的Seq对象比这里展示的要灵活得多。有一些很好的例子可以参考 Biopython 教程。然而，这个教程足以应付我们需要处理的 FASTQ 文件（请参见下一个教程）。
为了处理与链相关的问题，正如预期的那样，有一些序列函数，比如reverse_complement。
我们开始的 GenBank 记录包含了关于序列的大量元数据，所以一定要深入探索。

另请参见

Biopython 已知的遗传密码是由 NCBI 指定的，可以在www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi找到。
如同前面的教程一样，Biopython 教程是你的主要参考资料，可以在biopython.org/DIST/docs/tutorial/Tutorial.xhtml找到。
确保也查看 Biopython 的 SeqIO 页面，biopython.org/wiki/SeqIO。

使用现代序列格式

在这里，我们将处理 FASTQ 文件，这是现代测序仪使用的标准格式输出。你将学习如何处理每个碱基的质量分数，并考虑不同测序仪和数据库的输出变化。这是第一个将使用来自人类 1000 基因组计划的真实数据（大数据）的教程。我们将首先简要介绍这个项目。

准备工作

人类 1,000 基因组计划旨在 catalog 世界范围的人类遗传变异，并利用现代测序技术进行全基因组测序（WGS）。这个项目使所有数据公开可用，包括来自测序仪的输出、序列比对、SNP 调用以及许多其他数据。这一名称“1,000 基因组”实际上是一个误称，因为它目前包含了超过 2,500 个样本。这些样本分布在数百个人群中，遍布全球。我们将主要使用四个人群的数据：非洲约鲁巴人（YRI），拥有北欧和西欧血统的犹他州居民（CEU），东京的日本人（JPT），以及北京的汉族人（CHB）。之所以选择这些特定的人群，是因为它们是最早来自 HapMap 的样本，HapMap 是一个具有类似目标的旧项目。它们使用基因分型阵列来进一步了解该子集的质量。我们将在第六章“种群遗传学”中重新回顾 1,000 基因组和 HapMap 项目。

提示

下一代数据集通常非常大。由于我们将使用真实数据，你下载的某些文件可能会非常大。尽管我已经尽量选择最小的真实示例，但你仍然需要一个良好的网络连接以及相当大的磁盘空间。等待下载可能是你在这个配方中遇到的最大障碍，但数据管理是 NGS 中一个严重的问题。在现实生活中，你需要为数据传输预留时间，分配磁盘空间（这可能会涉及财务成本），并考虑备份策略。NGS 的最常见初始错误是认为这些问题微不足道，但事实并非如此。像将一组 BAM 文件复制到网络，甚至复制到你的计算机上，都会变成一件头疼的事。要做好准备。下载大文件后，至少应检查文件大小是否正确。一些数据库提供消息摘要 5（MD5）校验和。你可以通过使用像 md5sum 这样的工具，将这些校验和与下载文件中的校验和进行对比。

下载数据的指令位于笔记本的顶部，正如Chapter03/Working_with_FASTQ.py文件的第一单元格中所指定的。这是一个相当小的文件（27 NA18489）。如果你参考 1,000 基因组计划，你会发现绝大多数的 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），该分数为每个读取分配一个错误概率。这个编码在不同平台间可能略有不同。请确保在您的特定平台上检查此内容。

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

让我们打开文件：

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)

我们将打开一个gzip模块，并指定fastq格式。请注意，这种格式的某些变体会影响 Phred 质量分数的解释。您可能需要指定一个稍有不同的格式。有关所有格式的详细信息，请参考biopython.org/wiki/SeqIO。

提示

通常，您应该将 FASTQ 文件存储为压缩格式。这样不仅可以节省大量磁盘空间（因为这些是文本文件），而且可能还可以节省一些处理时间。尽管解压缩是一个较慢的过程，但它仍然可能比从磁盘读取一个更大的（未压缩的）文件更快。

我们将从之前的例子中打印标准字段和质量分数到rec.letter_annotations。只要我们选择正确的解析器，Biopython 会将所有 Phred 编码字母转换为对数分数，我们很快就会使用这些分数。

目前，不要这样做：

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

尽管这种方法可能适用于某些 FASTA 文件（以及这个非常小的 FASTQ 文件），但如果您这样做，您将分配内存以便将完整文件加载到内存中。对于一个普通的 FASTQ 文件，这是让您的计算机崩溃的最佳方式。一般来说，始终遍历您的文件。如果您需要对文件执行多次操作，您有两个主要选择。第一种选择是一次性执行所有操作。第二种选择是多次打开文件并重复迭代。

现在，让我们看看核苷酸读取的分布：

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调用的空间分布方面更多的信息。

让我们根据读取位置绘制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的分布作为序列开始位置的距离函数绘制出来：

图 3.1 – N调用的数量与序列读取起始位置的距离函数

图 3.1 – N调用的数量与序列读取起始位置的距离函数

你会看到，直到位置25，没有错误。这并不是你从典型的测序仪输出中得到的结果。我们的示例文件已经经过过滤，1,000 基因组过滤规则要求在位置25之前不能有N调用。

虽然我们不能在位置25之前研究N在这个数据集中的行为（如果你有一个未过滤的 FASTQ 文件，可以使用这段代码查看N在读取位置的分布），但我们可以看到，在位置25之后，分布远非均匀。在这里有一个重要的教训，那就是未调用的碱基数量是与位置相关的。那么，读取的质量如何呢？

让我们研究 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 质量分数是准确调用概率的对数表示。这个概率可以表示为！。因此，Q 为 10 表示 90% 的调用准确率，20 表示 99% 的调用准确率，30 表示 99.9%。对于我们的文件，最大准确率将是 99.99%（40）。在某些情况下，60 的值是可能的（99.9999%的准确率）。

更有趣的是，我们可以根据读取位置绘制质量分布：

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 – Phred 得分与测序读取起始位置距离的关系分布

还有更多内容...

尽管无法讨论来自测序仪文件的所有输出变化，但双端读取值得一提，因为它们很常见并且需要不同的处理方式。在双端测序中，DNA 片段的两端会被测序，并且中间有一个间隙（称为插入）。在这种情况下，将生成两个文件：X_1.FASTQ和X_2.FASTQ。这两个文件的顺序相同，且包含相同数量的序列。X_1中的第一个序列与X_2中的第一个序列配对，以此类推。关于编程技巧，如果你想保持配对信息，你可以执行如下操作：

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。

最后，如果你正在测序人类基因组，可能想要使用 Complete Genomics 的测序数据。在这种情况下，请阅读下一个食谱中的*更多内容…*部分，我们会简要讨论 Complete Genomics 数据。

另见

这里有一些提供更多信息的链接：

关于 FASTQ 格式的维基百科页面信息非常丰富（en.wikipedia.org/wiki/FASTQ_format）。
你可以在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 等人在《自然方法》期刊上发表的论文《使用下一代测序发现结构变异的计算方法》(www.nature.com/nmeth/journal/v6/n11s/abs/nmeth.1374.xhtml)；请注意，这篇文章并非开放获取。

处理对齐数据

在您从测序仪获得数据后，通常会使用像bwa这样的工具将您的序列与参考基因组进行比对。大多数用户会有自己物种的参考基因组。您可以在第五章，与基因组工作中了解更多关于参考基因组的信息。

对齐数据最常见的表示方法是 SAMtools 的tabix工具。SAMtools 可能是最广泛使用的 SAM/BAM 文件操作工具。

准备工作

如前一篇食谱中所讨论的，我们将使用来自 1,000 基因组计划的数据。我们将使用女性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 数据，检查类似类型的测序数据。

让我们看看以下步骤：

让我们检查一下头文件：

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 可能有多个实例，字典的值是一个列表（其中每个元素再次是一个字典，或者有时是包含标签/值对的字符串）。

现在，我们将检查一个单一的记录。每个记录的数据量相当复杂。在这里，我们将重点关注配对末端读取的一些基本字段。有关更多详细信息，请查看 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 文件对象可以通过其记录进行迭代。我们将遍历它，直到找到一个 Concise Idiosyncratic Gapped Alignment Report（CIGAR）字符串包含匹配和软剪切的记录。

CIGAR 字符串给出了单个碱基的比对信息。序列中被剪切的部分是比对工具未能对齐的部分（但未从序列中删除）。我们还需要读取序列、其配对 ID 以及映射到参考基因组上的位置（由于我们有配对末端读取，因此是配对的位置信息）。

首先，我们打印查询模板名称，接着是参考 ID。参考 ID 是指向给定参考序列查找表中序列名称的指针。一个示例可以使这点更加清晰。对于该 BAM 文件中的所有记录，参考 ID 是19（一个没有实际意义的数字），但如果你应用bam.getrname(19)，你会得到20，这就是染色体的名称。所以，不要将参考 ID（此处是19）与染色体名称（20）混淆。接下来是参考开始和参考结束。pysam是基于 0 的，而非基于 1 的，所以在将坐标转换为其他库时要小心。你会注意到，在这个例子中，开始和结束的位置分别是 59,996 和 60,048，这意味着一个 52 碱基的比对。为什么当读取大小是 76（记住，这是该 BAM 文件中使用的读取大小）时，只有 52 个碱基？答案可以通过 CIGAR 字符串找到，在我们的例子中是52M24S，即 52 个匹配碱基，后面是 24 个软剪接的碱基。

然后，我们打印比对的开始和结束位置，并计算其长度。顺便说一句，你可以通过查看 CIGAR 字符串来计算这一点。它从 0 开始（因为读取的第一部分已被比对），并在 52 处结束。长度再次是 76。

现在，我们查询配对端（如果你有双端读取，才会做这个操作）。我们获取它的参考 ID（如前面的代码片段所示），它的开始位置，以及两个配对之间的距离度量。这个距离度量只有在两个配对都映射到同一染色体时才有意义。

接着，我们绘制序列的 Phred 分数（参见之前的配方，处理现代序列格式，关于 Phred 分数的部分），然后仅绘制已比对部分的 Phred 分数。最后，我们打印出该序列（别忘了这么做！）。这是完整的序列，而不是剪接过的序列（当然，你可以使用之前的坐标来进行剪接）。

现在，让我们在 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 的记录（tabix）进行此类抓取操作；执行速度将会完全不同。

我们遍历所有位于 10 Mbp 边界内的记录。对于每个边界，我们获取比对的开始和结束，并增加在这些已比对位置中的映射性计数器。最后，我们将其转换为频率，然后进行绘制，如下图所示：

图 3.3 – 映射调用的百分比与测序器读取开始位置的函数关系

很明显，映射性分布远非均匀；在极端位置更差，中间位置则出现下降。

最后，让我们获取映射部分的 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 – Phred 分数在读取位置上的分布；底部蓝色表示从 0 到 5 百分位；绿色表示中位数，红色表示 95 百分位，紫色表示最大值

还有更多...

虽然我们将在本章的学习基因组可达性与过滤 SNP 数据配方中讨论数据过滤，但我们的目标并不是详细解释 SAM 格式或给出数据过滤的详细课程。这项任务需要一本专门的书籍，但凭借pysam的基础，你可以浏览 SAM/BAM 文件。不过，在本章的最后一个配方中，我们将探讨如何从 BAM 文件中提取全基因组度量（通过表示 BAM 文件度量的 VCF 文件注释），目的是了解我们数据集的整体质量。

你可能会有非常大的数据文件需要处理。有可能某些 BAM 处理会花费太多时间。减少计算时间的第一个方法是抽样。例如，如果你以 10% 进行抽样，你将忽略 10 条记录中的 9 条。对于许多任务，比如一些 BAM 文件质量评估的分析，以 10%（甚至 1%）进行抽样就足够获得文件质量的大致情况。

如果你使用的是人类数据，你可能会在 Complete Genomics 上进行测序。在这种情况下，比对文件将会有所不同。尽管 Complete Genomics 提供了将数据转换为标准格式的工具，但如果使用其自己的数据，可能会更适合你。

另见

额外的信息可以通过以下链接获得：

SAM/BAM 格式的描述可以在samtools.github.io/hts-specs/SAMv1.pdf找到。
你可以在 Abecasis 小组的维基页面上找到关于 SAM 格式的介绍：genome.sph.umich.edu/wiki/SAM。
如果你真的需要从 BAM 文件中提取复杂的统计数据，Alistair Miles 的pysamstats库是你的首选，网址：github.com/alimanfoo/pysamstats。
要将原始序列数据转换为比对数据，你需要一个比对工具；最常用的是 bwa，网址：bio-bwa.sourceforge.net/。
Picard（显然是指星际迷航：下一代中的人物）是最常用的工具来清理 BAM 文件；参考：broadinstitute.github.io/picard/。
序列分析的技术论坛叫做SEQanswers，网址：seqanswers.com/。
我想在这里重复一下 Biostars 上的推荐（在前面的“处理现代序列格式”一节中提到过）；这是一个信息宝库，且拥有一个非常友好的社区，网址：www.biostars.org/。
如果你有 Complete Genomics 数据，可以查看常见问题解答（FAQs），网址：www.completegenomics.com/customer-support/faqs/。

从 VCF 文件中提取数据

在运行基因型调用工具（例如，GATK 或 SAMtools）后，你将得到一个 VCF 文件，报告基因组变异信息，如 SNPs，cyvcf2模块。

准备工作

尽管 NGS 主要涉及大数据，但我不能要求你为本书下载过多的数据集。我认为 2 到 20 GB 的数据对于教程来说太多了。虽然 1,000 基因组计划的 VCF 文件和实际注释数据在这个数量级，但我们这里将处理更少的数据。幸运的是，生物信息学社区已经开发了工具，允许部分下载数据。作为 SAMtools/htslib包的一部分（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

第一行将部分下载来自 1,000 基因组计划的 22 号染色体 VCF 文件（最多 17 Mbp），然后，bgzip会进行压缩。

第二行将创建一个索引，这是我们直接访问基因组某一部分所需要的。像往常一样，你可以在笔记本中找到执行此操作的代码（Chapter03/Working_with_VCF.py文件）。

你需要安装cyvcf2：

conda install –c bioconda cyvcf2

提示

如果你遇到冲突解决问题，可以尝试改用 pip。这是一个最后的解决方案，当你使用 conda 时，往往会因为其无法解决软件包依赖问题而不得不这样做，你可以执行 pip install cyvcf2。

如何操作……

请查看以下步骤：

让我们从检查每条记录可以获取的信息开始：

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、INDEL 等，以及该变异在每个样本中的状态）。在变异（记录）级别，我们会找到 AC —— 在调用基因型中 ALT 等位基因的总数，AF —— 估算的等位基因频率，NS —— 有数据的样本数，AN —— 在调用基因型中的等位基因总数，以及 DP —— 总读取深度。还有其他信息，但它们大多数是 1000 基因组计划特有的（在这里，我们将尽量保持通用）。你自己的数据集可能有更多的注释（或没有这些注释）。

在样本级别，这个文件中只有两个注释：GT —— 基因型，和 DP —— 每个样本的读取深度。你有每个变异的总读取深度和每个样本的读取深度，请确保不要混淆两者。

现在我们知道了有哪些信息可用，让我们查看单个 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，例如只接受双等位基因 SNP），质量（如你所预期，采用 Phred 扩展评分），以及过滤状态。关于过滤状态，请记住，不论 VCF 文件中如何表示，你可能仍然需要应用额外的过滤器（如下一个教程所述，研究基因组可访问性和筛选 SNP 数据）。

然后，我们打印附加的变异级别信息（AC、AS、AF、AN、DP 等），接着是样本格式（在此案例中为 DP 和 GT）。最后，我们统计样本数并检查单个样本，以确认它是否针对该变异进行了调用。同时，还包括报告的等位基因、杂合性和相位状态（该数据集恰好是相位的，虽然这并不常见）。

让我们一次性检查变异的类型和非双等位基因 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 默认字典。我们发现该数据集中包含 INDEL、CNV 和——当然——SNP（大约三分之二是转换突变，一三分之一是倒位突变）。还有一个剩余的数量（79）为三等位基因 SNP。

还有更多内容……

本教程的目的是让你快速熟悉 cyvcf2 模块。此时，你应该已经能熟练使用该 API。我们不会花太多时间讲解使用细节，因为下一个教程的主要内容是：使用 VCF 模块研究变异调用的质量。

虽然cyvcf2非常快速，但处理基于文本的 VCF 文件仍然可能需要很长时间。有两种主要策略来处理这个问题。一种策略是并行处理，我们将在最后一章第九章，生物信息学管道中讨论。第二种策略是转换为更高效的格式；我们将在第六章，群体遗传学中提供示例。请注意，VCF 开发人员正在开发二进制变异调用格式（BCF）版本，以解决这些问题的部分内容（www.1000genomes.org/wiki/analysis/variant-call-format/bcf-binary-vcf-version-2）。

另见

一些有用的链接如下：

VCF 的规格可以在samtools.github.io/hts-specs/VCFv4.2.pdf找到。
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/howto-apply-hard-filters-to-a-call-set。

在第三种情况下，你有一个包含所有必要注释的 VCF 文件，但你可能希望应用更灵活的过滤器（例如，“如果读取深度 > 20，则映射质量 > 30 时接受；否则，映射质量 > 40 时接受”）。

在第四种情况下，你的 VCF 文件缺少所有必要的注释，你必须重新检查 BAM 文件（甚至其他信息来源）。在这种情况下，最好的解决方案是找到任何额外的信息，并创建一个带有所需注释的新 VCF 文件。一些基因型调用器（如 GATK）允许你指定需要哪些注释；你可能还需要使用额外的程序来提供更多的注释。例如，SnpEff (snpeff.sourceforge.net/) 会为你的 SNPs 注释其效应预测（例如，如果它们位于外显子中，它们是编码区还是非编码区？）。

提供一个明确的配方是不可能的，因为它会根据你的测序数据类型、研究物种以及你对错误的容忍度等变量有所不同。我们能做的是提供一套典型的高质量过滤分析方法。

在这个配方中，我们不会使用来自人类 1000 基因组计划的数据。我们想要的是脏的、未过滤的数据，其中有许多常见的注释可以用来进行过滤。我们将使用来自按蚊 1000 基因组计划的数据（按蚊是传播疟疾寄生虫的蚊子载体），该计划提供了过滤和未过滤的数据。你可以在www.malariagen.net/projects/vector/ag1000g上找到有关此项目的更多信息。

我们将获取大约 100 只蚊子的染色体3L的部分着丝粒区域，接着获取该染色体中间某部分（并索引两者）：

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/PacktPublishing/Bioinformatics-with-Python-Cookbook-third-edition/blob/main/Datasets.py以获取更新。像往常一样，用于下载该数据的代码在Chapter02/Filtering_SNPs.ipynb笔记本中。

最后，关于这个配方的一个警告：这里的 Python 难度会比平时稍微复杂一些。我们编写的越通用的代码，你就越容易将其重用于你的特定情况。我们将广泛使用函数式编程技术（lambda函数）和partial函数应用。

如何做...

看一下以下步骤：

让我们从绘制两个文件中基因组变异分布的图表开始：

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 Notebook，请记得删除第一行）。在我解释功能之前，请注意我们正在做什么。

对于这两个文件，我们将计算窗口统计数据。我们将我们的数据文件（包含 200,000 bp 的数据）分成大小为 2,000 的窗口（100 个窗口）。每次找到一个双等位基因 SNP，我们将在window函数相关的列表中添加一个 1。

window函数将获取一个 VCF 记录（rec.is_snp表示不是双等位基因的 SNP，长度为(rec.ALT) == 1），确定该记录所属的窗口（通过将rec.POS整除大小来执行），并通过作为fun参数传递给它的函数（在我们的情况下，只是 1）扩展该窗口结果列表。

因此，现在我们有了一个包含 100 个元素的列表（每个代表 2,000 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 个 1；因此，将它们求和将返回 200。

因此，我们能够以一种显然复杂的方式计算每个窗口中的 SNP 数量。为什么我们用这种策略执行事情很快就会显而易见。但是，现在，让我们检查这两个文件的计算结果，如下屏幕截图所示：

图 3.5 – 位于染色体 3L 附近 200 千碱基对（kbp）区域的 2000 bp 大小的双等位基因 SNP 分布窗口（橙色），以及染色体中部（蓝色）的情况；这两个区域来自约 100 只乌干达按蚊的 Anopheles 1,000 基因组计划

图 3.5 – 位于染色体 3L 附近 200 千碱基对（kbp）区域的 2000 bp 大小的双等位基因 SNP 分布窗口（橙色），以及染色体中部（蓝色）的情况；这两个区域来自约 100 只乌干达按蚊的 Anopheles 1,000 基因组计划。

小贴士

注意，在中心粒的 SNP 数量比染色体中部少。这是因为在染色体中调用变异体比在中部更困难。此外，中心粒的基因组多样性可能较少。如果您习惯于人类或其他哺乳动物，您会发现变异体的密度非常高——这就是蚊子的特点！

让我们来看看样本级别的注释。我们将检查映射质量零（请参阅 www.broadinstitute.org/gatk/guide/tooldocs/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 即可。如果您不习惯这种编程风格，它可能比较困难，但您很快就会发现它的好处。

现在，让我们打印每个窗口的中位数和第 75 百分位数（在本例中，窗口大小为 5,000）：

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 百分位数（虚线），这些 SNP 分布在每个窗口大小为 5,000 bp 的区域内，覆盖了位于着丝粒附近（蓝色）和染色体中部（绿色）的 200 kbp 区域；这两个区域来自 100 只左右乌干达蚊子的 3L 染色体，数据来自蚊子基因组 1,000 项目

图 3.6 – 样本 SNP 的 MQ0 中位数（实线）和第 75 百分位数（虚线），这些 SNP 分布在每个窗口大小为 5,000 bp 的区域内，覆盖了位于着丝粒附近（蓝色）和染色体中部（绿色）的 200 kbp 区域；这两个区域来自 100 只左右乌干达蚊子的 3L 染色体，数据来自蚊子基因组 1,000 项目。

对于 standard.vcf.gz 文件，中位数 MQ0 为 0（它绘制在图表底部，几乎不可见）。这是好的，表明大部分涉及变异调用的序列都能清晰地映射到基因组的这个区域。对于 centro.vcf.gz 文件，MQ0 的质量较差。此外，还有一些区域，基因型调用器无法找到任何变异（因此图表不完整）。

让我们将杂合性与DP（样本级注释）进行比较。在这里，我们将绘制每个 SNP 的杂合性调用比例与样本读取深度（DP）的关系图。首先，我们将解释结果，然后是生成它的代码。

下一张截图显示了在某一深度下杂合体调用的比例：

图 3.7 – 连续线表示在某一深度计算的杂合体调用比例；橙色区域为着丝粒区域；蓝色区域为“标准”区域；虚线表示每个深度的样本调用数；这两个区域来自于安哥拉疟蚊 1000 基因组计划中大约 100 只乌干达蚊子的 3L 号染色体

图 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 DataFrame，或者考虑使用 PyTables。这里的关键是要有一个足够通用的框架来计算样本注释之间的关系。

此外，注意多个注释的维度空间。例如，如果您的注释是浮动类型，您可能需要对其进行四舍五入（如果不进行处理，数据结构的大小可能会变得过大）。

现在，让我们来看一下绘图代码。我们分两部分来进行。以下是第一部分：

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 数量。

现在，让我们调用这个函数：

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

在这里，我们将使用两个坐标轴。左边是杂合 SNP 的比例，右边是 SNP 的数量。然后，我们为两个数据文件调用plot_hz_rel。其余部分是标准的matplotlib代码。

最后，让我们将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)

我们现在将遍历文件；样式现在应该对你清晰可见：

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

最后，我们使用 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')

在这里，我们仅为非着丝粒文件打印一个箱线图，如下面的图所示。结果符合预期：编码区域中的 SNP 可能会有更高的深度，因为它们位于更复杂的区域，这些区域比基因间的 SNP 更容易被调用：

图 3.8 – 不同 SNP 效应下变异读取深度的箱线图

还有更多内容...

关于过滤 SNP 和其他基因组特征的问题，需要一本书来详细讲解。这个方法将取决于你所拥有的测序数据类型、样本数量以及潜在的附加信息（例如，样本之间的家系关系）。

这个方法本身非常复杂，但其中一些部分非常简单（在一个简单的配方中，我无法强行加上过于复杂的内容）。例如，窗口代码不支持重叠窗口。并且，数据结构相对简单。然而，我希望它们能给你提供一种处理基因组高通量测序数据的整体策略。你可以在第四章，高级 NGS 处理中阅读更多内容。

另见

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

有许多过滤规则，但我想特别提醒你注意需要有足够好的覆盖度（显然需要超过 10x）。请参考*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 格式允许为注释轨道指定特征。它有许多用途，但常见的用途是将 BED 文件加载到基因组浏览器中以可视化特征。每行包含关于至少位置（染色体、起始和结束）的信息，也包括可选字段，如名称或链。关于该格式的详细信息可以在genome.ucsc.edu/FAQ/FAQformat.xhtml#format1找到。

准备工作

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

我们将从 Ensembl 获取这些信息。请访问uswest.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000115850，并选择LCT.bed文件，该文件位于Chapter03目录中。

这段代码的笔记本文件名为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）。

你需要安装 HTSeq：

conda install –c bioconda htseq

现在，我们可以开始了。

如何做到...

看一下以下步骤：

我们将首先为文件设置一个读取器。记住，这个文件已经提供给你，并且应该在你的当前工作目录中：
```
from collections import defaultdict
import re
import HTSeq
lct_bed = HTSeq.BED_Reader('LCT.bed')
```

现在我们将通过它们的名称提取所有类型的特征：

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

我们保存了最后一条记录，以便检查它：

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

让我们深入研究这个区间：

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

还有更多...

BED 格式可能比这更复杂。此外，前面的代码是基于我们的文件内容的特定前提。不过，这个例子应该足够让你入门。即使在最糟糕的情况下，BED 格式也不是特别复杂。

HTSeq 的功能远不止这些，但这个例子主要作为整个包的起点。HTSeq 具有可以替代我们到目前为止介绍的大部分配方的功能。

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