R-生物信息学秘籍-一-R 生物信息学秘籍（一）前言在《R 生物信息学宝典》中，您将通过实际的案例遇到生物信息学领域

R 生物信息学秘籍（一）

原文：annas-archive.org/md5/38cc3d9f06ec3561cf365e2fe73c5d5a

译者：飞龙

协议：CC BY-NC-SA 4.0

前言

在《R 生物信息学宝典》中，您将通过实际的案例遇到生物信息学领域中常见的以及不太常见的挑战。

本书将采用基于配方的方法，帮助您使用 R 在计算生物学领域进行实际的研究和分析。您将通过分析 Bioconductor、ggplot 和 tidyverse 库，深入理解您的数据，学习如何在 RNAseq、系统发育学、基因组学和序列分析等领域进行分析。您还将了解如何在生物信息学领域运用机器学习技术。您将发展关键的计算技能，例如在 R Markdown 中开发工作流，并设计您自己的包以便高效且可复现地重用代码。

到本书结束时，您将对生物信息学分析中最重要和最广泛使用的技术有扎实的理解，并掌握处理真实生物数据所需的工具。

本书适合的人群

本书适合数据科学家、生物信息学分析师、研究人员和 R 开发人员，他们希望通过基于配方的方法解决中级到高级的生物学和生物信息学问题。具备 R 编程语言的工作知识以及对生物信息学的基本了解是必须的。

如何最大程度地利用本书

本书要求具备良好的 R 语言和其各类库的知识。

下载示例代码文件

您可以从您的账户下载本书的示例代码文件，网址为 www.packt.com。如果您是从其他地方购买本书，可以访问 www.packtpub.com/support，并注册以便将文件直接通过电子邮件发送给您。

您可以按照以下步骤下载代码文件：

登录或在 www.packt.com 注册。
选择支持选项卡。
点击代码下载。
在搜索框中输入书名并按照屏幕上的指示操作。

下载文件后，请确保使用最新版本的以下工具解压文件：

适用于 Windows 的 WinRAR/7-Zip
适用于 Mac 的 Zipeg/iZip/UnRarX
适用于 Linux 的 7-Zip/PeaZip

本书的代码包也托管在 GitHub 上，地址为 github.com/PacktPublishing/R-Bioinformatics-Cookbook。如果代码有更新，将会在现有的 GitHub 仓库中进行更新。

我们还提供了来自丰富书籍和视频目录的其他代码包，您可以访问 github.com/PacktPublishing/。快去看看吧！

下载彩色图像

我们还提供了一份 PDF 文件，包含本书中使用的截图/图表的彩色图像。您可以在此下载： www.packtpub.com/sites/default/files/downloads/9781789950694_ColorImages.pdf。

使用的规范 packtpub.com/…/9781789950694_ColorImages.pdf

本书中使用了多种文本规范。

CodeInText：表示文本中的代码词汇、数据库表名、文件夹名称、文件名、文件扩展名、路径名、虚拟 URL、用户输入和 Twitter 账号。例如：“我们将讨论如何使用 ape 和 treeio 包将树形数据导入和导出 R。”

代码块如下设置：

if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install()

粗体：表示一个新术语、一个重要的词汇，或是您在屏幕上看到的单词。例如，菜单或对话框中的单词在文本中会这样显示。示例：“一些依赖项依赖于封装的 Java 代码，因此您需要为您的系统安装一个Java 运行时环境（JRE）。”

警告或重要说明将以这样的形式出现。

小贴士和技巧将以这样的形式出现。

章节

本书中有几个常见的标题（准备工作、如何操作…、它是如何工作的…、还有更多…、另见）。

为了清晰地说明如何完成一个食谱，请按照以下方式使用这些章节：

准备工作

本节告诉您在食谱中需要做的工作，并描述如何设置所需的任何软件或前期设置。

如何操作…

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

它是如何工作的…

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

还有更多…

本节提供有关食谱的附加信息，以帮助您更好地理解食谱。

另见

本节提供指向其他有用信息的链接，帮助您了解食谱。

联系我们

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

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

勘误：虽然我们已经尽一切努力确保内容的准确性，但错误确实会发生。如果您在本书中发现错误，我们将非常感激您能向我们报告。请访问 www.packtpub.com/support/err…，选择您的书籍，点击“勘误提交表单”链接并输入相关详情。

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

如果您有兴趣成为作者：如果您在某个领域有专长，并且有兴趣撰写或贡献书籍，请访问 authors.packtpub.com。

评审

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

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

第一章：执行定量 RNAseq

RNAseq 技术彻底改变了转录本丰度的研究，带来了高灵敏度检测和高通量分析。使用 RNAseq 数据的生物信息学分析管道通常从读取质量控制步骤开始，接着是将序列比对到参考基因组或将序列读取组装成更长的转录本de novo。之后，通过读取计数和统计模型来估算转录本丰度，并评估样本间的差异表达。自然，整个管道的各个步骤都有多种技术可供选择。质量控制和读取比对步骤通常会在 R 之外进行，因此在 R 中的分析将从转录本或基因注释文件（如 GFF 和 BED 文件）和比对后的读取文件（如 BAM 文件）开始。

R 中的分析工具功能强大且灵活。它们中的许多是 Bioconductor 套件的一部分，因此它们可以很好地集成在一起。研究人员希望通过 RNAseq 回答的核心问题通常是：哪些转录本有差异表达？在本章中，我们将探讨一些标准情况下的分析方法，假设我们已经知道感兴趣基因的基因组位置，并且在需要找到未注释转录本的情况下进行分析。我们还将探讨其他重要的分析方法，帮助回答问题多少重复实验足够？以及哪个等位基因表达更多？

在本章中，我们将涵盖以下分析方法：

使用 edgeR 估算差异表达
使用 DESeq2 估算差异表达
使用 powsimR 进行功效分析
使用 GRanges 对象查找未注释的转录区域
使用 bumphunter 查找初始显示高表达的区域
差异峰值分析
使用 SVA 估算批次效应
使用 AllelicImbalance 查找等位基因特异性表达
绘制和展示 RNAseq 数据

技术要求

你需要的示例数据可以从本书的 GitHub 仓库获取：github.com/PacktPublishing/R-Bioinformatics_Cookbook.如果你希望按照书中的代码示例使用数据，那么你需要确保这些数据在你工作目录的子目录中。

这里是你需要的 R 包。大多数可以通过install.packages()安装*；*其他一些则稍微复杂一些：

Bioconductor
- AllelicImbalance
- bumphunter
- csaw
- DESeq
- edgeR
- IRanges
- Rsamtools
- rtracklayer
- sva
- SummarizedExperiment
- VariantAnnotation
dplyr
extRemes
forcats
magrittr
powsimR
readr

Bioconductor非常庞大，并拥有自己的安装管理器。你可以通过以下代码安装它：

if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install()

更多信息可访问www.bioconductor.org/install/。

通常，在 R 中，用户会加载一个库并直接按名称使用其中的函数。这在交互式会话中很有用，但当加载了许多包时，它可能会导致混淆。为了澄清我在某个时刻使用的是哪个包和函数，我偶尔会使用packageName::functionName()这种约定。

有时，在一个代码片段的中间，我会暂停代码执行，以便你能看到一些中间输出或一个对象的结构，这是理解时很重要的内容。每当这种情况发生时，你会看到一个代码块，其中每一行都以##（双井号）开头。请看以下命令：

letters[1:5]

这将给我们以下输出：

## a b c d e

请注意，输出行前缀是##。

使用 edgeR 估计差异表达

edgeR 是一个广泛使用且功能强大的包，它实现了负二项模型，适用于稀疏计数数据（例如 RNAseq 数据），并且是一般线性模型框架中描述和理解计数关系以及多组实验的精确检验的强大工具。它使用了一种加权风格的标准化方法，称为 TMM，这是样本与对照之间对数比值的加权平均值，在去除具有高计数和异常对数比值的基因之后。TMM 值应该接近 1，但它可以作为一个修正因子，应用到整体文库大小的调整中。

在这个食谱中，我们将查看一些从准备注释区域的读取计数到识别基因组中差异表达特征的选项。通常，有一个上游步骤要求我们将高通量序列读取对齐到参考序列，并生成描述这些对齐的文件，例如.bam文件。准备好这些文件后，我们将启动 R 并开始分析。为了集中精力处理差异表达分析部分，我们将使用一个已经准备好的数据集，其中所有数据都已准备就绪。第八章，与数据库和远程数据源的工作，展示了如果你想了解如何做这一步，可以从原始数据开始到达这一阶段的方法。

由于有许多不同的工具和方法来获取这些读取的对齐信息，我们将从两种常见的输入对象类型开始这个过程。我们将使用一个计数表，就像我们从文本文件加载时会用到的那样，并且我们将使用一个ExpressionSet（eset）对象，这是 Bioconductor 中常用的对象类型。

我们准备好的数据集将是来自 NHGRI DNA 元素百科全书项目的modencodefly数据集，研究的模式生物是果蝇（Drosophila melanogaster）。你可以在www.modencode.org上阅读关于这个项目的更多信息。该数据集包含 147 个D. melanogaster的不同样本，果蝇的基因组约为 110 Mbp，注释了大约 15,000 个基因特征。

准备就绪

数据提供了计数矩阵和 ExpressionSet 对象，你可以查看本书末尾的附录，了解更多关于这些对象类型的信息。数据存储在本书的代码和数据仓库中，网址是 github.com/PacktPublishing/R_Bioinformatics_Cookbook，位于 datasets/ch1/modencodefly_eset.RData、datasets/ch1/modencodefly_count_table.txt 和 datasets/ch1/modencodelfy_phenodata.txt。我们还将使用 edgeR（来自 Bioconductor）、readr 和 magrittr 库。

如何操作...

我们将展示使用 edgeR 估计差异表达的两种方法。

从计数表中使用 edgeR

对于从计数表（例如文本文件）使用 edgeR 估计差异表达，我们将使用以下步骤：

加载计数数据：

count_dataframe <- readr::read_tsv(file.path(getwd(), "datasets", "ch1", "modencodefly_count_table.txt" ))
genes <- count_dataframe[['gene']]
count_dataframe[['gene']] <- NULL
count_matrix <- as.matrix(count_dataframe)
rownames(count_matrix) <- genes
pheno_data <- readr::read_table2(file.path(getwd(), "datasets", "ch1", "modencodefly_phenodata.txt"))

指定感兴趣的实验：

experiments_of_interest <- c("L1Larvae", "L2Larvae")
columns_of_interest <- which( pheno_data[['stage']] %in% experiments_of_interest )

构建分组因子：

library(magrittr)
grouping <- pheno_data[['stage']][columns_of_interest] %>% 
forcats::as_factor()

构建计数数据的子集：

counts_of_interest <-  count_matrix[,columns_of_interest]

创建 DGE 对象：

library(edgeR)
count_dge <- edgeR::DGEList(counts = counts_of_interest, group = grouping)

执行差异表达分析：

design <- model.matrix(~ grouping)
eset_dge <- edgeR::estimateDisp(eset_dge, design)
fit <- edgeR::glmQLFit(eset_dge, design)
result <- edgeR::glmQLFTest(fit, coef=2)
topTags(result)

从 ExpressionSet 对象中使用 edgeR

使用 eset 对象通过 edgeR 进行估计可以按照以下步骤完成：

加载 eset 数据：

load(file.path(getwd(), "datasets/ch1/modencodefly_eset.RData"))

指定感兴趣的实验：

experiments_of_interest <- c("L1Larvae", "L2Larvae")
columns_of_interest <- which( phenoData(modencodefly.eset)[['stage']] %in% experiments_of_interest )

构建分组因子：

grouping <- droplevels(phenoData(modencodefly.eset)[['stage']][columns_of_interest] )

构建计数数据的子集：

counts_of_interest <- exprs(modencodefly.eset)[, columns_of_interest]

创建 DGE 对象：

eset_dge <- edgeR::DGEList(
 counts = counts_of_interest,
 group = grouping 
 )

执行差异表达分析：

design <- model.matrix(~ grouping)
eset_dge <- edgeR::estimateDisp(eset_dge, design)

fit <- edgeR::glmQLFit(eset_dge, design)
result <- edgeR::glmQLFTest(fit, coef=2)
topTags(result)

如何工作...

我们展示了两种使用 edgeR 估计差异表达的方法。在本教程的前半部分，我们从文本文件中的数据开始，使用了 edgeR。

从计数表中使用 edgeR

在步骤 1中，我们使用 read_tsv() 函数（来自 readr 包）将制表符分隔的计数文本文件加载到一个名为 count_dataframe 的数据框中。然后，从中提取 'gene' 列到一个新变量 genes，并通过将其赋值为 NULL 从 count_dataframe 中删除。这一切的目的是为了便于使用基础的 as.matrix() 函数将其转换为 count_matrix 矩阵，并将基因信息作为 rownames 添加回来。最后，我们使用 readr read_table2() 函数从文件中加载所需的表型数据。

步骤 2 主要是确定在 count_matrix 中我们想要使用的列。我们定义了一个变量 experiments_of_interest，它包含了我们想要的列名，然后使用 %in% 运算符和 which() 函数创建一个与列数匹配的二进制向量。如果例如 columns_of_interest 向量的第三列为 TRUE，则表示该名称在 experiments_of_interest 变量中。

步骤 3 从加载 magrittr 包开始，以获取 %>% 操作符，这将允许管道传输。接着，我们使用 R 索引和二进制 columns_of_interest 因子来选择我们想要的列名，并将其传递给 forcats as_factor() 函数，得到一个因子对象作为我们的分组变量。样本分组信息本质上是一个因子，告诉我们哪些样本是同一样本的重复，这对于实验设计的描述非常重要。我们需要创建一个分组向量，其中每个索引都对应于计数表中的一列。因此，在下面的示例中，数据中的前三列将是同一样本的重复，计数表中的第二组三列将是另一组重复，以此类推。我们可以在分组向量中使用任何符号来表示组。分组向量越复杂，实验设计就可以越复杂。在这里的示例中，我们将使用一个简单的测试/对照设计：

numeric_groups <- c(1,1,1,2,2,2)
letter_groups <- c("A","A","A", "B","B","B")

像这样的简单向量就足够了，但你也可以使用一个因子对象。因子是 R 的分类数据类型，它作为一个整数向量实现，并且有相关联的名称标签，称为水平。当显示因子时，使用的是名称标签，而不是整数。因子对象有某种“记忆”，即使只使用了某些水平，所有可能使用的水平都会被保留，这样在例如将水平作为类别使用时，空的水平仍然可以显示出来。

在 步骤 4 中，我们使用索引提取我们实际想要分析的数据列。

到了 步骤 5，我们的准备工作完成，可以构建我们需要进行差异分析的 DGEList 对象。首先，我们加载 edgeR 库，并对 counts_of_interest 和我们的分组对象使用 DGEList() 函数。

在 步骤 6 中，通过 DGEList，我们可以进行 edgeR 过程。首先，我们使用基础的 model.matrix() 函数创建实验设计描述符设计对象。模型设计是必要的，它告诉函数如何比较样本；这是 R 中常见的操作，因此有一个基础函数。我们使用我们创建的 grouping 变量。我们必须使用 estimateDisp() 函数估算每个基因的离散度，然后可以在测试中使用这个变异性度量。最后，拟合一个广义线性模型，并使用两次 glmQLFTest() 应用准似然 F 检验，第一次使用离散度估算，eset_dge，第二次使用结果的 fit 对象。

我们可以使用 topTags() 函数查看差异表达基因的详细信息。我们得到以下输出：

 ## Coefficient: groupingL2Larvae
 ## logFC logCPM F PValue FDR
 ## FBgn0027527 6.318665 11.14876 42854.72 1.132951e-41 1.684584e-37
 ## [ reached 'max' / getOption("max.print") -- omitted 9 rows ]

列中显示了基因名、基因的 logFC 值、F 值、P 值和 假阳性率 (FDR)。通常，我们想要从中做出统计结论的列是 FDR。

使用 `edgeR` 从 ExpressionSet 对象中

在步骤 1中，我们将使用从准备好的eset对象中的edgeR。我们首先使用基础 R 函数加载该对象，因为它是以标准的 Rdata 格式文件存储的。

在步骤 2中，我们准备感兴趣实验的向量。与步骤 2中类似，唯一不同的是，我们不需要查看从文件创建的pheno_data对象；相反，我们可以使用eset函数phenoData()直接从eset对象中提取表型数据（请注意，这是eset与计数矩阵之间的主要区别之一——请参阅本书的附录以获取更多信息）。

在步骤 3中，我们创建分组因子。与之前一样，这可以通过使用phenoData()提取函数完成，但由于它返回的是一个因子，我们需要使用droplevels()函数删除未被选中的级别（有关因子对象的简要讨论，请参见前述方法中估计差异表达部分的步骤 3）。

在步骤 4中，我们将提取我们感兴趣的列的数据到标准的矩阵对象中。同样，我们有一个专用函数exprs()，用于从eset中提取表达值，并且我们可以使用column_names进行列索引来对子集进行操作。

在步骤 5中，我们使用DGEList()构造函数构建edgeR的数据结构，在步骤 6中，执行分析。这个步骤与第一种方法的步骤 6完全相同。

使用 DESeq2 估计差异表达

DESeq2包是用于计数数据差异分析的方法，因此它非常适合 RNAseq（以及其他计数类型数据，如ChIPSeq）。它使用离散度估计和相对表达变化来加强估计和建模，重点提高结果表中的基因排序。DESeq2与edgeR的不同之处在于，它使用几何风格的标准化方法，其中每条通道的缩放因子是通过基因计数与其几何均值比率的比值的中位数来计算的，而edgeR使用加权的标准化因子。两种标准化策略并不互相排斥，并且它们对数据有不同的假设。像任何RNAseq或大规模实验一样，永远没有“现成的”最佳答案。你将需要测试这些方法，甚至其他方法，并仔细检查控制基因和交叉验证实验的结果，以查看哪种方法表现最佳。性能将很大程度上取决于具体的数据集，因此我们在这里学习的灵活方法将为你如何自行测试不同的解决方案提供一个良好的思路。

我们将在本食谱中讨论的过程与之前食谱 1中的edgeR过程有些相似。我们可以将 ExpressionSets 和计数表作为输入提供给DESeq2，当我们准备好这些数据时，我们将有一组不同的函数来将数据转换为DESeqDataSet，而不是像edgeR那样使用DGEList。

准备工作

如同食谱 1中所示，数据以计数矩阵和ExpressionSet对象的形式提供，你可以在本书的附录部分找到有关这些对象类型的更多信息。数据位于本书的代码和数据仓库：github.com/PacktPublishing/R_Bioinformatics_Cookbook，路径为datasets/ch1/modencodefly_eset.RData，datasets/ch1/modencodefly_count_table.txt，以及datasets/ch1/modencodelfy_phenodata.txt。再次提醒，我们将使用readr和magrittr，以及 Bioconductor 中的SummarizedExperiement和 DESeq2。

如何操作……

使用 DESeq2 估计差异表达可以通过两种方式进行，如下所示。

从计数矩阵中使用 DESeq2

从计数表（例如，文本文件）中使用 DESeq2 估计差异表达，我们将使用以下步骤：

加载计数数据：

count_dataframe <- readr::read_tsv(file.path(getwd(), "datasets", "ch1", "modencodefly_count_table.txt" ))
genes <- count_dataframe[['gene']]
count_dataframe[['gene']] <- NULL
count_matrix <- as.matrix(count_dataframe)
rownames(count_matrix) <- genes
pheno_data <- readr::read_table2(file.path(getwd(), "datasets", "ch1", "modencodefly_phenodata.txt"))

指定感兴趣的实验：

experiments_of_interest <- c("L1Larvae", "L2Larvae")
columns_of_interest <- which( pheno_data[['stage']] %in% experiments_of_interest )

形成分组因子：

library(magrittr)
grouping <- pheno_data[['stage']][columns_of_interest] %>% 
 forcats::as_factor()

形成计数数据的子集：

counts_of_interest <- count_matrix[,columns_of_interest]

构建 DESeq 对象：

library("DESeq2")
dds <- DESeqDataSetFromMatrix(countData = counts_of_interest,
 colData = grouping,
 design = ~ stage)

执行分析：

dds <- DESeq(dds)

提取结果：

res <- results(dds, contrast=c("stage","L2Larvae","L1Larvae"))

从 ExpressionSet 对象中使用 DESeq2

要从 ExpressionSet 对象中使用 DESeq2 估计差异表达，我们将使用以下步骤：

加载eset数据并转换为DESeqDataSet()：

library(SummarizedExperiment)
load(file.path(getwd(), "datasets/ch1/modencodefly_eset.RData"))
summ_exp <- makeSummarizedExperimentFromExpressionSet(modencodefly.eset)
ddsSE <- DESeqDataSet(summ_exp, design= ~ stage)

执行分析并提取结果：

ddsSE <- DESeq(ddsSE)
resSE <- results(ddsSE, contrast=c("stage","L2Larvae","L1Larvae"))

它是如何工作的……

在本食谱的第一部分，我们使用 DESeq1 从文本文件中的数据开始；正如你会注意到的，步骤 1到步骤 4与前一部分完全相同。

从计数矩阵中使用 DESeq2

在步骤 1中，我们使用readr包中的read_tsv()函数加载带有制表符分隔的计数文本文件到名为count_dataframe的数据框中。然后，我们从中提取'gene'列到一个新变量genes，并通过赋值NULL将其从count_dataframe中删除。这些操作的目的是为了方便我们使用基础的as.matrix()函数将其转换为count_matrix矩阵，并将基因信息作为rownames重新添加进去。最后，我们使用readr的read_table2()函数从文件中加载所需的表型数据。

步骤 2的任务是确定我们想要在count_matrix中使用哪些列。我们定义一个变量experiments_of_interest，保存我们需要的列名，然后使用%in%操作符和which()函数创建一个与列数匹配的二进制向量。如果比如columns_of_interest向量的第三列为'TRUE'，这表示该列名存在于experiments_of_interest变量中。

第 3 步 从加载 magrittr 包开始，以获取 %>% 操作符，这将允许管道操作。然后我们使用 R 索引与二进制 columns_of_interest 因子选择我们需要的列名，并将其传递给 forcats 包中的 as_factor() 函数，以获得分组变量的因子对象。样本分组信息基本上是一个因子，告诉我们哪些样本是相同事物的重复，这对于实验设计描述非常重要。你可以在 步骤 3 中的 食谱 1 中看到对这些分组/因子对象的扩展描述。

在 第 4 步 中，我们使用索引提取我们实际需要分析的数据列。

到了 第 5 步，我们进入了实际的分析部分。首先，我们将计数矩阵转换为 DESeqDataSet 对象；这可以通过转换函数 DESeqDataSetFromMatrix() 来完成，传入计数、分组和设计。设计采用 R 公式的形式，因此使用了 ~ stage 注解。

在 第 6 步 中，我们使用 DESeq() 函数对 dds DESeqDataSet 对象进行实际分析，在 第 7 步 中，使用 results() 函数将结果存储到 res 变量中。输出有以下六列：

baseMean log2FoldChange lfcSE stat pvalue padj

这显示了基因的平均计数、样本之间的 log2 fold change、log2 fold change 的标准误差、Wald 统计量，以及原始和调整后的 P 值。padj 列表示调整后的 P 值，通常用于得出显著性的结论。

使用来自 ExpressionSet 对象的 DESeq2

第 1 步 和 第 2 步 显示了如何从 eset 对象开始进行相同的操作。只需两步，因为 DESeq2 在处理 Bioconductor 对象时比 edgeR 更加高效。在 第 8 步 中，我们使用 load() 函数加载 eset 数据。然后我们使用来自 SummarizedExperiment Bioconductor 包的 makeSummarizedExperimentFromExpressionSet() 函数，将 eset 转换为 SummarizedExperiment，这可以直接用于 第 9 步 中的 DESeq() 函数。此步骤与 第 6 步 和 第 7 步 完全相同。

使用 powsimR 进行功效分析

任何实验的一个重要前提是评估实验设计的功效，以优化统计敏感性。从本质上讲，功效分析可以告诉我们，在给定实验变异性下，为确定特定效应大小所需的样本重复次数。

我们将使用 powsimR 包，它不属于 Bioconductor，用来进行两种类型的功效分析。这两种分析都使用一个小的真实数据集，但首先，我们将用两种处理方法——测试组和对照组——进行分析，然后，只用一个处理组进行分析。在每个分析中，我们将估计需要多少个重复实验才能发现基因表达的特定差异——如果这些差异存在的话。powsimR 采用基于仿真的方法，实际生成多个数据集并评估每个数据集中的检测功效，最终形成一个功效分布。因此，第一步是为这些仿真估计一些参数——为此，我们需要一些样本数据或初步数据。之后，我们可以运行仿真并评估功效。

准备工作

这个配方的数据集将是来自 Arabidopsis 的一个测试或对照 RNAseq 实验，每组有三个重复。这些数据可以在本书的数据仓库中的 datasets/ch1/arabidopsis.RDS 中找到，作为一个预处理过的计数矩阵。在本节中，我们将使用一个来自 Arabidopsis thaliana 的简单测试或对照实验的计数数据集。该矩阵有六列（三个 mock 处理组和三个 hrcc 处理组），共 26,222 行，每行代表一个基因特征。我们需要 dplyr、extRemes 和 powsimR 包来运行这段代码。

我们感兴趣的包 powsimR 不在 CRAN 上，它作为源代码托管在 GitHub 上，网址为 github.com/bvieth/powsimR。你需要使用 devtools 来安装它，以下代码可以完成安装：

install.packages("devtools")
devtools::install_github("bvieth/powsimR")

如果你这样做，仍然有可能会安装失败。它有很多依赖项，你可能需要手动安装这些依赖；在该包的 GitHub 仓库上有更多信息，建议你查看最新的内容。写这篇文档时，你需要做以下两个步骤。首先，创建这里描述的 ipak 函数，然后使用 ipak 函数运行三个不同的包安装步骤：

ipak <- function(pkg, repository = c("CRAN", "Bioconductor", "github")) {
    new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
    # new.pkg <- pkg
    if (length(new.pkg)) {
        if (repository == "CRAN") {
            install.packages(new.pkg, dependencies = TRUE)
        }
        if (repository == "Bioconductor") {
            if (strsplit(version[["version.string"]], " ")[[1]][3] > "3.5.0") {
                if (!requireNamespace("BiocManager")) {
                  install.packages("BiocManager")
                }
                BiocManager::install(new.pkg, dependencies = TRUE, ask = FALSE)
            }
            if (strsplit(version[["version.string"]], " ")[[1]][3] < "3.5.0") {
                source("https://bioconductor.org/biocLite.R")
                biocLite(new.pkg, dependencies = TRUE, ask = FALSE)
            }
        }
        if (repository == "github") {
            devtools::install_github(new.pkg, build_vignettes = FALSE, force = FALSE, 
                dependencies = TRUE)
        }
    }
}

# CRAN PACKAGES
cranpackages <- c("broom", "cobs", "cowplot", "data.table", "devtools", "doParallel", 
    "dplyr", "drc", "DrImpute", "fastICA", "fitdistrplus", "foreach", "gamlss.dist", 
    "ggExtra", "ggplot2", "ggthemes", "grDevices", "glmnet", "grid", "gtools", 
    "Hmisc", "kernlab", "MASS", "MBESS", "matrixStats", "mclust", "methods", 
    "minpack.lm", "moments", "msir", "NBPSeq", "nonnest2", "parallel", "penalized", 
    "plyr", "pscl", "reshape2", "Rmagic", "rsvd", "Rtsne", "scales", "Seurat", 
    "snow", "stats", "tibble", "tidyr", "VGAM", "ZIM")

ipak(cranpackages, repository = "CRAN")

# BIOCONDUCTOR
biocpackages <- c("AnnotationDbi", "bayNorm", "baySeq", "Biobase", "BiocGenerics", 
    "BiocParallel", "DEDS", "DESeq2", "EBSeq", "edgeR", "IHW", "iCOBRA", "limma", 
    "Linnorm", "MAST", "monocle", "NOISeq", "qvalue", "ROTS", "RUVSeq", "S4Vectors", 
    "scater", "scDD", "scde", "scone", "scran", "SCnorm", "SingleCellExperiment", 
    "SummarizedExperiment", "zinbwave")
ipak(biocpackages, repository = "Bioconductor")

# GITHUB
githubpackages <- c("nghiavtr/BPSC", "cz-ye/DECENT", "mohuangx/SAVER", "statOmics/zingeR")
ipak(githubpackages, repository = "github")

完成这些步骤后，你应该能够用以下代码安装我们需要的包：

devtools::install_github("bvieth/powsimR", build_vignettes = TRUE, dependencies = FALSE)
library("powsimR")

目前，为了使其正常工作，你还需要手动加载 dplyr。

如何做...

我们将通过以下步骤进行功效分析：

估计仿真参数值：

arab_data <- readRDS(file.path(getwd(), "datasets", "ch1", "arabidopsis.RDS" ))
means_mock <- rowMeans(arab_data[, c("mock1", "mock2", "mock3")])
means_hrcc <- rowMeans(arab_data[, c("hrcc1", "hrcc2", "hrcc3")])
log2fc <- log2(means_hrcc / means_mock)
prop_de <- sum(abs(log2fc) > 2) / length(log2fc)

检查 log2 变动比率的分布：

finite_log2fc <-log2fc[is.finite(log2fc)]
plot(density(finite_log2fc))
extRemes::qqnorm(finite_log2fc )

设置仿真运行的参数值：

library(powsimR)
 library(dplyr)

 params <- estimateParam(
 countData = arab_data,
 Distribution = "NB",
 RNAseq = "bulk",
 normalization = "TMM" # edgeR method, can be others
 )

 de_opts <- DESetup(ngenes=1000,
 nsims=25,
 p.DE = prop_de,
 pLFC= finite_log2fc,
 sim.seed = 58673
 )

 sim_opts <- SimSetup(
 desetup = de_opts,
 params = params
 )

 num_replicates <- c(2, 3, 5, 8, 12,15)

运行仿真：

 simDE <- simulateDE(n1 = num_replicates,
 n2 = num_replicates,
 sim.settings = sim_opts,
 DEmethod = "edgeR-LRT",
 normalization = "TMM",
 verbose = FALSE)

运行仿真评估：

 evalDE <- evaluateDE(simRes = simDE,
 alpha.type = 'adjusted',
 MTC = 'BH',
 alpha.nominal = 0.1,
 stratify.by = 'mean',
 filter.by = 'none',
 strata.filtered = 1,
 target.by = 'lfc',
 delta = 0)

绘制评估图：

 plotEvalDE(evalRes = evalDE,
 rate='marginal',
 quick=FALSE, annot=TRUE)

它是如何工作的...

powsimR 中的功效分析要求我们做一些预分析，以便估算一些重要参数。为了执行基于仿真的功效分析，我们需要估计处理组之间的 log fold change 分布以及差异表达特征的比例。

在步骤 1中，我们将获得每个特征在两种处理条件下的平均计数。加载表达数据后，我们使用readRDS()函数，然后在某些列上使用rowMeans()函数，获得mock和hrcc1两种处理下每个基因的平均表达计数。接着，我们可以获得这些值的 log2 比率（只需将这两个向量相除，并在最后一行使用标准算术运算符计算出那些 log2 折叠变化大于 2 的值）。检查最终的prop_de变量时，会得到以下输出：

prop_de
## [1] 0.2001754

因此，大约 0.2 的特征在计数上变化了 log2 两倍。

步骤 2 查看基因表达比率的分布。我们首先从log2fc变量中移除非有限值。我们必须这样做，因为在计算比率时，会在 R 中生成Inf值；这种情况发生在分母（对照样本）的均值为零时。我们可以通过使用来自is.finite()函数的二进制向量对向量进行索引来移除这些值。移除Inf值后，我们可以进行绘图。首先，我们使用density()函数做一个正态密度图，展示比率的分布。然后，我们使用extRemes包中的qqnorm()函数，该函数将数据与来自理想正态分布、均值相同的采样数据进行对比。强烈的线性相关性表示原始数据符合正态分布。我们可以在以下截图中看到输出：

它们看起来呈对数正态分布，因此我们可以假设它们符合对数正态分布。

这里最长的步骤，步骤 3，实际上只有四行代码。我们基本上是在为模拟设置参数，这需要我们指定很多值。第一组params是通过estimateParam()函数创建的，需要提供数据源（countData）、要使用的分布（我们设置Distribution = "NB"，选择负二项分布）；RNAseq 实验的类型——我们的实验是批量 RNAseq 实验（RNAseq = "bulk"），以及归一化策略——我们使用 edgeR 风格的 TMM（normalization = "TMM"）。第二组desetup是通过DESetup()函数创建的；在这里，我们选择与模拟差异表达基因数量相关的参数。我们设置总共模拟 1,000 个基因（ngenes）和 25 次模拟运行（nsims）。我们将差异表达的比例设置为步骤 1中估算的值，存储在prop_de中。我们使用finite_log2fc的向量作为pLFC参数的输入。设置sim.seed不是必须的，但它可以确保不同运行之间的可重复性。第三行使用SimSetup()函数将params和desetup合并为一个单一对象sim_opts。最后，我们创建一个num_replicates向量，指定要模拟的生物学重复数（RNA 样本数）。

步骤 4相对简单：我们使用前几步创建的sim_opts参数运行差异表达仿真，选择"edgeR-LRT"作为差异表达方法，选择"TMM"作为标准化方法。仿真数据存储在simDE变量中。

在步骤 5中，我们创建了一个仿真评估——它分析并提取了各种统计数据。我们将simDE仿真数据传递给evaluateDE()函数，并附带与分组、过滤和显著性相关的值。

最后，在步骤 6中，我们可以绘制来自步骤 5的evalDE对象，并查看仿真结果。我们得到以下图表，其中可以看到不同重复次数下的不同功效。请注意，x轴表示使用的重复 RNA 样本数量，指标包括 FDR、假阴性/假阳性率 (FNR/FPR)，以及真阴性/真阳性率 (TNR/TPR)：

还有更多...

当我们只有一个样本（或者可能只有一个重复）时，我们很难估算 log2 倍数变化分布以及差异表达基因的数量。我们可以使用回调函数代替估算来生成所需的数字。该函数的主体只需返回一个指定分布中的数字，具体分布由你决定。在这里，我们将构建一个返回均值为 0，标准差为 2 的正态分布数字的函数。这反映了我们认为 log 倍数变化分布是正态分布，且具有这些参数。构建完函数后，它将在DESetup()函数中替代 log2 倍数变化的向量。至于差异表达基因的比例，我们只需要猜测或从我们已知的实验系统中获取估算值：

log2fc_func <- function(x){ rnorm(x, 0, 2)} 
prop_de = 0.1
de_opts <- DESetup(ngenes=1000,
 nsims=25,
 p.DE = prop_de,
 pLFC= log2fc_func,
 sim.seed = 58673
 )

寻找未注释的转录区域

一个常见的挑战是找到并计算对齐在注释区域之外的读取值。在 RNAseq 实验中，这些读取值可能代表未注释的基因和新型转录本。实际上，我们有一些已知的基因，并且可以看到它们被转录，因为它们有对齐的读取覆盖，但其他转录区域并未包含在任何注释中，我们希望知道这些读取代表的对齐位置。在这个配方中，我们将探讨一种看似简单的技术来寻找这些区域。

准备就绪

我们的数据集将是一个合成数据集，包含一个小的 6,000 bp 基因组区域和两个具有读取值的基因特征，以及一个未注释的区域，具有对齐的读取值，如下图所示：

我们将需要 Bioconductor 的csaw、IRanges、SummarizedExperiment和rtracklayer库，以及其他一些来自 Bioconductor 基本包的函数。数据位于本书的数据仓库中的datasets/ch1/windows.bam和datasets/ch1/genes.gff。

如何操作...

使用powsimR进行功效分析可以按照以下步骤进行：

设置加载函数：

get_annotated_regions_from_gff <- function(file_name) { 
    gff <- rtracklayer::import.gff(file_name) 
    as(gff, "GRanges") 
}

获取整个基因组的窗口计数：

whole_genome <- csaw::windowCounts( 
    file.path(getwd(), "datasets", "ch1", "windows.bam"),
    bin = TRUE,
    filter = 0,
    width = 500,
    param = csaw::readParam(
        minq = 20,
        dedup = TRUE,
        pe = "both"
    )
)
colnames(whole_genome) <- c("small_data")

annotated_regions <- get_annotated_regions_from_gff(file.path(getwd(), "datasets", "ch1", "genes.gff"))

查找注释与窗口之间的重叠，并对子集进行操作：

library(IRanges)
library(SummarizedExperiment)
windows_in_genes <-IRanges::overlapsAny( SummarizedExperiment::rowRanges(whole_genome), annotated_regions )

将窗口子集分为注释区域和非注释区域：

annotated_window_counts <- whole_genome[windows_in_genes,] 
non_annotated_window_counts <- whole_genome[ ! windows_in_genes,]

将数据提取到计数矩阵中：

assay(non_annotated_window_counts)

它是如何工作的...

在步骤 1中，我们创建一个函数，加载 GFF 文件中的基因区域信息（请参见本书的附录以了解 GFF 的描述），并使用rtracklayer包将其转换为 Bioconductor 的GRanges对象。此方法有效，因为GRanges对象可以像普通的 R 矩阵或数据框一样进行子集选择。它们在这方面表现得像一个“矩阵”，尽管GRanges比矩阵复杂得多，但它的行为却非常相似。这使得一些操作和提取变得非常简单。我们在本食谱中广泛使用GRanges以及相关的类RangedSummarizedExperiment。

在步骤 2中，我们使用csaw windowCounts()函数在 500 bp 的窗口中获取整个基因组的计数。width参数定义了窗口大小，param参数决定了什么是通过的读数；在这里，我们将最小读数质量（minq）设置为 PHRED 评分 20，去除 PCR 重复（dedup = TRUE），并要求每个读数的两个配对都对齐（pe="both"）。返回的whole_genome对象是RangedSummarizedExperiment。我们将whole_genome中单一数据列的名称设置为small_data。最后，我们使用自定义函数get_annotated_regions_from_gff()，根据 GFF 文件中的基因信息创建GRanges对象annotated_regions。

在步骤 3中，我们使用IRanges overlapsAny()函数检查窗口位置是否与基因区域有重叠。此函数需要GRanges对象，因此我们使用SummarizedExperiment的rowRanges()函数从whole_genome变量中提取，并将该对象与现有的GRanges对象的annotated_regions一起传递给overlapsAny()。此函数返回一个二进制向量，我们可以用它来进行子集选择。

在步骤 4中，我们只需使用二进制向量windows_in_genes来对子集whole_genome对象，从而提取注释窗口（存入annotated_window_counts）作为GRanges对象。然后，我们可以使用相同的代码，通过逻辑反转二进制向量（使用!操作符）来获取非注释窗口，从而得到non_annotated_window_counts。

最后，在步骤 5中，我们可以使用assay()函数从GRanges对象中提取实际计数。

还有更多...

我们可能需要从 GFF 以外的其他文件格式获取注释区域。rtracklayer支持多种格式——这里是处理 BED 文件的函数：

get_annotated_regions_from_bed <- function(file_name){ 
    bed <- rtracklayer::import.bed(file_name) 
    as(bed, "GRanges") 
}

使用 bumphunter 方法发现高表达区域的初步步骤

查找所有来自同一、可能未注释的基因组特征的读数比对区域是一项常见任务。这里的目标是将读数比对区域分组，以便标记出具有显著覆盖度的区域，并随后对样本之间的表达水平差异进行比较。

准备工作...

我们将使用相同的windows数据集，其中有一个实验包含三个峰值，并将其传入我们在食谱 4中使用的函数——这样我们就知道我们要查找三个突起。数据存储在本书的数据仓库中的datasets/ch1/windows.bam文件下。我们需要Rsamtools和bumphunter库。

如何进行操作...

加载数据并获取每位置的覆盖度：

library(Rsamtools) 
library(bumphunter) 
pileup_df <- Rsamtools::pileup(file.path(getwd(), "datasets", "ch1", "windows.bam"))

查找初步簇：

clusters <- bumphunter::clusterMaker(pileup_df$seqnames, pileup_df$pos, maxGap = 100)

查找具有最小阈值的突起：

bumphunter::regionFinder(pileup_df$count, pileup_df$seqnames, pileup_df$pos, clusters, cutoff=1)

它是如何工作的...

在步骤 1中，我们使用 Rsamtools 的pileup()函数并使用默认设置，获取每碱基的覆盖度数据框。每一行代表参考中的单个核苷酸，计数列给出该点的覆盖深度。结果存储在pileup_df数据框中。

在步骤 2中，我们使用bumphunter的clusterMaker()函数在pileup_df上进行操作，简单地将距离较近的读数分组到一起形成簇。我们传入序列名称、位置和最大距离参数（maxGap）。该函数返回一个与数据框长度相等的簇编号向量，指示每一行在数据框中的簇归属。如果我们用table进行统计，就能看到每个簇中包含的行数（簇大小）：

table(clusters)
## clusters
##    1    2    3 
## 1486 1552 1520

在步骤 3中，我们改进了方法；我们使用regionFinder()，它应用了一个读深度阈值，以确保簇的最低读深度。我们传入与步骤 2类似的数据，添加簇归属向量 clusters 以及最小读数阈值——在这个非常小的数据集中，我们设置为 1。步骤 3的结果是被聚集在一起的区域，并以有用的表格形式呈现：

##    chr start  end     value  area cluster indexStart indexEnd    L
## 3 Chr1  4503 5500 10.401974 15811       3       3039     4558 1520
## 1 Chr1   502 1500  9.985868 14839       1          1     1486 1486
## 2 Chr1  2501 3500  8.657216 13436       2       1487     3038 1552

在这些区域预测中，我们可以清楚地看到这三个位于数据中的区域，读数大致会有几个核苷酸的偏差。

还有更多...

如果你有多个实验需要分析，可以尝试bumphunter()函数。它会在矩阵中的多个数据列上操作，并进行线性建模，以评估来自重复样本的位点和存在的不确定性；它的操作方式与regionFinder()非常相似。

差异峰值分析

当你发现未注释的转录本时，可能想要看看它们在不同实验之间是否有差异表达。我们已经了解了如何使用edgeR和DESeq来进行这一分析，但一个问题是如何将一个如RangedSummarizedExperiment的对象（它由数据和描述峰值区域的GRanges对象组成）转换为内部的DESeq对象。在这个食谱中，我们将探讨如何将这些对象中的数据汇总并将其转换为正确的格式。

准备工作

对于这个食谱，你需要本书仓库中的datasets/ch1/arabidopsis_rse.RDS文件，这里面包含了拟南芥 RNAseq 的RangedSummarizedExperiment版本。我们还将继续使用之前使用过的DESeq和SummarizedExperiment Bioconductor 包。

如何做...

加载数据并设置一个创建区域标签的函数：

library(SummarizedExperiment) 
arab_rse <- readRDS(file.path(getwd(), "datasets", "ch1", "arabidopsis_rse.RDS") ) 

 make_tag <- function(grange_obj){
     paste0( 
        grange_obj@seqnames, 
        ":", 
        grange_obj@ranges@start, 
        "-", 
        (grange_obj@ranges@start + grange_obj@ranges@width) 
    ) 
}

提取数据并注释行：

counts <- assay(arab_rse)

if ( ! is.null(names(rowRanges(arab_rse))) ){
  rownames(counts) <- names(rowRanges(arab_rse))
} else {
  rownames(counts) <- make_tag(rowRanges(arab_rse))
}

它是如何工作的...

步骤 1开始时加载我们预先准备好的RangedSummarized实验数据；请注意，GRanges对象中的names槽没有填充。接下来，我们创建一个自定义函数make_tag()，它通过拼接seqnames、starts和通过传递的GRanges对象计算出的结束位置（start + width）来工作。请注意@符号语法：这是因为GRange是一个 S4 对象，槽是通过@而不是熟悉的$来访问的。

在步骤 2中，代码通过assay()函数从RangedSummarizedExperiment中提取实际数据。返回的矩阵没有行名，这样不太有用，因此我们使用if语句来检查names槽——如果它可用，我们将其作为行名；如果不可用，我们使用GRanges对象中的位置信息，在我们创建的make_tag()函数中生成行名标签。这样就会得到以下输出——一个计数矩阵，其中位置标签作为行名，可以在本章的食谱 1和食谱 2中描述的 DESeq 和 edgeR 中使用：

head(counts)
##                  mock1 mock2 mock3 hrcc1 hrcc2 hrcc3
## Chr1:3631-5900      35    77    40    46    64    60
## Chr1:5928-8738      43    45    32    43    39    49
## Chr1:11649-13715    16    24    26    27    35    20
## Chr1:23146-31228    72    43    64    66    25    90
## Chr1:31170-33154    49    78    90    67    45    60
## Chr1:33379-37872     0    15     2     0    21     8

使用 SVA 估计批次效应

高通量数据，如 RNAseq，常常受到实验设计中未明确建模的技术误差的影响，这些误差可能会混淆差异表达的检测。来自不同日期、不同地点或不同机器上运行的样本的计数差异就是一个常见的技术误差示例，这种误差通常存在，我们应该尽量在实验设计中进行建模。解决这个问题的一种方法是将代理变量纳入我们的实验设计中，解释批次效应，并在建模和差异表达分析阶段考虑它。我们将使用SVA包来实现这一点。

准备工作

我们需要 Bioconductor 中的 SVA 包和我们的拟南芥计数数据，位于datasets/ch1/arabidopsis.RDS。

如何做...

使用 SVA 估计批次效应，请执行以下步骤：

加载库和数据：

library(sva)
arab <- readRDS(file.path(getwd(), "datasets", "ch1", "arabidopsis.RDS"))

在某些实验中过滤掉计数过少的行：

keep <- apply(arab, 1, function(x) { length(x[x>3])>=2 } )
arab_filtered <- arab[keep,]

创建初步设计：

groups <- as.factor(rep(c("mock", "hrcc"), each=3))

设置测试模型和空模型并运行 SVA：

test_model <- model.matrix(~groups)
null_model <- test_model[,1]
svar <- svaseq(arab_filtered, test_model, null_model, n.sv=1)

提取代理变量以用于后续的设计：

design <- cbind(test_model, svar$sv)

它是如何工作的...

在步骤 1中，脚本首先加载拟南芥 RNAseq 数据的计数矩阵，你会记得这是一个简单的实验，包含三次mock处理和三次hrcc处理的重复。

在步骤 2中，我们创建了一个包含我们希望保留的行索引的向量，这通过测试每一行是否至少有两列计数大于 3 来完成——这是通过使用apply()和匿名函数在计数矩阵的行上进行操作实现的。

在步骤 3中，我们创建了一个groups因子，描述了实验样本的类型。

步骤 4是实际操作的部分；我们在model.matrix()中使用groups因子来创建我们想要使用的模型设计。我们还需要一个空模型，在此实验设计中，它等同于第一列，因此我们从test_model设计对象中提取它。然后，我们可以使用关键的svaseq()函数来估计需要添加到设计中的代理变量。我们将test_model和null_model结合，并使用n.sv来选择使用的代理变量数量，对于像这样的简单设计，通常设置为 1。

最后一步，步骤 5，是将代理变量添加到设计模型中，我们通过将test_model和svar中的sv列（svsar$sv）结合在一起完成此操作。最终的设计对象可以像其他任何设计对象一样用于edgeR和DESeq2等包，且这些方法将使用代理变量来处理批次效应。

使用 AllelicImbalance 寻找等位基因特异性表达

等位基因特异性表达是指在转录本的不同等位基因变异之间存在差异丰度的情况。RNAseq 有助于为具有转录变异的基因提供等位基因特异性表达的定量估计——即转录本中可能导致不同蛋白质的变异。在本教程中，我们将探讨一种方法，用于确定转录本的哪些变异可能在不同样本中表现出偏向性表达。这些读取数据将来自不同的.bam文件，且变异必须已经被识别。这意味着你已经完成了读取比对和变异调用步骤，并且拥有每个样本的.bam和.vcf文件。我们将在此教程中使用AllelicImbalance和VariantAnnotation包。

准备工作

你需要从 Bioconductor 中安装AllelicImbalance和VariantAnnotation。AllelicImbalance包提供了一个小但信息丰富的数据集，包含人类基因组 hg19 版本上 17 号染色体的三个 SNP。文件已被提取到本书的数据仓库中的datasets/ch1/allele_expression目录下。

如何操作...

加载库并设置导入文件夹：

library(AllelicImbalance)
library(VariantAnnotation)

region_of_interest <- GRanges(seqnames = c("17"), ranges = IRanges(79478301, 79478361)) 
bam_folder <- file.path(getwd(), "datasets", "ch1", "allele_expression")

加载目标区域的读取数据和变异信息：

reads <- impBamGAL(bam_folder, region_of_interest, verbose = FALSE)

vcf_file <-file.path( getwd(), "datasets", "ch1", "allele_expression","ERP000101.vcf" )
variant_positions <- granges(VariantAnnotation::readVcf(vcf_file), "hg19" )

allele_counts <- getAlleleCounts(reads, variant_positions, verbose=FALSE)

构建 ASEset 对象：

ase.vcf <- ASEsetFromCountList(rowRanges = variant_positions, allele_counts)

reference_sequence <- file.path(getwd(), "datasets", "ch1", "allele_expression", "hg19.chr17.subset.fa")

ref(ase.vcf) <- refAllele(ase.vcf,fasta=reference_sequence)
alt(ase.vcf) <- inferAltAllele(ase.vcf)

对所有变异进行测试：

binom.test(ase.vcf, n="*")

它是如何工作的...

在步骤 1中，脚本首先创建了一个熟悉的GRanges对象，描述了我们感兴趣的区域以及包含.bam文件的文件夹。

在步骤 2中，impBamGAL()函数加载了目标区域的读取数据。变异信息被加载到variant_positions中——另一个GRanges对象，读取数据和变异用于使用getAlleleCounts()生成等位基因计数。*

完成这一部分后，在第三步，我们可以构建ASESet对象ase.vcf（一个继承自RangedSummarizedExperiment的类），使用构造函数ASEsetFromCountList()；然后我们使用设置函数ref()和alt()，应用参考和替代碱基身份。

最后，在第四步，我们可以应用测试。binom.test()对每个位置每个样本（.bam文件）进行二项分布检验，测试参考和替代等位基因的计数是否有偏差。n参数告诉测试考虑哪个链——在这个例子中，我们没有设置每条链的信息，因此使用"*"来忽略链特性。

这将给出以下输出：

##               chr17_79478287 chr17_79478331 chr17_79478334
## ERR009113.bam          0.500   1.000000e+00   1.000000e+00
## ERR009115.bam         0.125   6.103516e-05   3.051758e-05

还有更多...

如果需要，前面的分析可以扩展为每条链和每种表型的测试。脚本需要进行修改，以在ASESet对象构造步骤中引入链信息。通常，这需要 RNAseq 实验和对齐步骤在考虑链特性的前提下执行，并且生物信息学流程需按此配置。可以在构造步骤中使用colData参数和一个表型或样本类型的向量来添加表型信息，这些信息对应ASESet对象中的列。

绘制和展示 RNAseq 数据

绘制 RNAseq 数据总览或单个基因或特征的图表是质量控制和理解的重要步骤。在这个步骤中，我们将看到如何绘制样本中基因计数图、如何创建一个 MA 图（它绘制计数与倍数变化，并帮助我们发现与表达相关的样本偏差），以及如何创建一个火山图（它绘制显著性与倍数变化，并帮助我们轻松发现最有意义的变化）。

准备工作

在这个步骤中，我们将使用DESeq2包、ggplot2包、magrittr和dplyr。我们将使用为modencodefly数据创建的DESeqDataSet对象，第 2 步配方中保存的版本位于本书数据仓库中的datasets/ch1/modencode_dds.RDS文件中。

如何做到这一点...

加载库并创建 RNAseq 结果的数据框：

library(DESeq2)
library(magrittr)
library(ggplot2)

dds <- readRDS("~/Desktop/r_book/datasets/ch1/modencode_dds.RDS")

创建一个单一基因的箱线图，基于"stage"进行条件化：

plotCounts(dds, gene="FBgn0000014", intgroup = "stage", returnData = TRUE) %>%
  ggplot() + aes(stage, count) + geom_boxplot(aes(fill=stage)) + scale_y_log10() + theme_bw()

创建一个基于显著性条件的 MA 图：

result_df <- results(dds, contrast=c("stage","L2Larvae","L1Larvae"), tidy= TRUE) %>%
  dplyr::mutate(is_significant=padj<0.05)

ggplot(result_df) + aes(baseMean, log2FoldChange) + geom_point(aes(colour=is_significant)) + scale_x_log10() + theme_bw()

创建一个基于显著性条件的火山图：

ggplot(result_df) + aes(log2FoldChange, -1 * log10(pvalue))  + geom_point(aes(colour=is_significant)) + theme_bw()

它是如何工作的...

第一步是简要加载我们需要的数据集和库。

在第二步，我们利用了DESeq2中plotCounts()和results()函数中的一些有用参数。plotCounts()中的returnData标志将可选地返回给定条件下某个基因的计数信息的整洁数据框，因此我们可以将数据传递给ggplot()来制作单个基因的箱线图。magrittr的%>%操作符允许我们将plotCounts()的返回值直接传递给ggplot()的第一个位置参数，而无需先保存到中间变量中。

在步骤 3中，我们使用 DESeq2 中的results()函数来获取results数据框，并将其传递给dplyr的mutate()函数，以添加一个名为is_significant的新列，如果padj列的值小于 0.05，则该列值为TRUE。然后，我们在ggplot()命令中使用返回的result_df数据框，绘制baseMean（计数）与 log2 折叠变化之间的散点图，点的颜色由is_significant变量决定，实际上是根据 P 值是否小于 0.05 来着色。

在步骤 4中，我们使用相同的result_df数据框，绘制 log2 折叠变化与'pvalue'的负 log10 值的关系，得到一个'火山图'，展示 P 值与差异表达水平之间的关系：

前面三个图是这三个`ggplot``()命令的合成结果输出。

第二章：使用 HTS 数据查找遗传变异

高通量测序（HTS）使得在短时间内发现遗传变异并进行全基因组基因分型和单倍型分析成为可能。这项技术所释放的数据洪流为生物信息学家和计算机科学家提供了一些独特的机会，创造了很多创新的新数据存储和分析流程。变异调用的基本流程从 HTS 读取的质量控制开始，然后是将这些读取比对到参考基因组。这些步骤通常会在 R 分析之前进行，并通常会生成一个包含读取比对的 BAM 文件或包含变异位置的 VCF 文件（有关这些文件格式的简要讨论，请参见本书的附录），我们将在 R 代码中对其进行处理。

由于变异调用和分析是生物信息学中的基本技术，Bioconductor 配备了构建软件和进行分析所需的工具。研究人员想要问的关键问题可能从 我的基因组上遗传变异的位置在哪里？ 到 有多少个变异？ 再到 如何对它们进行分类？ 我们将查看一些解决这些问题的方案，并且还将探讨一些重要的通用技术，这些技术使我们能够在基因组上可视化变异和标记，并评估变异与基因型的关联。我们还将了解遗传变异的其他定义，并探索如何评估单个位点的拷贝数。

本章将涵盖以下方案：

使用 VariantTools 查找序列数据中的 SNP 和 indel
预测长参考序列中的开放阅读框
使用 karyoploteR 在遗传图谱上绘制特征
寻找替代转录本同种型
使用 VariantAnnotation 选择和分类变异
提取感兴趣基因组区域的信息
查找表型和基因型与 GWAS 的关联
估计感兴趣位点的拷贝数

技术要求

以下是你将需要的 R 包。一些可以通过 install.packages() 安装。列在 Bioconductor 下的包需要通过专用安装器安装。相关说明请参见此处。如果你需要进行其他操作，包的安装将在使用这些包的方案中进行描述：

Bioconductor：以下是这些包：
- Biostrings
- GenomicRanges
- gmapR
- karyoploteR
- rtracklayer
- systemPipeR
- SummarizedExperiment
- VariantAnnotation
- VariantTools
rrBLUP

Bioconductor 非常庞大，并且拥有自己的安装管理器。你可以通过以下代码安装这些包（更多信息请参考 www.bioconductor.org/install/）：

if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install()

通常，在 R 中，用户会加载一个库并直接通过名称使用其中的函数。这在交互式会话中非常方便，但当加载多个包时，可能会导致混淆。为了明确当前使用的是哪个包和函数，我有时会使用packageName::functionName()的约定。

有时，在一个过程的中途，我会暂停代码，以便你可以看到一些中间输出或理解某个对象的结构。每当发生这种情况时，你会看到一个代码块，每行的开头都会有双井号（##）符号，如下所示：

letters[1:5]
## a b c d e

本章所有代码和数据都可以在本书的 GitHub 仓库找到：github.com/danmaclean/R_Bioinformatics_Cookbook。

使用 VariantTools 从序列数据中找到 SNP 和插入缺失变异（indels）

一个关键的生物信息学任务是处理高通量序列读取的比对结果，通常存储在 BAM 文件中，并计算变异位置列表。当然，这个过程可以通过许多外部命令行程序和工具来完成，通常会生成一个 VCF 格式的变异文件，但 Bioconductor 中有一些非常强大的包能够完成整个流程，并且以快速高效的方式，通过利用 BiocParallel 的并行计算功能——这一功能旨在加速处理 Bioconductor 对象中的大规模数据集——实现这一目标。使用 Bioconductor 工具可以让我们将所有处理步骤都保持在 R 环境中，在这一部分中，我们将使用纯粹的 R 代码和一些 Bioconductor 包，逐步完成一个从读取数据到变异基因列表的完整流程。

准备工作

在本节中，我们将使用一组合成的读取数据，覆盖人类基因组 17 号染色体前 83KB 左右的区域。这些读取数据是使用samtools中的wgsim工具生成的——一个外部命令行程序。它们包含 64 个由wgsim引入的 SNP，这些 SNP 可以在datasets/ch2/snp_positions.txt中的样本数据中看到。你会看到，随着程序的进展，默认参数会找到比实际更多的 SNP——你需要仔细调节参数，以精细调整 SNP 发现过程。

如何操作...

使用VariantTools从序列数据中找到 SNP 和插入缺失变异（indels）可以通过以下步骤完成：

导入所需的库：

library(GenomicRanges)
library(gmapR)
library(rtracklayer)
library(VariantAnnotation)
library(VariantTools)

然后，加载数据集：

bam_folder <- file.path(getwd(), "datasets", "ch2")
bam_folder_contents <- list.files(file.path(getwd(), "datasets", "ch2" ) )
bam <- file.path( bam_folder, "hg17_snps.bam")
fasta_file <- file.path(bam_folder,"chr17.83k.fa")

设置基因组对象和参数对象：

fa <- rtracklayer::FastaFile(fasta_file)

genome <- gmapR::GmapGenome(fa, create=TRUE)

qual_params <- TallyVariantsParam(
                   genome = genome,
                   minimum_mapq = 20)

var_params <- VariantCallingFilters(read.count = 19,
                                    p.lower = 0.01
                                    )

调用变异：

called_variants <- callVariants(bam, qual_params, 
                                calling.filters = var_params
                                )

head(called_variants)

现在，我们进入注释阶段并加载来自.gff或.bed文件的特征位置信息：

get_annotated_regions_from_gff <- function(file_name) {
  gff <- rtracklayer::import.gff(file_name) 
  as(gff, "GRanges")
}

get_annotated_regions_from_bed <- function(file_name){
  bed <- rtracklayer::import.bed(file_name)
  as(bed, "GRanges")
}

genes <- get_annotated_regions_from_gff(file.path( bam_folder, "chr17.83k.gff3"))

现在我们计算哪些变异与哪些基因重叠：

overlaps <- GenomicRanges::findOverlaps(called_variants, genes) 
overlaps

最后，我们通过重叠列表来筛选基因。

genes[subjectHits(overlaps)]

它是如何工作的...

这是一个复杂且涉及多个步骤的管道。加载库后，前四行设置了我们需要的来自数据集目录的文件。请注意，我们需要一个.bam 文件和一个 fasta 文件。接下来，我们使用 gmapR::GmapGenome() 函数与 fasta 对象一起创建一个 GmapGenome 对象——这个对象描述了基因组以供后续的变异调用函数使用。接下来的两个函数，TallyVariantParams() 和 VariantCallingFilters()，对于正确调用和筛选候选 SNP 至关重要。这些函数是你可以设置定义 SNP 或 indel 参数的地方。这里的选项故意设置得很少。正如从输出中看到的，虽然我们创建了 64 个 SNP，但只有 6 个被成功调用。

一旦参数定义完毕，我们使用 callVariants() 函数，将我们设置的所有信息作为输入，获得一个变异的 vranges 对象。

这将产生以下输出：

 VRanges object with 6 ranges and 17 metadata columns:
##           seqnames    ranges strand         ref              alt
##              <Rle> <IRanges>  <Rle> <character> <characterOrRle>
##   [1] NC_000017.10        64      *           G                T
##   [2] NC_000017.10        69      *           G                T
##   [3] NC_000017.10        70      *           G                T
##   [4] NC_000017.10        73      *           T                A
##   [5] NC_000017.10        77      *           T                A
##   [6] NC_000017.10        78      *           G                T

然后，我们可以设置 GFF 文件注释的 GRanges 对象（我还提供了一个从 BED 文件获取注释的函数）。

这将产生以下输出：

## Hits object with 12684 hits and 0 metadata columns:
##           queryHits subjectHits
##           <integer>   <integer>
##       [1]     35176           1
##       [2]     35176           2
##       [3]     35176           3
##       [4]     35177           1

最后的步骤是使用 XRanges 对象强大的重叠和子集功能。我们使用 GenomicRanges::findOverlaps() 查找实际的重叠——返回的 overlaps 对象实际上包含了每个输入对象中重叠对象的索引。

这将产生以下输出：

## GRanges object with 12684 ranges and 20 metadata columns:
##               seqnames      ranges strand |   source       type     score
##                  <Rle>   <IRanges>  <Rle> | <factor>   <factor> <numeric>
##       [1] NC_000017.10 64099-76866      - |   havana ncRNA_gene      <NA>
##       [2] NC_000017.10 64099-76866      - |   havana    lnc_RNA      <NA>
##       [3] NC_000017.10 64099-65736      - |   havana       exon      <NA>

因此，我们可以使用 subjectHits(overlaps) 直接筛选出包含 SNP 的基因，并获得一个非常简洁的列表。

还有更多内容...

当我们对筛选条件和我们调用的变异集感到满意时，可以使用以下代码保存变异的 VCF 文件：

VariantAnnotation::sampleNames(called_variants) <- "sample_name"
vcf <- VariantAnnotation::asVCF(called_variants)
VariantAnnotation::writeVcf(vcf, "hg17.vcf")

另见

虽然我们的方案使步骤和代码清晰易懂，但我们需要更改的实际参数和数值不能像这样直接描述，因为这些值将极大依赖于数据集。VariantTools 文档包含了如何正确设置参数的详细讨论：bioconductor.org/packages/release/bioc/vignettes/VariantTools/inst/doc/VariantTools.pdf。

在长参考序列中预测开放阅读框

先前未测序基因组的草图组装可以是一个丰富的生物学知识来源，但当基因组学资源如基因注释不可用时，继续进行可能会变得棘手。在这里，我们将看看一个第一阶段的流程，用来寻找潜在的基因和基因组位点，完全是de novo的，且没有超出序列的信息。我们将使用一组非常简单的规则来寻找开放阅读框——这些序列从起始密码子开始，并以终止密码子结束。实现这一点的工具被封装在Bioconductor包中的一个函数systemPipeR内。我们最终会得到另一个GRanges对象，可以将其整合到下游的处理流程中，帮助我们交叉引用其他数据，例如 RNAseq，正如我们在第一章《执行定量 RNAseq》中的查找未注释的转录区食谱中看到的那样。最后一步，我们将看看如何使用基因组模拟来评估哪些开放阅读框实际上可能是有效的，而不是仅仅偶然出现的。

做好准备

在这个食谱中，我们只需要输入Arabidopsis 叶绿体基因组的短 DNA 序列；它位于datasets/ch2/arabidopsis_chloroplast.fa.我们还需要Bioconductor包中的Biostrings和systemPipeR。

如何操作...

在长参考序列中预测开放阅读框可以通过以下步骤完成：

加载库并输入基因组：

library(Biostrings)
library(systemPipeR)

dna_object <- readDNAStringSet(file.path(getwd(), "datasets","ch2", "arabidopsis_chloroplast.fa"))

预测ORFs（开放阅读框）：

predicted_orfs <- predORF(dna_object, n = 'all', type = 'gr', mode='ORF', strand = 'both', longest_disjoint = TRUE) 
predicted_orfs

计算参考基因组的特性：

bases <- c("A", "C", "T", "G")
raw_seq_string <- strsplit(as.character(dna_object), "")

seq_length <- width(dna_object[1])
counts <- lapply(bases, function(x) {sum(grepl(x, raw_seq_string))}  )
probs <- unlist(lapply(counts, function(base_count){signif(base_count / seq_length, 2) }))

创建一个函数，找到模拟基因组中最长的 ORF：

get_longest_orf_in_random_genome <- function(x,
  length = 1000, 
  probs = c(0.25, 0.25, 0.25, 0.25), 
  bases = c("A","C","T","G")){

  random_genome <- paste0(sample(bases, size = length, replace = TRUE, prob = probs), collapse = "")
  random_dna_object <- DNAStringSet(random_genome)
  names(random_dna_object) <- c("random_dna_string")
  orfs <- predORF(random_dna_object, n = 1, type = 'gr', mode='ORF', strand = 'both', longest_disjoint = TRUE)
  return(max(width(orfs)))
}

在 10 个模拟基因组上运行该函数：

random_lengths <- unlist(lapply(1:10, get_longest_orf_in_random_genome, length = seq_length, probs = probs, bases = bases))

获取最长随机 ORF 的长度：

longest_random_orf <- max(random_lengths)

只保留预测的长度大于最长随机 ORF 的 ORF：

keep <- width(predicted_orfs) > longest_random_orf
orfs_to_keep <- predicted_orfs[keep]
orfs_to_keep

它是如何工作的...

这个食谱的第一部分是我们实际预测 ORF 的地方。首先，我们使用Biostrings包中的readDNAStringSet()函数将 DNA 序列加载为DNAStringSet对象。systemPipeR包中的predORF()函数使用这个对象作为输入，并根据设置的选项实际预测开放阅读框。在这里，我们返回两个链上的所有 ORF。

这将产生以下输出：

## GRanges object with 2501 ranges and 2 metadata columns:
##           seqnames        ranges strand | subject_id inframe2end
##              <Rle>     <IRanges>  <Rle> |  <integer>   <numeric>
##      1 chloroplast   86762-93358      + |          1           2
##   1162 chloroplast     2056-2532      - |          1           3
##      2 chloroplast   72371-73897      + |          2           2
##   1163 chloroplast   77901-78362      - |          2           1
##      3 chloroplast   54937-56397      + |          3           3

我们返回一个GRanges对象，其中描述了 2,501 个开放阅读框。这个数字太多了，所以我们需要过滤掉其中一些；特别是，我们可以找出哪些 ORF 是由序列中的偶然事件引起的。为此，我们需要做一点模拟，这也正是接下来代码部分要做的事情。

为了估计随机 ORF 的最大长度，我们将创建一系列长度与输入序列相等且具有相同碱基比例的随机基因组，并观察可以预测的最长 ORF 是什么。我们进行几次迭代，并大致了解由偶然事件引起的最长 ORF 的长度。这个长度将作为一个阈值，用于拒绝真实序列中预测的 ORF。

实现这一目标需要一些设置和自定义函数。首先，我们定义将作为简单字符向量使用的碱基。然后，我们通过分割dna_object的as.character版本来获取原始 DNA 序列的字符向量。我们利用这些信息计算输入序列中每个碱基的比例，首先计算每个碱基的数量（得到counts），然后将其除以序列的长度，得到probs。在这两个步骤中，我们使用lapply()遍历bases向量和counts列表，应用匿名函数来计算这些变量的结果列表。unlist()用于将最终列表简化为一个简单的向量。

一旦我们完成设置，就可以构建我们的get_longest_orf_in_random_genome()模拟函数。该函数通过从bases中的选择中按给定的probs概率采样字符生成一个随机基因组。该向量通过paste0()合并为一个字符串，然后转换为DNAStringSet对象，以供predORF()函数使用。这一次，我们只请求最长的 ORF，使用n = 1并返回其长度。

这将产生以下输出：

## GRanges object with 10 ranges and 2 metadata columns:
##         seqnames        ranges strand | subject_id inframe2end
##            <Rle>     <IRanges>  <Rle> |  <integer>   <numeric>
##    1 chloroplast   86762-93358      + |          1           2
##    2 chloroplast   72371-73897      + |          2           2
##    3 chloroplast   54937-56397      + |          3           3
##    4 chloroplast   57147-58541      + |          4           1

现在，我们可以运行这个函数，我们使用lapply()运行 10 次，并使用我们之前计算的length、probs和bases信息。unlist()将结果转换为一个简单的向量，我们使用max()提取 10 次运行中的最长一个。我们可以通过子集化原始的predicted_orfs GRanges对象，保留那些比随机生成的 ORFs 更长的 ORF。

还有更多...

一旦你获得了一组你满意的 ORF，你可能希望将它们保存到文件中。你可以使用BSgenome包中的getSeq()函数，通过传递原始序列对象——dna_object——以及orfs_to_keep中的范围来完成。然后，使用names()为结果命名，最后，你可以使用writeXStringSet()函数将它们保存到文件中：

extracted_orfs <- BSgenome::getSeq(dna_object, orfs_to_keep) 
names(extracted_orfs) <- paste0("orf_", 1:length(orfs_to_keep))
writeXStringSet(extracted_orfs, "saved_orfs.fa")

使用 karyoploteR 绘制遗传图谱上的特征

我们能做的最有价值且富有洞察力的事情之一就是可视化数据。我们常常希望在染色体或遗传图谱上看到某些感兴趣的特征与其他特征的关系。这些图有时称为染色体图，有时称为示意图，在本节中，我们将看到如何使用karyoploteR包创建这些图。该包将熟悉的GRanges对象作为输入，并根据配置生成详细的图表。我们将快速浏览几种不同的图表样式，以及一些配置选项，以解决图表中的问题，比如标签溢出或重叠。

准备工作

对于这个配方，你需要安装karyoploteR，但我们将使用的所有数据将在配方本身中生成。

如何进行操作...

使用karyoploteR绘制遗传图谱上的特征，可以通过以下步骤完成：

首先，我们加载库：

library(karyoploteR)
library(GenomicRanges)

然后，设置基因组对象，这将作为我们核型图的基础：

genome_df <- data.frame( 
    chr = paste0("chr", 1:5), 
    start = rep(1, 5), 
    end = c(34964571, 22037565, 25499034, 20862711, 31270811) 
    ) 
genome_gr <- makeGRangesFromDataFrame(genome_df)

设置我们将绘制为标记的 SNP 位置：

snp_pos <- sample(1:1e7, 25)
snps <- data.frame(
  chr = paste0("chr", sample(1:5,25, replace=TRUE)),
  start = snp_pos,
  end = snp_pos
)
snps_gr <- makeGRangesFromDataFrame(snps)

为标记创建一些标签：

snp_labels <- paste0("snp_", 1:25)

设置绘图边距：

plot.params <- getDefaultPlotParams(plot.type=1)
plot.params$data1outmargin <- 600

创建基础图表并添加轨道：

kp <- plotKaryotype(genome=genome_gr, plot.type = 1, plot.params = plot.params)
kpPlotMarkers(kp, snps_gr, labels = snp_labels)

它是如何工作的...

代码首先加载我们需要的库，然后我们构建一个data.frame，描述我们要绘制的基因组，并相应地设置名称和长度。然后，使用makeGRangesFromDataFrame()转换函数将data.frame转换为genome_gr——一个GRanges对象。接着，我们使用sample()函数生成一个包含 25 个随机 SNP 的data.frame，并选择一个位置和染色体。同样，这些数据也被转换为GRanges。现在我们可以设置我们的绘图。首先，我们通过getDefaultPlotParams()从包内部获取默认的绘图参数对象。我们可以修改此对象，以便对绘图的默认设置进行更改。

注意，我们选择了plot.type = 1——这是一个简单的绘图，具有直接位于每个染色体区域上方的数据轨道。我们需要调整数据轨道的边距高度，以防止标记标签溢出到顶部——这可以通过plot.params$data1outmargin <- 600来完成。最后，我们可以绘制我们的图表；我们通过调用plotKaryotype()并传入genome_gr对象、plot.type以及在修改后的plot.params对象中的参数来创建基本的绘图对象kp。

这将产生以下输出：

我们的标记是使用kpPlotMarkers()函数绘制的，传入新的kp绘图对象、snps_gr数据和 SNP 标签。

还有更多...

我们可以通过karyoploteR将多种类型的数字数据添加到数据轨道中。以下示例演示了如何将一些数字数据绘制成简单的线条。第一步是准备我们的数据。在这里，我们创建了一个包含 100 个随机数的data.frame，这些随机数映射到染色体 4 的 100 个窗口，并且像之前一样，我们创建了一个GRanges对象。这次，我们将在染色体的上方和下方分别添加一个数据轨道——一个用于 SNP 标记，另一个用于新数据（注意，这里的plot.type = 2）。接下来，我们需要设置绘图的参数——特别是边距，以防标签和数据重叠；但是之后，绘图的调用与之前相同，这次加入了kpLines()的调用。这里的关键参数是y，它描述了数据在每个绘图点的y值（请注意，这来自我们的numeric_data对象中的单列数据）。现在我们得到了一个包含数字数据轨道的绘图，沿着染色体 4 进行显示。以下是此示例所需执行的步骤：

创建一些数字数据：

numeric_data <- data.frame(
 y = rnorm(100,mean = 1,sd = 0.5 ),
 chr = rep("chr4", 100),
 start = seq(1,20862711, 20862711/100),
 end = seq(1,20862711, 20862711/100)
)
numeric_data_gr <- makeGRangesFromDataFrame(numeric_data)

设置绘图边距：

plot.params <- getDefaultPlotParams(plot.type=2)
plot.params$data1outmargin <- 800
plot.params$data2outmargin <- 800
plot.params$topmargin <- 800

创建一个绘图并添加轨道：

kp <- plotKaryotype(genome=genome_gr, plot.type = 2, plot.params = plot.params)
kpPlotMarkers(kp, snps_gr, labels = snp_labels)
kpLines(kp, numeric_data_gr, y = numeric_data$y, data.panel=2)

这将产生以下输出：

参见：

这里有许多未涵盖的轨道类型和绘图布局。请尝试查看 karyoploteR 手册，以获取完整列表：bioconductor.org/packages/release/bioc/vignettes/karyoploteR/inst/doc/karyoploteR.html。

karyoploteR的一个特点是它仅绘制水平的染色体图。对于垂直图，Bioconductor 中还有chromPlot包。

使用 VariantAnnotation 选择和分类变异

在我们调用变异的管道中，我们通常希望进行后续分析步骤，这些步骤需要根据单个变异的特征（例如，备选等位基因的覆盖深度）进一步进行过滤或分类。这最好从 VCF 文件中进行，常见的做法是保存一个包含所有变异的 VCF 文件，之后再进行过滤。在这一部分中，我们将讨论如何获取输入 VCF 并进行过滤，以保留那些备选等位基因在样本中为主要等位基因的变异。

准备工作

我们需要一个tabix索引的 VCF 文件；我提供了一个在datasets/ch2/sample.vcf.gz文件中。我们还需要 Bioconductor 包VariantAnnotation。

如何操作...

使用VariantAnnotation选择和分类变异可以通过以下步骤完成：

创建一个预过滤函数：

is_not_microsat <- function(x){ !grepl("microsat", x, fixed = TRUE)}

将预过滤函数加载到FilterRules对象中：

prefilters <- FilterRules(list(microsat = is_not_microsat) )

创建一个过滤函数，以保留那些参考等位基因出现在少于一半读取中的变异：

major_alt <- function(x){
 af <- info(x)$AF 
 result <- unlist(lapply(af, function(x){x[1] < 0.5}))
 return(result)
}

将过滤函数加载到一个FilterRules对象中：

filters <- FilterRules(list(alt_is_major = major_alt))

加载输入的 VCF 文件并应用过滤：

vcf_file <- file.path(getwd(), "datasets", "ch2", "sample.vcf.gz")
filterVcf(vcf_file, "hg17", "filtered.vcf", prefilters = prefilters, filters = filters)

它是如何工作的...

在这个非常简短的脚本中，实际上有大量的内容。总体思路是我们需要定义两组过滤规则——prefilter和filter。这通过定义函数来实现，函数接受解析后的 VCF 记录并返回TRUE，如果该记录通过了过滤。预过滤通常是基于未解析的 VCF 记录行的文本过滤——即记录的原始文本。我们的第一行代码定义了一个is_not_microsat()函数，该函数接收一个字符串，并使用grepl()函数判断该行是否包含单词microsat，如果不包含，则返回TRUE。预过滤函数被打包到我们称之为prefilters的FilterRules对象中。

过滤器更为复杂。这些过滤器作用于解析后的 VCF 记录（作为VCF类对象）。我们的major_alt()函数使用info()VCF访问器函数提取 VCF 记录中的info数据。它返回一个数据框，其中每一列都是info部分的一个独立部分。我们提取AF列，它返回一个列表，其中每个 VCF 都有一个元素。为了遍历这些元素，我们使用lapply()函数应用一个匿名函数，如果参考等位基因的比例低于 0.5（即替代等位基因是主要等位基因），则返回TRUE。然后我们使用unlist()将结果转换为向量。major_alt()函数然后被打包到一个名为filters的FilterRules对象中。

最后，完成所有设置后，我们可以加载输入 VCF 文件并使用filterVCF()函数进行过滤。此函数需要FilterRules对象和输出的过滤后的 VCF 文件名。我们使用filtered.vcf作为要写入的文件。

另请参见

在过滤函数中，我们可以利用其他访问器函数来获取 VCF 记录的不同部分。有geno()和fixed()函数，它们将返回描述这些 VCF 记录部分的数据结构。你可以像使用info()一样使用它们来创建过滤器。

提取感兴趣的基因组区域信息

很多时候，你可能需要更详细地查看落在特定基因组区域内的数据，无论是基因中的 SNP 和变异，还是特定基因座中的基因。这个极其常见的任务可以通过功能强大的GRanges和SummarizedExperiment对象来很好地处理，尽管这些对象的设置可能有些繁琐，但它们具有非常灵活的子集操作，付出的努力是值得的。我们将探讨几种设置这些对象的方法以及几种操作它们以获取有趣信息的方法。

准备工作

在这个示例中，我们需要GenomicRanges、SummarizedExperiment和rtracklayer Bioconductor 包。我们还需要两个输入数据文件：Arabidopsis染色体 4 的 GFF 文件，位于datasets/ch2/arabidopsis_chr4.gff文件中，以及同一染色体的基因仅特征的较小文本版本，位于datasets/ch2/arabidopsis_chr4.txt。

如何实现……

提取感兴趣的基因组区域信息可以通过以下步骤完成：

加载包并定义一些从常见文件创建GRanges的函数：

library(GenomicRanges)
library(rtracklayer)
library(SummarizedExperiment)

get_granges_from_gff <- function(file_name) {
  gff <- rtracklayer::import.gff(file_name)
  as(gff, "GRanges")
}

get_granges_from_bed <- function(file_name){
  bed <- rtracklayer::import.bed(file_name)
  as(bed, "GRanges")
}

get_granges_from_text <- function(file_name){
  df <- readr::read_tsv(file_name, col_names = TRUE ) 
  GenomicRanges::makeGRangesFromDataFrame(df, keep.extra.columns = TRUE)
}

实际上使用这些函数创建一些GRanges对象：

gr_from_gff <- get_annotated_regions_from_gff(file.path(getwd(), "datasets", "ch2", "arabidopsis_chr4.gff"))
gr_from_txt <- get_granges_from_text(file.path(getwd(), "datasets", "ch2", "arabidopsis_chr4.txt"))

通过对属性进行筛选来提取区域；在此案例中——seqnames和metadata列：

genes_on_chr4 <- gr_from_gff[ gr_from_gff$type == "gene" & seqnames(gr_from_gff) %in% c("Chr4") ]

手动创建一个感兴趣的区域：

region_of_interest_gr <- GRanges(
    seqnames = c("Chr4"), 
    IRanges(c(10000), width= c(1000))
)

使用感兴趣的区域对子集进行筛选：

overlap_hits <- findOverlaps(region_of_interest_gr, gr_from_gff)
features_in_region <- gr_from_gff[subjectHits(overlap_hits) ]
features_in_region

它是如何工作的……

这里的第一步是创建一个GRanges对象，该对象描述了你感兴趣的基因组特征。我们创建的三个函数分别加载来自不同文件类型的信息，即.gff、.bed和一个制表符分隔的.txt文件，并返回所需的GRanges对象。在第 2 步中，我们利用 GFF 和文本函数创建了两个GRanges对象：gr_from_gff和gr_from_txt。这些对象随后用于子集操作。首先，在第 3 步中，我们对特征属性进行了子集操作。代码在 4 号染色体上找到了基因类型的特征。请注意，在Chr4中寻找基因和特征的语法区别。GRanges对象中的基本列—即seqnames、width和start—都有访问器函数来返回向量。

因此，我们在条件的第二部分使用了它。所有其他列—在GRanges术语中称为元数据—可以使用标准的$语法访问，因此我们在条件的第一部分使用了它。

在第 4 步中，我们在自定义的最小GRanges对象中创建了一个特定的区域。这个对象只包含一个区域，但可以通过在手动指定的向量中添加更多的seqnames、start和width来添加更多区域。最后，在第 5 步中，我们使用findOverlaps()函数获取在gr_from_gff对象中与手动创建的region_of_interest重叠的特征的索引，并使用这些索引来子集较大的gr_from_gff对象。

这将产生以下输出：

## GRanges object with 1 range and 10 metadata columns:
##       seqnames     ranges strand |   source     type     score     phase
##          <Rle>  <IRanges>  <Rle> | <factor> <factor> <numeric> <integer>
##   [1]     Chr4 2895-10455      - |   TAIR10     gene      <NA>      <NA>
##                ID        Name                Note          Parent
##       <character> <character>     <CharacterList> <CharacterList>
##   [1]   AT4G00020   AT4G00020 protein_coding_gene            <NA>
##             Index Derives_from
##       <character>  <character>
##   [1]        <NA>         <NA>
##   -------
##   seqinfo: 1 sequence from an unspecified genome; no seqlengths

请注意，我们需要使用subjectHits()访问器来提取 subject hits 列。

还有更多...

通过利用GRanges，也可以以相同的方式提取数据框或矩阵的子集，GRanges是其他对象的一部分。在以下示例中，我们创建了一个随机数据矩阵，并用它构建了一个SummarizedExperiment对象，该对象使用GRanges对象来描述其行：

set.seed(4321)
experiment_counts <- matrix( runif(4308 * 6, 1, 100), 4308)
sample_names <- c(rep("ctrl",3), rep("test",3) )
se <- SummarizedExperiment::SummarizedExperiment(rowRanges = gr_from_txt, assays = list(experiment_counts), colData = sample_names)

然后，我们可以像之前一样进行子集操作，并返回数据的子集以及范围的子集。assay()函数返回实际的数据矩阵：

overlap_hits <- findOverlaps(region_of_interest_gr, se)
data_in_region <- se[subjectHits(overlap_hits) ]
assay(data_in_region)

这将给出结果输出：

##          [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
## [1,] 69.45349 90.44524 88.33501 60.87932 86.24007 45.64919

通过 GWAS 寻找表型和基因型关联

使用高通量测序技术在许多样本中找到数千种基因变异的强大应用之一是全基因组关联研究（GWAS）——基因型和表型的关系。GWAS 是一项基因组分析，研究不同个体或遗传群体中的遗传变异，查看是否某个特定变异与某一性状相关。进行这种研究的方法有很多，但都依赖于收集特定样本中的变异数据，并在某种方式下与表型进行交叉参考。在本例中，我们将研究 2006 年 Yu et al 提出的复杂混合线性模型（Nature Genetics，38:203-208）。描述统一混合线性模型的工作原理超出了本例的范围，但它是一个适用于样本量大且具有广泛等位基因多样性的数据的模型，且可用于植物和动物数据。

准备工作

在本例中，我们将构建运行分析所需的数据结构，这些数据结构来自输入的 VCF 文件。我们将使用 GWAS() 函数，该函数位于 rrBLUP 包中。我们的样本数据文件包含三个 SNP——出于教学目的，这有助于我们的编程任务，但对于 GWAS 研究来说，数量显然过于微小。虽然代码可以运行，但结果在生物学上没有实际意义。

我们需要 rrBLUP，它不是 Bioconductor 的一部分，因此需要通过 install.packages() 安装，同时还需要 VariantAnnotation 和 datasets/ch2/small_sample.vcf 文件。

如何实现...

通过 GWAS 查找表型和基因型关联可以通过以下步骤完成：

加载库并获取 VCF 文件：

library(VariantAnnotation)
library(rrBLUP)
set.seed(1234)
vcf_file <- file.path(getwd(), "datasets", "ch2", "small_sample.vcf")
vcf <- readVcf(vcf_file, "hg19")

提取基因型、样本和标记位置信息：

gts <- geno(vcf)$GT

samples <- samples(header(vcf))
markers <- rownames(gts)
chrom <- as.character(seqnames(rowRanges(vcf)))
pos <- as.numeric(start(rowRanges(vcf)))

创建一个自定义函数，将 VCF 基因型转换为 GWAS 函数使用的约定：

convert <- function(v){
  v <- gsub("0/0", 1, v)
  v <- gsub("0/1", 0, v)
  v <- gsub("1/0", 0, v)
  v <- gsub("1/1",-1, v)
  return(v)
}

调用函数并将结果转换为数值矩阵：

gt_char<- apply(gts, convert, MARGIN = 2)

genotype_matrix <- matrix(as.numeric(gt_char), nrow(gt_char) )
colnames(genotype_matrix)<- samples

构建描述变异的数据框：

variant_info <- data.frame(marker = markers,
                           chrom = chrom,
                           pos = pos)

构建一个合并的变异/基因型数据框：

genotypes <-  cbind(variant_info, as.data.frame(genotype_matrix))
genotypes

构建一个 phenotype 数据框：

phenotypes <- data.frame(
  line = samples,
  score = rnorm(length(samples))
                         )

phenotypes

运行 GWAS：

GWAS(phenotypes, genotypes,plot=FALSE)

它是如何工作的...

本例中的大部分代码都是设置代码。在加载库并使用 set.seed() 固定随机数生成器以确保结果可重复之后，第一步是加载有用变异的 VCF 文件，第二步是提取一些有用的信息：我们获取一个基因型矩阵，使用 geno(vcf)$GT 调用，它返回一个矩阵，其中每一行是一个变异，每一列是一个样本，基因型记录在交点处。然后，我们使用一些访问器函数提取样本和标记名称，以及每个变异的参考序列（chrom）和位置（pos）。在步骤 3 中，我们定义一个转换函数（convert()），将 VCF 风格的杂合和纯合标注映射到 GWAS() 使用的格式。简而言之，在 VCF 中，"0/0" 表示 AA（纯合），在 GWAS() 中编码为 1，"0/1" 和 "1/0" 表示杂合 Aa 或 0，而 "1/1" 表示纯合 aa 或 -1。

在步骤 4中，我们将convert()应用于gts矩阵。烦人的是，返回值是一个字符矩阵，必须转换为数值型并重新包装成矩阵，这就是步骤 4中最后几行的目的。

在步骤 5中，我们构建一个数据框，描述我们之前创建的样本、标记和序列信息中的变异信息，而在步骤 6中，我们实际上将变异信息与基因型编码结合起来。

这将产生以下输出：

##           marker chrom     pos NA00001 NA00002 NA00003
## 1      rs6054257    20   14370       1       0      -1
## 2   20:17330_T/A    20   17330       1       0       1
## 3 20:1230237_T/G    20 1230237       1       1       0

请注意，列的顺序很重要。GWAS()函数期望我们按照这里指定的顺序提供这些信息。

在步骤 7中，我们构建了表型信息。第一列必须命名为line，并包含与基因型矩阵列相同顺序的样本名称。其余列可以是表型分数，并具有固定效应。

这将产生类似以下的输出（如果你在脚本顶部省略了set.seed()调用，由于随机化过程和示例数据中的小样本量，你的实际数字可能会有所不同）：

##      line     score
## 1 NA00001 -1.2070657
## 2 NA00002 0.2774292
## 3 NA00003 1.0844412

最后，在步骤 8中，我们运行GWAS()函数。

这将产生以下输出（同样，你的数字可能会有所不同）：

## [1] "GWAS for trait: score"
## [1] "Variance components estimated. Testing markers."

##           marker chrom     pos      score
## 1      rs6054257    20   14370 0.3010543
## 2   20:17330_T/A    20   17330 0.3010057
## 3 20:1230237_T/G    20 1230237 0.1655498

默认情况下，函数尝试生成一个图形。由于数据点太少，无法生成图形，因此我们在此通过plot = FALSE关闭该功能。

估计一个感兴趣位点的拷贝数

通常，我们会关心一个序列在样本中的出现频率——也就是估计在你的特定样本中，某个位点是否发生了重复，或其拷贝数是否增加。该位点可以是 Kbp 尺度上的基因，也可以是 Mbp 尺度上的较大 DNA 片段。我们在这个食谱中的方法是，在比对后的 HTS 读取覆盖度基础上估算背景覆盖度，然后检查我们感兴趣区域的覆盖度。我们感兴趣区域的覆盖度与背景水平的比率将为我们提供该区域的拷贝数估计。这个食谱是第一步。我们使用的背景模型非常简单——仅计算一个全局均值，但稍后我们将讨论一些替代方法。此外，这个食谱不涉及倍性——即细胞中整个基因组的拷贝数。通过类似的数据（尤其是 SNP 的主要/次要等位基因频率）可以估算倍性，但这是一个非常复杂的流程。请查看另请参阅部分，了解关于用于该长期分析的包的建议。

准备工作

对于这个食谱，我们需要csaw Bioconductor 包和样本hg17人类基因组.bam文件，文件位于datasets/ch2/hg17_snps.bam。

如何实现...

估计一个感兴趣位点的拷贝数可以通过以下步骤完成：

加载库并获取基因组窗口中的计数：

library(csaw)
whole_genome <- csaw::windowCounts( 
    file.path(getwd(), "datasets", "ch2", "hg17_snps.bam"), 
    bin = TRUE, 
    filter = 0, 
    width = 100, 
    param = csaw::readParam( minq = 20, dedup = TRUE, pe = "both" ) 
) 
colnames(whole_genome) <- c("h17")

从SummarizedExperiment中提取数据：

counts <- assay(whole_genome)[,1]

计算一个低计数阈值，并将计数较低的窗口设置为NA：

min_count <- quantile(counts, 0.1)[[1]]
counts[counts < min_count] <- NA

在中间的一组窗口中翻倍计数——这些将作为我们的高拷贝数区域：

n <- length(counts)
doubled_windows <- 10

left_pad <- floor( (n/2) - doubled_windows )
right_pad <- n - left_pad -doubled_windows
multiplier <- c(rep(1, left_pad ), rep(2,doubled_windows), rep(1, right_pad) )
counts <- counts * multiplier

计算每个窗口的平均覆盖度以及该窗口与平均覆盖度的比值，并通过绘图检查该比值向量：

 mean_cov <- mean(counts, na.rm=TRUE) 
 ratio <- matrix(log2(counts / mean_cov), ncol = 1)
 plot(ratio)

使用新数据和旧数据的行数据构建SummarizedExperiment：

se <- SummarizedExperiment(assays=list(ratio), rowRanges= rowRanges(whole_genome), colData = c("CoverageRatio"))

创建感兴趣区域并提取其中的覆盖数据：

region_of_interest <- GRanges(
  seqnames = c("NC_000017.10"),
  IRanges(c(40700), width = c(1500) )
)

overlap_hits <- findOverlaps(region_of_interest, se)
data_in_region <- se[subjectHits(overlap_hits)]
assay(data_in_region)

它是如何工作的...

在步骤 1中，本配方以熟悉的方式开始，使用csaw包在我们的人类第 17 号染色体小段上获取 100 bp 窗口的读取计数。读取过滤选项在param参数中设置。在步骤 2中，我们提取数据的第一列（也是唯一一列），使用assay()函数和子集提取来获得一个简单的计数向量。接下来，在步骤 3中，我们使用quantile()函数从counts向量中获取低于 10^(百分位)的min_count值。需要使用双括号子集提取，以便从quantile()函数返回的命名向量中获得单一数字。min_count值将作为切断点。counts向量中低于该值的所有值都将被设置为NA，以便从分析中移除——这充当了一个低覆盖阈值，并且可以根据需要在您自己修改的配方中调整使用的百分位数。

在步骤 4中，我们添加了一些覆盖度翻倍的区域——以便我们能检测到它们。我们选择了一些窗口，在其中翻倍计数，然后创建一个与计数长度相等的multiplier向量，其中1表示我们不希望更改计数，2表示我们希望将其翻倍。然后我们应用乘法。步骤 4可能会在您的分析中省略，因为它是一个合成数据生成步骤。

在步骤 5中，我们实际计算背景覆盖度。我们这里的函数是一个简单的全局平均，保存在mean_cov中——但您可以使用许多其他函数。有关此的讨论，请参见另见部分。我们还计算每个窗口计数与全局mean_cov的比值的log2()，并将其保存在一个名为ratio的单列矩阵对象中——因为我们需要将结果保存为矩阵形式，才能用于最终的SummarizedExperiment对象。我们快速使用plot()来检查ratio，并可以清晰地看到数据中间窗口计数翻倍的情况。

这将产生以下输出：

在步骤 6中，我们构建一个新的SummarizedExperiment对象se，用于保存窗口范围和新的比值数据。我们从window_counts中提取GRanges和colData对象，并添加我们新的ratio矩阵。现在我们可以开始对其进行子集提取，查看感兴趣区域的覆盖情况。

在第 7 步中，我们构造了一个手动的GRanges对象，表示我们感兴趣的任意区域，称为region_of_interest，并使用它通过findOverlaps()找到se对象中重叠的窗口。然后我们使用结果中的overlap_hits向量对子集se对象，并通过assay()函数查看该兴趣区域的计数。

这将产生以下输出：

##              [,1]
##  [1,]  0.01725283
##  [2,]  0.03128239
##  [3,] -0.05748994
##  [4,]  0.05893873
##  [5,]  0.94251006
##  [6,]  0.88186246
##  [7,]  0.87927929
##  [8,]  0.63780103
##  [9,]  1.00308550
## [10,]  0.75515798
## [11,]  0.80228189
## [12,]  1.05207419
## [13,]  0.82393626
## [14,]          NA
## [15,]          NA
## [16,] -0.16269298

在输出中，我们可以看到该区域相对于背景大约有一个 log2 比例为 1（二倍）的覆盖度，我们可以将其解读为拷贝数为 2。

另见

该配方中的背景水平计算非常简单——这对于学习配方是很好的，但在你自己真实数据中可能很快就不够用了。你可以采用许多选项来修改计算背景水平的方式，以适应你自己的数据。查看zoo包中的rollmeans()和rollmedians()函数——这些函数可以在任意步长的滚动窗口中给出均值和中位数，能为你提供一个移动窗口的背景平均值，可能会更加合适。

与拷贝数相关的分析之一是从 SNP 等位基因频率估计倍性。你可以查看vcfR包中的freq_peaks()函数，作为从BAM文件中的变异信息估计倍性的起点。

第三章：搜索基因和蛋白质的结构域与模体

基因、蛋白质和整个基因组的序列包含了它们功能的线索。重复的子序列或相互之间高度相似的序列，可能是进化保守性或功能相关性的线索。因此，针对模体和结构域的序列分析是生物信息学的核心技术。Bioconductor 包含了许多有用的包，用于分析基因、蛋白质和基因组。在本章中，你将学习如何使用 Bioconductor 分析序列中具有功能意义的特征，例如从广泛使用的数据库中提取的 de novo DNA 模体和已知结构域。你将学习一些基于核函数的机器学习包，用于发现蛋白质序列特征。你还将学习一些大规模比对技术，用于处理非常多或非常长的序列。你将使用 Bioconductor 和其他统计学习包。

本章将涉及以下配方：

使用 universalmotif 查找 DNA 模体
使用 PFAM 和 bio3d 查找蛋白质结构域
查找 InterPro 结构域
执行基因或蛋白质的多重比对
使用 DECIPHER 对基因组长度序列进行比对
用于蛋白质新特征检测的机器学习
使用 bio3d 进行蛋白质的 3D 结构比对

技术要求

你需要的样本数据可以从本书的 GitHub 仓库获取：github.com/danmaclean/R_Bioinformatics_Cookbook[.]（github.com/danmaclean/…

以下是你将需要的 R 包。大多数包可以通过 install.packages() 安装；其他包则稍微复杂一些：

ape
Bioconductor：
- Biostrings
- biomaRt
- DECIPHER
- EnsDb.Rnorvegicus.v79
- kebabs
- msa
- org.At.tair.db
- org.Eck12.db
- org.Hs.eg.db
- PFAM.db
- universalmotif
bio3d
dplyr
e1071
seqinr

Bioconductor 是一个庞大的工具集，具有自己的安装管理器。你可以使用以下代码进行安装：

if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install()

更多信息请访问 www.bioconductor.org/install/。

通常，在 R 中，用户会加载一个库并直接按名称使用函数。这在交互式会话中非常方便，但当加载了许多包时，可能会引起混淆。为了澄清我在某一时刻使用的是哪个包和函数，我会偶尔使用 packageName::functionName() 的约定。

有时，在一个过程的中途，我会中断代码，方便你查看一些中间输出或是需要理解的对象结构。每当发生这种情况时，你会看到一个代码块，其中每行前面都有 ## 双井号符号。请参考以下命令：

letters[1:5]

这将给我们如下的输出——请注意，输出行以##为前缀：

## a b c d e

本章中我们要使用的一些软件包依赖于必须单独安装的第三方软件。在 Windows、Linux 或 macOS 上安装和管理生物信息学软件的一种极好的方式是使用 conda 包管理器，并结合 bioconda 包频道。你可以通过一些简单的命令安装大量软件。要安装这两者，首先阅读bioconda.github.io/上的当前安装说明。

使用 universalmotif 查找 DNA motifs

在处理 DNA 序列时，一个非常常见的任务是查找 motif 实例——即在较长序列中定义的短序列。这些可能代表蛋白质- DNA 结合位点，例如基因启动子中的转录因子结合位点或增强区域。在此分析中有两种起点：要么你有一个 motif 数据库，想用它来扫描目标 DNA 序列并提取出现的 motif，要么你只有感兴趣的序列，想要查找其中是否有任何重复的 motif。在本食谱中，我们将讨论如何执行这两项操作。我们将在这两种情况下都使用universalmotif包。

准备工作

本食谱中，我们需要datasets/ch3/simple_motif.txt和datasets/ch3/promoters.fa文件，一个描述简单 motif 的简单矩阵，格式为位置特异性权重矩阵（PSWM）（简要描述请参见附录），以及一组来自转录起始位点上游的序列。

本食谱还需要在你的系统中有一个有效的MEME副本。MEME是一个用于在序列集中过度表现的统计序列 motif 的程序。当用于启动子或基因上游区域时，这些 motif 可能代表转录因子结合位点。MEME的网页地址为alternate.meme-suite.org/，如果你已经安装了 conda，可以通过conda install -c bioconda meme来安装。MEME包无法在 Windows 系统上使用。如果你希望在 Windows 上运行它，可以考虑在 Cygwin（Linux 模拟层）下运行它：www.cygwin.com/。你可能还需要在 Cygwin 下安装新的 R 版本。

如何操作...

使用universalmotif查找 DNA motifs 可以按照以下步骤进行：

首先，加载所需的库和感兴趣的 motif：

library(universalmotif)
library(Biostrings)

motif <- read_matrix(file.path(getwd(), "datasets", "ch3","simple_motif.txt"))

然后，加载要用该 motif 扫描的序列：

sequences <- readDNAStringSet(file.path(getwd(), "datasets", "ch3", "promoters.fa"))

执行序列扫描：

motif_hits <- scan_sequences(motif, sequences = sequences)
motif_hits

请注意，motif_hits包含关于 motif 在每个目标序列中位置的信息。

计算该结构是否在序列中富集：

motif_info <- enrich_motifs(motif, sequences, shuffle.k = 3, verbose = 0, progress = FALSE, RC = TRUE)
motif_info

请注意，motif 信息包含关于在一组序列中统计富集的信息。

运行MEME以查找新型 motif：

meme_path = "/Users/macleand/miniconda2/bin/meme"
meme_run <- run_meme(sequences, bin = meme_path, output = "meme_out", overwrite.dir = TRUE)
motifs <- read_meme("meme_out/meme.txt")
view_motifs(motifs)

工作原理...

这段代码真是太棒了！仅用几行代码，我们就完成了整个分析。我们首先加载了一个基序的矩阵描述和一些我们希望找到启动子的序列——这发生在步骤 1 和 2 中，结果我们得到了一个universalmotif对象和一个DNAStringSet对象。接下来的真正工作发生在步骤 3 和 4 中。scan_sequences()函数会搜索每个序列并报告找到基序的位置——查看motif_hits对象，看看它们在哪里。

这将产生以下输出：

##               motif  sequence start stop  score max.score score.pct 
## 1   YTTTYTTTTTYTTTY AT4G28150    73   87  7.531  22.45824  33.53335 
## 2   YTTTYTTTTTYTTTY AT4G28150    75   89 10.949  22.45824  48.75270

在判断一个基序是否显著时，universalmotifs包中的enrich_motifs()函数会在第 4 步为我们完成这项工作，结果会产生以下输出：

##             motif total.seq.hits num.seqs.hit num.seqs.total
## 1 YTTTYTTTTTYTTTY            916           50             50
##   total.bkg.hits num.bkg.hit num.bkg.total   Pval.hits   Qval.hits
## 1            265          48            50 4.75389e-85 4.75389e-85
##     Eval.hits
## 1 9.50778e-85

它会搜索序列，寻找可能出现基序的地方，并对其进行计数，执行 Fisher 精确检验，比较我们的序列集中的基序频率与自动生成的背景集中的基序频率。最终的motif_info输出包含一个p值报告。为了找到新的基序，我们在第 5 步运行外部软件MEME。run_meme()函数需要知道MEME包在你系统中的位置，因此我们需要在meme_path变量中定义该路径。

请注意，系统中meme_path的值与此处提到的值会有所不同——这是我系统中的一个示例。

我们将这些信息传递给函数，并提供包含序列的DNAStringSet对象。该函数还需要一个输出目录，用于写入MEME的结果，因为它不会返回任何有用的信息给 R。run_meme()函数会在后台执行MEME，一旦运行完成，我们可以通过read_meme()函数和文件名从meme.txt文件中加载结果。它会返回一个universalmotif对象。最后，在这里，我们通过view_motifs()函数快速查看motifs对象：

这会给我们一个非常直观的基序可视化。

还有更多内容...

从现有的数据库（如 JASPAR 和 TRANSFAC）中加载基序非常简单，使用universalmotif时，只需对read_matrix()函数做直接替换。查看以下函数，它们可以从不同格式中加载基序：read_cisbp()、read_homer()、read_jaspar()、read_matrix()、read_meme()、read_motifs() 和 read_uniprobe()。

使用 PFAM 和 bio3d 查找蛋白质结构域

发现蛋白质序列的功能是一个关键任务。我们可以通过多种方式来完成这项任务，包括使用 BLAST 等工具对已知蛋白质数据库进行全序列相似性搜索。如果我们希望获得更具信息量和更精细的资料，可以选择在序列中查找个体功能域。Pfam等数据库和hmmer等工具使得这一过程成为可能。Pfam将蛋白质结构域编码为配置文件隐马尔可夫模型，而hmmer则使用这些模型扫描序列并报告任何可能存在的结构域。通常，基因组注释项目会为我们完成这些搜索，这意味着要在我们的序列中找到Pfam结构域，关键就在于搜索数据库。Bioconductor 在这些特定数据库的数据包装方面做得非常出色，通常是以**.db**为后缀的包。在本配方中，我们将展示如何判断一个包是否包含 Pfam 结构域信息，如何为特定基因提取这些信息，以及如果没有现成的信息时，如何自行运行 Pfam 搜索。

正在准备中

对于这个示例，我们需要一些 Bioconductor 的Annotationdbi数据库包—具体来说，org.Hs.eg.db、org.EcK12.eg.db和org.At.tair.db。

你还需要bio3d包，当前版本—在写这篇文档时—只有在使用开发版时，bio3d才会连接到 Pfam 服务器。你可以通过devtools包从 BitBucket 安装该版本：

install.packages("devtools")
library(devtools)
install_bitbucket("Grantlab/bio3d", subdir = "ver_devel/bio3d/")

如何操作...

使用PFAM.db和bio3d查找蛋白质结构域，可以通过以下步骤完成：

加载数据库包并检查数据库中的键类型：

library(org.Hs.eg.db)
keytypes(org.Hs.eg.db)

注意这个输出中的ENSEMBL键—我们可以用它来查询数据库。

使用keys()函数获取键的向量：

k <- head(keys(org.Hs.eg.db, keytype = "ENSEMBL"), n = 3 )

查询数据库：

result <- select(org.Hs.eg.db, keys = k, keytype="ENSEMBL", columns = c("PFAM"))
result

加载PFAM数据库以提取描述信息：

library(PFAM.db)
descriptions <- PFAMDE

从PFAM数据库中获取所有键：

all_ids <- mappedkeys(descriptions)

获取PFAM ID 的所有描述：

id_description_mapping <- as.data.frame(descriptions[all_ids])

将描述信息与PFAM关联：

dplyr::left_join(result, id_description_mapping, by = c("PFAM" = "ac") )

工作原理...

这种方法的关键在于找出我们使用的数据库是否确实包含 PFAM 域信息。这就是我们在步骤 1 中做的—我们使用keytypes()函数列出可用的搜索键。PFAM 出现在结果中。一旦我们确认可以使用此数据库来获取我们需要的信息，就可以遵循一个相当标准的流程：

获取查询用的键列表—比如基因名。这里我们直接从数据库中提取它们，但它们也可以来自其他地方。这将产生以下输出：

##  [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT" 
##  [5] "ENSEMBLTRANS" "ENTREZID"     "ENZYME"       "EVIDENCE"    
##  [9] "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"       
## [13] "IPI"          "MAP"          "OMIM"         "ONTOLOGY"    
## [17] "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"        
## [21] "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"     
## [25] "UNIGENE"      "UNIPROT"

使用select()函数查询数据库，该函数根据提供的键拉取数据。columns参数告诉它要拉取哪些数据。这里的表达式将获取我们感兴趣基因的 PFAM ID。
制作所有 PFAM ID 和描述的列表。我们加载 PFAM.db 包，并使用它提供的 PFAMDE 对象来获取 ID 和描述之间的映射。这将产生以下输出。请注意，因为我们从外部数据库中提取数据，数据库的更改可能会在这里反映出来：

##           ENSEMBL    PFAM
## 1 ENSG00000121410 PF13895
## 2 ENSG00000175899 PF01835
## 3 ENSG00000175899 PF07678
## 4 ENSG00000175899 PF10569
## 5 ENSG00000175899 PF07703
## 6 ENSG00000175899 PF07677
## 7 ENSG00000175899 PF00207
## 8 ENSG00000256069    <NA>

然后，我们可以使用 mappedkeys() 函数获取实际的描述信息对象。
接下来，提取和转换 all_ids 对象的描述为数据框。
最后，我们将 PFAM 域的描述与之前获得的 PFAM ID 进行连接，使用具有共同数据的列——PFAM 和 ac。这将产生以下输出：

##           ENSEMBL    PFAM                                       de
## 1 ENSG00000121410 PF13895                    Immunoglobulin domain
## 2 ENSG00000175899 PF01835                               MG2 domain
## 3 ENSG00000175899 PF07678               A-macroglobulin TED domain
## 4 ENSG00000175899 PF10569                                     <NA>
## 5 ENSG00000175899 PF07703 Alpha-2-macroglobulin bait region domain
## 6 ENSG00000175899 PF07677  A-macroglobulin receptor binding domain
## 7 ENSG00000175899 PF00207             Alpha-2-macroglobulin family
## 8 ENSG00000256069    <NA>                                     <NA>

还有更多...

我提到食谱的关键部分——特别是第 6 步的连接——是确保数据库包含正确的键，特别是 PFAM，以便继续进行。根据生物体和数据库的不同，PFAM 注释可能不存在。以下是如何检查你感兴趣的数据库中是否存在它，使用了两个示例数据库 org.At.tair.db 和 org.Eck12.eg.db，一个拟南芥数据库：

library(org.At.tair.db)
columns(org.At.tair.db)

还需要一个 E.coli 数据库：

library(org.EcK12.eg.db)
columns(org.EcK12.eg.db)

只需使用 columns() 函数报告数据库中的数据列。如果 PFAM 出现，你可以按照以下步骤进行。如果未出现，则可以作为替代方案，在服务器上运行 PFAM 并自行进行注释。以下代码使用 bio3d 函数 hmmer() 在 EBI 服务器上对输入的蛋白质序列进行 PFAM 搜索。返回的对象在 hit.tbl 槽中包含 PFAM 输出的数据帧：

sequence <- read.fasta(file.path(getwd(), "datasets", "ch3", "ecoli_hsp.fa") )
# run pfamseq on protein
result <- hmmer(sequence, type="hmmscan", db="pfam")
result$hit.tbl

这将产生以下输出：

##   name        acc bias dcl desc  evalue flags hindex ndom nincluded
## 1 GrpE PF01025.19  3.3 272 GrpE 1.4e-46     3   8846    1         1
##   nregions nreported    pvalue score taxid     pdb.id bitscore mlog.evalue
## 1        1         1 -115.4076 158.2     0 PF01025.19    158.2    105.5824

查找 InterPro 域

InterPro 是由多个蛋白质数据库提供的预测模型或签名的数据库。InterPro 聚合来自多个来源的信息，以减少注释的冗余性并帮助解释性。在这个示例中，我们将扩展我们用于 PFAM 域的方法，查看感兴趣序列上的 InterPro 域的注释。我们将从 Ensembl 核心数据库开始。

准备工作

我们需要 ensembldb、Ensdb.Rnorvegicus.v79 和 biomaRt 这几个 Bioconductor 包。

如何操作...

可以通过以下步骤找到 InterPro 蛋白质域：

加载库并仔细检查我们的数据库包是否携带所需的蛋白质数据：

library(ensembldb)
library(EnsDb.Rnorvegicus.v79)
hasProteinData(EnsDb.Rnorvegicus.v79)

构建要查询的基因列表——请注意，这里我需要的 keytype 是 GENEID：

e <- EnsDb.Rnorvegicus.v79
k <- head(keys(e, keytype = "GENEID"), n = 3 )

使用 select() 函数提取相关数据：

select(e, keys = GeneIdFilter(k),
       columns = c("TXBIOTYPE", "UNIPROTID", "PROTEINID","INTERPROACCESSION"))

工作原理...

代码是通过相关包在特定的 Rattus norvegicus Ensembl 数据库中进行数据库查找。该过程与 PFAM 域搜索类似：

我们使用 EnsemblDB 包特定的 hasProteinData() 函数来检查数据库是否包含我们需要的信息。如果输出是 TRUE，则很好：

## [1] TRUE

我们再次构建感兴趣的基因列表——这里我从数据库中提取列表，但这些 ID 可以来自任何地方。
最后，我们用感兴趣的基因作为关键字搜索数据库。请注意，我们需要GeneIdFilter()函数包装器和columns参数来选择要返回的数据。这将导致一个包含以下信息的数据框架：

##         TXBIOTYPE UNIPROTID          PROTEINID INTERPROACCESSION  GENEID
## 1  protein_coding    Q32KJ7 ENSRNOP00000052495         IPR017850  ENSRNOG00000000001
## 2  protein_coding    Q32KJ7 ENSRNOP00000052495         IPR000917  ENSRNOG00000000001
## 3  protein_coding    C9E895 ENSRNOP00000000008         IPR015424  ENSRNOG00000000007

这里还有更多内容……

我们在这个配方中使用的方法非常适合 Ensembl 核心数据库，但还有其他非 Ensembl 核心数据库，我们可能想要搜索；对此，有 biomaRt。 biomaRt 允许我们定义到我们可能知道的其他数据库的连接。许多这些数据库公开了我们可以用来查询它们的 API。为此，加载biomaRt库并使用useMart()函数来定义与适当主机和数据集的连接。然后，使用连接和列以及基因 ID 使用getBM()函数进行查询。如果您的查询是interpro，则会返回 InterPro 的搜索结果。以下示例在plants.ensembl.org上为两个拟南芥基因进行搜索：

library(biomaRt)
biomart_athal <- useMart(biomart = "plants_mart", host = "plants.ensembl.org", dataset = "athaliana_eg_gene")
getBM( c("tair_locus", "interpro"), filters=c("tair_locus"), values = c("AT5G40950", "AT2G40510"), mart = biomart_athal)

这将返回以下输出：

##   tair_locus  interpro
## 1  AT2G40510 IPR000892
## 2  AT2G40510 IPR038551
## 3  AT5G40950 IPR001684
## 4  AT5G40950 IPR018261

另请参阅……

如果您在确定 marts 和 columns 的名称时遇到问题，请尝试来自bioMart的listMarts()和listDatasets()函数，它们将提供当前可用 marts 和它们包含的数据列表。

执行基因或蛋白质的多重比对

在构建系统发育树之前或作为确定保守和分散区域的结束任务中，对序列进行对齐是生物信息学分析的重要组成部分，并且在 R 和 Bioconductor 中使用ape和DECIPHER包广泛涵盖了从序列到对齐的极为简单的程序。在这个配方中，我们将看看如何实现这一过程。

请注意，不同的序列长度有不同的技术。在这个第一个配方中，我们将看看如何对齐一些 Kbp 长度的序列，例如代表基因和蛋白质的序列。

准备就绪

此配方需要msa包。这是一个相当庞大的包，包括外部软件：Clustal、Clustal Omega 和 Muscle。还需要ape和seqinR包。作为测试数据集，我们将使用存储在书籍数据和代码存储库中的一些血红蛋白蛋白质序列，路径为datasets/ch3/hglobin.fa。您的系统中还需要 PDFLatex。您可以在此处找到安装信息：www.latex-project.org/get/。

如何做到这一点……

可以通过以下步骤执行基因或蛋白质的多重比对：

加载库和序列：

library(msa)
seqs <- readAAStringSet(file.path(getwd(), "datasets", "ch3", "hglobin.fa"))
seqs

执行多序列比对：

alignment <- msa(seqs, method = "ClustalOmega")
alignment

这将返回一个如下所示的对齐对象：

## ClustalOmega 1.2.0 
## 
## Call:
##    msa(seqs, method = "ClustalOmega")
## 
## MsaAAMultipleAlignment with 3 rows and 142 columns
##     aln                                               names
## [1] MVLSPADKTNVKAAWGKVGAHAG...PAVHASLDKFLASVSTVLTSKYR HBA_HUMAN
## [2] MVLSGEDKSNIKAAWGKIGGHGA...PAVHASLDKFLASVSTVLTSKYR HBA_MOUSE
## [3] MSLTRTERTIILSLWSKISTQAD...ADAHAAWDKFLSIVSGVLTEKYR HBAZ_CAPHI
## Con MVLS??DKTNIKAAWGKIG?HA?...PAVHASLDKFLASVSTVLTSKYR Consensus

使用以下代码查看结果：

msaPrettyPrint(alignment, output="pdf", showNames="left",
showLogo="none", askForOverwrite=FALSE, verbose=FALSE, file="whole_align.pdf")

使用以下代码查看放大区域：

msaPrettyPrint(alignment, c(10,30), output="pdf", showNames="left",
file = "zoomed_align.pdf", showLogo="top", askForOverwrite=FALSE, verbose=FALSE)

如何运行它……

这里的配方简单而明了——使用msa执行 MSA 非常简单。在步骤 1 中，我们使用常见的readAAStringSet()函数加载包和序列，以给出seqs——一个AAStringSet对象，我们可以检查并获取以下输出：

## A AAStringSet instance of length 3
##     width seq                                          names               
## [1]   142 MVLSPADKTNVKAAWGKVGAH...HASLDKFLASVSTVLTSKYR HBA_HUMAN
## [2]   142 MVLSGEDKSNIKAAWGKIGGH...HASLDKFLASVSTVLTSKYR HBA_MOUSE
## [3]   142 MSLTRTERTIILSLWSKISTQ...HAAWDKFLSIVSGVLTEKYR HBAZ_CAPHI

接下来，在步骤 2中，将msa()函数传入seqs对象和对齐方法的名称。这里我们使用ClustalOmega（你可以选择ClustalOmega、ClustalW或Muscle）。method参数指定了用于实际对齐的外部程序的名称。对齐器运行后，你将获得一个MsaMultipleAlignment对象——这是一个容器，包含了对齐后的序列，格式如下：

## ClustalOmega 1.2.0 
## 
## Call:
##    msa(seqs, method = "ClustalOmega")
## 
## MsaAAMultipleAlignment with 3 rows and 142 columns
##     aln                                               names
## [1] MVLSPADKTNVKAAWGKVGAHAG...PAVHASLDKFLASVSTVLTSKYR HBA_HUMAN
## [2] MVLSGEDKSNIKAAWGKIGGHGA...PAVHASLDKFLASVSTVLTSKYR HBA_MOUSE
## [3] MSLTRTERTIILSLWSKISTQAD...ADAHAAWDKFLSIVSGVLTEKYR HBAZ_CAPHI
## Con MVLS??DKTNIKAAWGKIG?HA?...PAVHASLDKFLASVSTVLTSKYR Consensus

在步骤 3 中，我们使用msaPrettyPrint()函数将对齐的可视化结果写入 PDF 文件。该函数有多个参数，用于描述对齐图的布局。可视化结果必须写入文件，不能像普通的绘图对象那样发送到交互式会话的绘图窗口中。图片存放的文件由file参数指定。生成的图片如下：

在步骤 4 中，我们使用第二个位置参数，通过c(10,30)向量将视图限制在 10 到 30 之间。我们得到以下图片：

不幸的是，由于图片制作过程在后台使用了 Latex，因此我们无法将图片强制转换为比 PDF 更有用的格式，或像其他绘图对象一样进行渲染。

还有更多...

在这个阶段，树形可视化的序列相似性常常非常有用。我们可以使用ape和seqinr包生成这种可视化。我们可以将对齐对象转换为描述序列距离的seqinr distance对象，并从中使用ape创建一个简单的邻接树，然后进行绘制：

library(ape)
library(seqinr)
alignment_seqinr <- msaConvert(alignment, type="seqinr::alignment")
distances <- seqinr::dist.alignment(alignment_seqinr, "identity")
tree <- ape::nj(distances)
plot(tree, main="Phylogenetic Tree of HBA Sequences")

这将产生以下输出：

使用 DECIPHER 对齐基因组长度序列

对齐比基因和蛋白质更长的序列，例如来自组装项目的 contigs、染色体或整个基因组，是一个棘手的任务，并且需要不同于短序列的技术。序列越长，越难对齐。长序列的对齐尤其在计算时间上非常耗费，因为对于短序列有效的算法，随着序列长度的增加，其计算时间呈指数级增长。执行长序列的对齐通常从寻找短的锚对齐开始，然后从那里进行扩展。我们通常会得到同源块——在不同基因组对齐之间匹配较好的区域。

在这个教程中，我们将介绍DECIPHER包用于基因组长度的对齐。我们将使用一些叶绿体基因组——约 150 Kbp 长的小型细胞器基因组，这些基因组在进化上相对保守，作为我们的目标序列。

准备就绪

确保你已经安装了DECIPHER包。我们将使用datasets/ch3/plastid_genomes.fa文件作为示例。

如何进行...

使用DECIPHER对齐基因组长度序列可以按照以下步骤进行：

加载库和基因组序列：

library(DECIPHER)
long_seqs <- readDNAStringSet(file.path(getwd(), "datasets", "ch3", "plastid_genomes.fa"))
long_seqs

在本地数据库中准备序列：

Seqs2DB(long_seqs, "XStringSet", "long_db", names(long_seqs))

查找 synteny 区块：

synteny <- FindSynteny("long_db")
pairs(synteny)

这将创建一个同源区块的点图。

绘制 syntenic 区块：

plot(synteny)

现在，进行实际的比对：

alignment <- AlignSynteny(synteny, "long_db")

然后逐一保存配对比对结果：

blocks <- unlist(alignment[[1]])
writeXStringSet(blocks, "genome_blocks_out.fa")

它是如何工作的...

DECIPHER 包非常强大，因此在进行实际分析之前，需要进行一些设置。在步骤 1 中，我们加载库并将序列加载到 long_seqs，一个 DNAStringSet 对象中；但在步骤 2 中，我们为后续步骤构建了一个额外的磁盘 SQLite 数据库。这是通过 Seqs2DB() 函数完成的，该函数接受 long_seqs（输入类型为 XStringSet，DNAStringSet 的父类）、数据库的名称（long_db）以及一个序列名称的向量，我们通过 names() 函数提取该向量。一旦数据库构建完成，我们就可以在以下工作流中使用它：

使用 FindSynteny() 函数在数据库中查找同源区块。此操作将产生以下输出：

## A DNAStringSet instance of length 5
##      width seq                                         names               
## [1] 130584 GGCATAAGCTATCTTCCCAA...GATTCAAACATAAAAGTCCT NC_018523.1 Sacch...
## [2] 161592 ATGGGCGAACGACGGGAATT...AGAAAAAAAAATAGGAGTAA NC_022431.1 Ascle...
## [3] 117672 ATGAGTACAACTCGAAAGTC...GATTTCATCCACAAACGAAC NC_022259.1 Nanno...
## [4] 154731 TTATCCATTTGTAGATGGAA...TATACACTAAGACAAAAGTC NC_022417.1 Cocos...
## [5] 156618 GGGCGAACGACGGGAATTGA...TTTTGTAGCGAATCCGTTAT NC_022459.1 Camel...

使用同源区块作为种子，并通过 AlignSynteny() 函数执行实际的比对。

这些操作在步骤 3 和 5 中完成。FindSynteny() 需要数据库的名称；AlignSynteny() 需要 synteny 对象和数据库名称。

最后，我们可以输出结果。带有 synteny 对象的 pairs() 函数将创建一个同源区块的点图：

使用带有 synteny 对象的 plot() 函数，可以创建一个有用的热图，显示同源区块相对于第一个基因组的位置。跨基因组的相同颜色区域表示同源序列区域：

最后一步，步骤 6，是稍微复杂的保存过程。alignment 对象是一个 R 列表，每个成员表示一个比对——它本身也是一个列表。通过提取并对每个返回的元素使用 unlist()，你将得到一个可以用 writeXStringSet() 保存为典型 FASTA 比对的对象（blocks）。记住，你需要分别为 blocks 对象的每个成员执行此操作。

用于蛋白质中新特征检测的机器学习

有时，我们会有一个来自某些分析或实验的蛋白质序列列表，这些序列在某种程度上是生物学相关的——例如，它们可能都与同一个靶标结合——我们将希望确定这些蛋白质中负责该功能的部分。正如我们在前面配方中所做的，域和基序的寻找可以有所帮助，但前提是我们之前见过这些域，或者该序列特别保守，或者在统计学上被过度表示。另一种方法是尝试使用机器学习，在这种方法中，我们构建一个可以准确分类我们感兴趣的蛋白质的模型，然后利用该模型的特性来告诉我们哪些部分的蛋白质导致了这个分类。在这个配方中，我们将采用这种方法；具体来说，我们将训练一个支持向量机（SVM）。

准备工作

对于这个配方，我们需要kebabs和Biostrings、e1071以及readr库，并且两个输入数据文件。机器学习在有大量训练样本时效果最佳，但它们运行需要时间，因此我们有一个相对较小的输入，即 170 个大肠杆菌蛋白质，根据STRING数据库（string-db.org/）的数据显示，这些蛋白质与pfo蛋白有实验性互作证据。这些是正向训练样本。我们还需要负向训练样本——这又是 170 个与pfo蛋白没有互作证据的大肠杆菌蛋白，它们是随机选择的。所有蛋白质序列都存储在datasets/ch3/ecoli_proteins.fa文件中。与该文件一起，还有一个文本文件记录了每个蛋白质的类别。datasets/ch3/ecoli_protein_classes.txt是一个单列文本文件，描述了每个蛋白质的类别——“1”表示正向pfo互作，“-1”表示负向pfo互作。类别文件中的行索引与序列文件中的蛋白质索引相匹配。

如何实现...

蛋白质的新特征检测的机器学习可以通过以下步骤完成：

加载库和输入文件：

library(kebabs)
library(Biostrings) 
seqs <- readAAStringSet(file.path(getwd(), "datasets", "ch3", "ecoli_proteins.fa")) 
classes <- readr::read_csv(file.path(getwd(), "datasets", "ch3", "ecoli_protein_classes.txt"), col_names = TRUE)
classes <- classes$class

将数据划分为训练集和测试集：

num_seqs <- length(seqs) 
training_proportion <- 0.75 
training_set_indices <- sample(1:num_seqs, training_proportion * num_seqs) 
test_set_indices <- c(1:num_seqs)[-training_set_indices]

使用训练集构建模型：

kernel <- gappyPairKernel(k=1, m=3) 
model <- kbsvm(x=seqs[training_set_indices], y=classes[training_set_indices], kernel=kernel, pkg="e1071", svm="C-svc", cost=15)

使用模型预测测试集的类别：

predictions <- predict(model, seqs[test_set_indices]) 
evaluatePrediction(predictions, classes[test_set_indices], allLabels=c(1,-1) )

这将给出以下输出：

##     1 -1
## 1  36 23
## -1 10 19
## 
## Accuracy:              62.500% (55 of 88)
## Balanced accuracy:     61.749% (36 of 46 and 19 of 42)
## Matthews CC:            0.250
## 
## Sensitivity:           78.261% (36 of 46)
## Specificity:           45.238% (19 of 42)
## Precision:             61.017% (36 of 59)

检查序列的预测图谱：

seq_to_test <- seqs[[1]][1:10]
seq_to_test

这将给出以下输出：

## 10-letter "AAString" instance ## seq: MSFDIAKYPT

然后，使用以下代码绘制prediction_profile：

prediction_profile <-getPredictionProfile(seq_to_test, kernel, featureWeights(model), modelOffset(model) ) plot(prediction_profile)

它是如何工作的...

这里的第一步很简单：我们加载我们感兴趣的序列以及它们所属的类别。因为我们将ecoli_protein_classes.txt文件加载到一个数据框中，当我们需要一个简单的向量时，我们使用$子集操作符从数据框中提取classes列。这样做会返回我们需要的那个单列向量对象。之后，工作流程很简单：

决定数据中有多少用于训练，多少用于测试：在第 1 步中，我们选择将数据的 75%作为训练集，在创建training_proportion变量时进行设置。它与num_seqs一起用于sample()函数，以随机选择序列的索引并将其放入训练集。training_set_indices变量包含我们稍后将用来子集化数据的整数。最初，我们通过使用方括号[]和否定运算符-，创建一个互补的索引列表test_set_indices。基本上，这种构造方法是创建一个包含所有不在training_set_indices中的索引的向量的惯用方式。
构建和训练支持向量机模型：在第 2 步中，我们构建了分类模型。首先，我们选择一个内核，将输入数据映射到一个支持向量机可以学习的矩阵空间。在这里，它来自gappyPairKernel()函数——请注意，内核类型有很多种；这个内核非常适合序列数据。我们将kernel传递给kbsvm()函数，并将training_set_indices子集的序列作为x参数传入，training_set_indices子集的类别作为y参数传入。此函数中的其他参数决定了模型的具体类型、包和训练参数。这些选项非常多，而且它们对最终模型的效果有很大的影响。了解并进行科学实验，找到最适合你数据的设置是非常值得的。最终模型保存在model变量中。
在未见数据上测试模型：现在我们有了一个模型，我们可以用它来预测未见蛋白质的类别。这个阶段将告诉我们模型的效果如何。在第 3 步中，我们使用predict()函数与模型以及我们没有用来训练的序列（即test_set_indices中的序列）进行预测，并返回一个预测对象。通过将预测结果传递给evaluatePrediction()函数，并结合类别向量中的实际类别和所有可能类别标签的向量（使用allLabels参数），我们可以得到模型的准确率和其他指标的摘要。这里模型的准确率是 62%，这还算可以；比随机猜测要好。但我们数据集相对较小，且模型没有经过优化；如果做更多的工作，模型的表现可能会更好。请注意，如果你运行这段代码，可能会得到不同的结果。由于训练集序列的选择是随机的，模型的表现可能会略微变好或变差，取决于输入数据的具体内容。
估计序列的预测特征：为了实际找到在分类中重要的区域，可能也涉及蛋白质的功能，我们使用getPredictionProfile()函数来分析一个序列。在第四步中，我们用列表和双括号索引提取第一个序列的小片段（10 个氨基酸），并用单括号索引获取范围；例如，seqs[[1]][1:10]。我们这样做是为了在最后一步中的可视化更清晰。你也可以直接使用完整的序列。getPredictionProfile()函数需要kernel和model对象才能正常工作。它将输出以下结果：

##     1 -1
## 1  36 23
## -1 10 19
## 
## Accuracy:              62.500% (55 of 88)
## Balanced accuracy:     61.749% (36 of 46 and 19 of 42)
## Matthews CC:            0.250
## 
## Sensitivity:           78.261% (36 of 46)
## Specificity:           45.238% (19 of 42)
## Precision:             61.017% (36 of 59)

最后，我们可以plot()预测特征：该特征展示了每个氨基酸对整体决策的贡献，并增加了学习结果的可解释性。在这里，第四个残基D对该蛋白质的决策贡献非常大。通过在多个序列中检查这种贡献模式，可以阐明影响决策的模式。值得注意的是，你可能会看到与以下示例略有不同的图像——这是由于算法中的随机过程——这是你应纳入分析的一部分：确保任何显著差异不是由于运行代码时的随机选择造成的。在这个例子中，最强的贡献依然应来自于"D"：

使用 bio3d 进行 3D 结构蛋白质比对

两个分子模型之间的三维结构比对可以揭示两种蛋白质共有或独特的结构特性。这些结构特性可能暗示着进化或功能上的共性。在本教程中，我们将展示如何在三维空间中对两个蛋白质序列进行比对，并在 3D 渲染软件中查看它们。

准备工作

对于本节内容，我们至少需要两个外部软件——PyMOL 和 MUSCLE——一个是 3D 结构渲染程序，另一个是序列比对工具。

可以通过 conda 安装 MUSCLE，命令如下：

conda install -c bioconda muscle

MUSCLE 的一个版本已经与msa包一起安装，并且 bio3d 可以引用这个安装版本。我们将在本教程中使用msa版本。

PyMOL 对于可视化至关重要，可以通过 conda 安装，命令如下：

conda install -c schrodinger pymol

要查找该软件的安装路径，可以使用which pymol命令。

除此之外，你还需要两个包含人类和果蝇硫氧还蛋白结构的文件来进行操作：datasets/ch3/1xwc.pdb和datasets/ch3/3trx.pdb。

如何操作...

使用 bio3d 进行 3D 结构蛋白质比对的步骤如下：

加载库和 PDB 结构文件：

library(bio3d)
a <- read.pdb(file.path(getwd(), "datasets", "ch3" ,"1xwc.pdb"))
b <- read.pdb(file.path(getwd(), "datasets", "ch3", "3trx.pdb"))

然后，进行结构比对：

pdbs <- pdbaln(list("1xwc"=a,"3trx"=b), fit=TRUE, exefile="msa")

在 PyMOL 中启动并渲染比对：

pymol_path = "/Users/macleand/miniconda2/bin/pymol"
pymol(pdbs, exefile = pymol_path, type = "launch", as="cartoon")

它是如何工作的...

和往常一样，第一步是加载库，然后加载输入数据。这里在第 1 步，我们使用 read.pdb() 函数加载两个独立的 PDB 文件。在第 2 步，我们使用 pdbaln() 函数进行对齐。所有我们想要对齐的 PDB 对象首先被放入一个适当命名的 list 对象中。fit 参数设置为 TRUE，以便基于所有对齐的序列执行结构叠加（执行的是基于序列的叠加）。

exefile 参数告诉 pdbaln() 在哪里执行基于序列的对齐部分，如这里所示；"msa" 的值使用 msa 包中的对齐工具，但你也可以使用一个指向替代可执行文件的路径，或者将 exefile 替换为你有效的电子邮件地址，使用 web.args="your.email@something.org" 这种形式在 EBI 上通过网络进行对齐。

一旦我们有了一个 pdbs 中的对齐对象，就可以在 PyMOL 中可视化它。我们在 pymol_path 变量中设置 PyMOL 的路径，并将其与类型 "launch" 一起传递给 pymol() 函数，这将创建一个交互式的 PyMOL 会话。或者，省略 type 参数将生成一个 PyMOL 脚本，稍后可以使用。PyMol 应该会显示接下来的图片。截图显示了两种对齐蛋白质的渲染效果：人类版本为红色，果蝇版本为黄色。

还有更多...

pdbaln() 函数对于长度相似的结构效果很好。对于那些 PDB 长度不相等的结构，你可以尝试 struct.aln() 函数。

R-生物信息学秘籍-一-