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

布客飞龙
109 阅读23分钟

R 生物信息学秘籍（二）

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

译者：飞龙

协议：CC BY-NC-SA 4.0

第四章：系统发育分析与可视化

比较序列以推断进化关系是生物信息学中的一项基本技术。这项技术在 R 中也有着悠久的历史。除了 Bioconductor 之外，还有许多用于进化分析的包。在本章的食谱中，我们将详细研究如何处理来自不同来源的树形格式。重点将放在如何操作树形结构，聚焦于特定部分，并使用基于新的ggplot的树形可视化包进行可视化，这对于查看和注释大型树形结构特别有用。

本章将涵盖以下食谱：

  • 使用 ape 和 treeio 读取和写入各种树形格式

  • 使用 ggtree 快速可视化多基因树

  • 使用 treespace 量化树之间的距离

  • 使用 ape 提取和处理子树

  • 为对齐可视化创建点图

  • 使用 phangorn 从比对中重建树形

技术要求

你需要的示例数据可以从本书的 GitHub 仓库获得，链接是github.com/danmaclean/R_Bioinformatics_Cookbook。如果你想按照代码示例直接使用这些数据，请确保这些数据位于你工作目录的一个子目录中。

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

  • ape

  • adegraphics

  • Bioconductor:

    • Biostrings

    • ggtree

    • treeio

    • msa

  • devtools

  • dotplot

  • ggplot2

  • phangorn

  • treespace

Bioconductor 非常庞大，并且有自己的安装管理器。你可以使用以下代码进行安装：

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

更多信息请参见www.bioconductor.org/install/。

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

有时，在食谱的中途，我会中断代码，这样你就能看到一些中间输出或一个重要的对象结构，帮助理解。每当这种情况发生时，你会看到一个代码块，每行都以##（双哈希）符号开头。请看以下命令：

letters[1:5]

这将给我们带来以下输出：

## a b c d e

注意，输出行的前面会有##作为前缀。

使用 ape 和 treeio 读取和写入各种树形格式

系统发育分析是生物学和生物信息学的基石。相关程序种类繁多且复杂，计算耗时，数据集通常庞大。许多程序是独立运行的，且拥有专有的输入和输出格式。这形成了一个复杂的生态系统，我们在处理系统发育数据时必须应对，意味着通常最简单的策略是使用多种工具组合来加载、转换并保存分析结果，以便能够在不同的软件包中使用它们。在这个操作中，我们将探讨如何在 R 中处理系统发育树数据。迄今为止，R 对各种树格式的支持较为有限，但一些核心包具有足够标准化的对象，工作流可以集中在少数几种类型上，且转换为这些类型的过程已经简化。我们将使用ape和treeio包来导入和导出树数据。

准备工作

对于本节内容，我们需要从本书数据仓库中的datasets/ch4/文件夹获取树和系统发育信息，特别是mammal_tree.nwk和mammal_tree.nexus文件，它们分别是哺乳动物系统发育树的 Newick 格式和 Nexus 格式（你可以在本书的附录中查看这些文件格式的简要描述）。我们还需要beast_mcc.tree，这是 BEAST 运行的树文件，和RAxML_bipartitionsBranchLabels.H3，它是 RAxML 的输出文件。这两个文件来自treeio包提供的大量数据。我们还需要 Bioconductor 包treeio和ape包。

操作步骤...

使用ape和treeio读取和写入树格式的步骤如下：

  1. 加载ape库并加载树数据：
library(ape)
newick <-ape::read.tree(file.path(getwd(), "datasets", "ch4", "mammal_tree.nwk"))
nexus <-ape::read.nexus(file.path(getwd(), "datasets", "ch4", "mammal_tree.nexus"))
  1. 加载treeio库并加载 BEAST/RAxML 输出：
library(treeio)
beast <- read.beast(file.path(getwd(), "datasets", "ch4", "beast_mcc.tree"))
raxml <- read.raxml(file.path(getwd(), "datasets", "ch4", "RAxML_bipartitionsBranchLabels.H3"))
  1. 检查不同函数返回的对象类型：
class(newick)
class(nexus)

class(beast)
class(raxml)
  1. 将tidytree转换为phylo，反之亦然：
beast_phylo <- treeio::as.phylo(beast)
newick_tidytree <- treeio::as.treedata(newick)
  1. 使用以下代码编写输出文件：
treeio::write.beast(newick_tidytree,file = "mammal_tree.beast")
ape::write.nexus(beast_phylo, file = "beast_mcc.nexus")

工作原理...

在步骤 1中，我们使用了ape包中非常简单的加载函数——我们使用read.tree()和read.nexus()函数，它们能够读取通用格式的树。在步骤 2中，我们使用treeio包的特定格式函数来加载 BEAST 和 RaXML 的输出。步骤 3只是确认这些函数返回的对象类型；注意，ape返回phylo对象，而treeio返回treedata对象。通过treeio中的as.phylo()和as.treedata()，我们可以相互转换这两种对象类型。在步骤 4中，通过这种转换，我们可以将多种格式的输入导入到 R 中进行后续分析。最后，在步骤 5中，我们将文件写出。

参见

我们在步骤 2中使用的加载函数只是其中的几个。请参阅treeio包的说明文档，了解完整的函数列表。

使用 ggtree 快速可视化多个基因的树

一旦你计算出了树，首先要做的就是查看它。虽然许多程序都能做到这一点，但 R 具有一种极其强大、灵活且快速的系统，形式为ggtree包。在这个配方中，我们将学习如何将数据导入ggtree，并仅用几条命令重新布局、突出显示和注释树图。

准备就绪

你需要ggplot2、ggtree和ape包。此外，还需要从本书的仓库中的datasets/ch4文件夹中获取itol.nwk文件，它是来自Interactive Tree of Life在线工具的公共数据集的 Newick 格式树，包含 191 个物种。

如何实现...

使用ggtree快速可视化多基因的树，可以按照以下步骤执行：

  1. 加载库并获取 Newick 树的phylo对象：
library(ggplot2)
library(ggtree)
itol <-ape::read.tree(file.path(getwd(), "datasets", "ch4", "itol.nwk"))
  1. 制作基础的树形图：
ggtree(itol)
  1. 制作圆形图：
ggtree(itol, layout = "circular")
  1. 旋转并反转树：
ggtree(itol) + coord_flip() + scale_x_reverse()
  1. 向树的末端添加标签：
ggtree(itol) + geom_tiplab( color = "blue", size = 2)
  1. 制作颜色条来注释特定的谱系：
ggtree(itol, layout = "circular") + geom_strip(13,14, color="red", barsize = 1)
  1. 制作颜色块来突出显示特定的谱系：
ggtree(itol, layout = "unrooted") + geom_hilight_encircle(node = 11, fill = "steelblue")

它是如何工作的...

这段代码能够非常快速地实现很多功能。它通过类似于ggplot的层语法来实现这一点。以下是每个步骤的作用及其输出：

  1. 从文件加载一棵树。这里的树有 191 个叶子节点，所以非常大。恰好它是 Newick 格式的，因此我们使用ape的read.tree()函数。请注意，在后续步骤中，我们不需要为ggtree创建treedata对象；从read.tree()返回的phylo对象完全可以传递给ggtree()。

  2. 使用ggtree()创建一棵基本的树。这是一个封装了较长ggplot样式语法的函数，具体来说是ggplot(itol) + aes(x,y) + geom_tree() + theme_tree()。因此，可以将所有常见的ggplot函数作为额外层叠加到图表中。此步骤中的代码将生成以下图：

  1. 更改图表的布局。将布局参数设置为圆形将得到一个圆形的树。通过此参数，还可以选择其他多种树形：

  1. 我们可以使用标准的ggplot函数coord_flip()和scale_x_reverse()将树的左右方向改为上下方向，图表将呈现如下效果：

  1. 我们可以使用geom_tiplab()在树的末端添加名称。size参数设置文本的大小。以下代码将生成如下输出：

  1. 通过添加geom_strip()层，我们可以用一块颜色条注释树中的谱系。第一个参数（在本例中为13）是树中的起始节点，第二个参数是树中颜色条的结束节点。barsize参数设置颜色块的宽度。结果如下所示：

  1. 我们可以使用geom_hilight_encircle()几何图形在无根树中突出显示某些分支。我们需要为node参数选择一个值，这告诉ggtree()哪个节点是颜色的中心。这里的代码输出以下结果：

还有更多...

步骤 6和步骤 7依赖于我们知道要操作树中的哪些节点。因为节点是通过数字而非名称来标识的，这并不总是显而易见的。如果我们使用MRCA()（最近共同祖先）函数，就可以找到我们想要的节点编号。只需传递一个节点名称的向量，它就会返回代表 MRCA 的节点 ID：

MRCA(itol, tip=c("Photorhabdus_luminescens", "Blochmannia_floridanus"))

这将输出以下结果：

## 206

使用 treespace 量化树之间的差异

比较树以区分或将它们分组，可以帮助研究人员观察进化模式。通过跨物种或菌株跟踪单一基因的多棵树，可以揭示该基因在不同物种中变化的差异。这些方法的核心是树之间距离的度量。在这个食谱中，我们将计算一个这样的度量，找到 15 个不同物种中 20 棵不同基因树的成对差异——因此，每棵树中有 15 个相同名称的树枝。树之间的这种相似性通常是进行比较和获取距离所必需的，只有满足这些条件，我们才能进行这样的分析。

准备就绪

对于这个食谱，我们将使用treespace包来计算距离和聚类。我们将使用ape和adegraphics来加载附加功能和可视化函数。这里的输入数据将是datasets/ch4/gene_trees中的所有 20 个文件，每个文件都是一个 Newick 格式的树，代表 15 个物种中每个物种的单个基因。

如何进行...

使用treespace量化树之间的差异可以通过以下步骤执行：

  1. 加载库：
library(ape)
library(adegraphics)
library(treespace)
  1. 将所有树文件加载到一个multiPhylo对象中：
treefiles <- list.files(file.path(getwd(), "datasets", "ch4", "gene_trees"), full.names = TRUE) 
tree_list <- lapply(treefiles, read.tree) 
class(tree_list) <- "multiPhylo"
  1. 计算 Kendall-Colijn 距离：
comparisons <- treespace(tree_list, nf = 3)
  1. 绘制成对距离：
adegraphics::table.image(comparisons$D, nclass=25)
  1. 绘制主成分分析（PCA）和聚类：
plotGroves(comparisons$pco, lab.show=TRUE, lab.cex=1.5) 
groves <- findGroves(comparisons, nclust = 4) 
plotGroves(groves)

它是如何工作的...

这里简短而强大的代码非常有效——在少数几条命令中就能提供大量的分析。

在步骤 1中，首先加载我们需要的库。

在步骤 2中，加载必要的库后，我们创建一个字符向量treefiles，它保存我们希望使用的 20 棵树的路径。我们使用的list.files()函数接受一个文件系统路径作为参数，并返回该路径下找到的文件名。由于treefiles是一个向量，我们可以将其作为lapply()的第一个参数。

如果你不熟悉，lapply()是一个迭代函数，它返回一个 R 列表（因此是lapply()）。简单来说，lapply()会将第二个参数中指定的函数应用于第一个参数中的列表。当前元素作为目标函数的第一个参数传递。因此，在第 2 步中，我们在treefiles中列出的每个文件上运行ape的read.tree()函数，返回一个phylo树对象的列表。最后一步是确保tree_list对象具有类multiPhylo，这样我们就满足了后续函数的要求。值得一提的是，multiPhylo对象本身就是一个类似列表的对象，因此我们只需要通过class()函数将multiPhylo字符串添加到类属性中即可。

在第 3 步中，来自同名包的treespace()函数做了大量的分析。首先，它对输入中的所有树进行成对比较，然后使用 PCA 进行聚类。结果以列表对象返回，成员D包含树的成对距离，pco包含 PCA。默认的距离度量，Kendall-Colijn 距离，特别适用于我们这里的有根基因树，尽管该度量可以更改。nf参数仅告诉我们保留多少个主成分。由于我们的目标是绘图，我们不需要超过三个主成分。

在第 4 步中，我们使用adegraphics中的table.image()函数绘制comparisons$D中的距离矩阵—这是一个方便的热图风格函数。nclass参数告诉我们要使用多少个颜色层级。我们得到的图如下所示：

在第 5 步中，plotGroves()函数直接绘制treespace对象，因此我们可以看到 PCA 的绘图：

我们可以使用findGroves()函数将树分组为由nclust参数给定的组数，并重新绘制以查看结果：

还有更多...

如果你有很多树，并且图表看起来拥挤，你可以使用以下代码创建一个可以缩放和平移的交互式图：

plotGrovesD3(comparisons$pco, treeNames=paste0("species_", 1:10) )

使用 ape 提取和操作子树

在这个简短的教程中，我们将看看操作树是多么简单；具体来说，如何将子树提取为一个新对象，以及如何将树组合成其他树。

准备工作

我们需要一棵示例树；datasets/ch4文件夹中的mammal_tree.nwk文件就可以。我们需要的所有函数都可以在ape包中找到。

如何操作...

使用ape提取和操作子树可以通过以下步骤执行：

  1. 加载ape库，然后加载树：
library(ape)
newick <-read.tree(file.path(getwd(), "datasets", "ch4", "mammal_tree.nwk"))
  1. 获取所有子树的列表：
l <- subtrees(newick)
plot(newick)
plot(l[[4]], sub = "Node 4")
  1. 提取特定的子树：
small_tree <- extract.clade(newick, 9)
  1. 合并两棵树：
new_tree <- bind.tree(newick, small_tree, 3)
plot(new_tree)

它是如何工作的...

本食谱中的函数非常简单，但非常有用。

第一步是一个常见的树加载步骤。我们需要一个 phylo 对象的树来进行后续操作。

第二步使用subtrees()函数，它提取所有非平凡的（超过一个节点的）子树，并将它们放入一个列表中。列表中的成员根据原始树中的节点编号进行编号，每个列表中的对象都是一个phylo对象，像父树一样。我们可以使用plot()函数查看原始树和节点 4 的子树，并生成以下图示：

在第三步中，我们使用extract.clade()函数获取一个特定的子树。该函数的第一个参数是树，第二个参数是将被提取的节点。实际上，所有该节点下游的节点都会被提取，返回一个新的phylo对象。

最后的示例展示了如何使用bind.tree()函数将两个phylo对象结合起来。第一个参数是主树，它将接收第二个参数中的树。在这里，我们将把small_tree接到 Newick 树上。第三个参数是主树中要加入第二棵树的节点。同样，返回的是一个新的phylo对象。当我们绘制新的树时，可以看到相对于原始树的重复部分：

还有更多...

上述函数的一个小问题是它们要求我们知道要操作的节点编号。一种简单的访问方式是使用交互式的subtreeplot()命令。subtreeplot(newick)代码会生成一个交互式的树图，如这里的例子。通过点击树中的特定节点，我们可以让查看器渲染该子树并打印节点 ID。然后我们可以在函数中使用这个节点 ID：

创建对齐可视化的点图

对齐序列对的点图可能是最古老的对齐可视化方法。在这些图中，两个序列的位置分别绘制在x轴和y轴上，对于该坐标系中的每个坐标点，如果该点处的字母（核苷酸或氨基酸）相对应，就会画出一个点。由于该图能够展示出两个序列中不一定在同一区域匹配的区域，这是一种视觉上快速发现插入、删除以及结构重排的好方法。在这个例子中，我们将展示如何使用dotplot包和一些代码，快速构建一个点图，并获取文件中所有序列对的点图网格。

准备工作

我们将需要datasets/ch4/bhlh.fa文件，其中包含豌豆、大豆和莲花的三种基本螺旋-环-螺旋（bHLH）转录因子序列。我们还需要dotplot包，它不在 CRAN 或 Bioconductor 上，因此你需要使用devtools包从 GitHub 安装它。以下代码应该可行：

library(devtools)
install_github("evolvedmicrobe/dotplot", build_vignettes = FALSE)

如何实现...

创建用于对齐可视化的点图可以通过以下步骤完成：

  1. 加载库和序列：
library(Biostrings)
library(ggplot2)
library(dotplot)
seqs <- readAAStringSet(file.path(getwd(), "datasets", "ch4", "bhlh.fa"))
  1. 创建一个基本的点图：
dotPlotg(as.character(seqs[[1]]), as.character(seqs[[2]] ))
  1. 更改点图并应用 ggplot2 主题和标签：
dotPlotg(as.character(seqs[[1]]), as.character(seqs[[2]] ), wsize=7, wstep=5, nmatch=4) + 
theme_bw() + 
labs(x=names(seqs)[1], y=names(seqs)[2] )
  1. 创建一个函数，从提供的序列和序列索引中生成点图：
make_dot_plot <- function(i=1, j=1, seqs = NULL){
 seqi <- as.character(seqs[[i]])
 seqj <- as.character(seqs[[j]])
 namei <- names(seqs)[i]
 namej <- names(seqs)[j]
 return( dotPlotg(seqi, seqj ) + theme_bw() + labs(x=namei, y=namej) )
}
  1. 设置数据结构以运行函数：
combinations <- expand.grid(1:length(seqs),1:length(seqs))
plots <- vector("list", nrow(combinations) )
  1. 对所有可能的序列对组合运行该函数：
for (r in 1:nrow(combinations)){
 i <- combinations[r,]$Var1[[1]]
 j <- combinations[r,]$Var2[[1]]
 plots[[r]] <- make_dot_plot(i,j, seqs)
}
  1. 绘制点图网格：
cowplot::plot_grid(plotlist = plots)

它是如何工作的...

本食谱的第一部分非常熟悉。我们加载库，并使用 Biostrings 加载我们的蛋白质序列。请注意，seqs 变量中的序列是 XStringSet 类的一个实例。

在 步骤 2 中，我们可以使用 dotplotg() 函数创建一个基本的点图。参数是我们希望绘制的序列。注意，我们不能直接传递 XStringSet 对象；我们需要传递字符向量，因此我们使用 as.character() 函数将序列转换为该格式。运行此代码会生成以下点图：

在 步骤 3 中，我们通过首先改变匹配的考虑方式来详细说明基本点图。通过 wsize=7 选项，我们表示一次查看七个残基（而不是默认的一个），wstep=5 选项告诉绘图程序每次跳过五个残基（而不是默认的一个），nmatch=4 选项告诉绘图程序当四个残基相同时，标记该窗口为匹配。然后，我们通过添加 ggplot2 主题并以通常的 ggplot 方式进行自定义，最后使用标签函数添加轴名称。由此，我们得到了与第一个不同的点图：

自定义函数 make_dot_plot()，在 步骤 4 中定义，接收两个数字变量 i 和 j 以及 seqs 参数中的 XStringSet 对象。然后，它将 seqs 对象中的第 i 个和第 j 个序列转换为字符，并将其存储在 seqi 和 seqj 变量中。同时，它提取这些序列的名称并分别存储在 namei 和 namej 中。最后，它创建并返回一个使用这些变量生成的点图。

要使用该函数，我们需要两个条件：要绘制的序列组合和一个用于存储结果的列表。在 步骤 4 中，使用 expand.grid() 函数创建所有可能的序列组合的数据框，并将其存储在 combinations 变量中。使用 vector() 函数创建的 plots 变量包含一个 list 对象，该对象有足够的槽位来存储结果的点图。

步骤 6 是一个循环，遍历组合数据框中的每一行，提取我们希望处理的序列编号，并将其存储在 i 和 j 变量中。然后，调用 make_dot_plot() 函数，传入 i、j 和 seqs，并将其结果存储在我们创建的 plots 列表中。

最后，在步骤 7中，我们使用cowplot库的plot_grid()函数，结合我们的图形列表，生成所有可能组合的主图，如下所示：

使用 phangorn 从对齐中重建树

到目前为止，在本章节中，我们假设树已经可用并准备好使用。当然，构建系统树有很多方法，在本食谱中，我们将看看一些可用的不同方法。

准备工作

在本章节中，我们将使用datasets/ch4/文件、酵母 ABC 转运蛋白序列的abc.fa文件、Bioconductor Biostrings包以及来自 CRAN 的msa和phangorn包。

操作步骤...

使用phangorn构建树可以通过以下步骤执行：

  1. 加载库和序列，并进行对齐：
library(Biostrings)
library(msa)
library(phangorn)

seqs <- readAAStringSet(file.path(getwd(), "datasets", "ch4", "abc.fa"))
aln <- msa::msa(seqs, method=c("ClustalOmega"))
  1. 将对齐转换为phyDat对象：
aln <- as.phyDat(aln, type = "AA")
  1. 从距离矩阵生成 UPGMA 和邻接法树：
dist_mat <- dist.ml(aln)

upgma_tree <- upgma(dist_mat)
plot(upgma_tree, main="UPGMA")

nj_tree <- NJ(dist_mat)
plot(nj_tree,"unrooted", main="NJ")
  1. 计算自助法并绘制：
bootstraps <- bootstrap.phyDat(aln,FUN=function(x) { NJ(dist.ml(x)) } , bs=100)
plotBS(nj_tree, bootstraps, p = 10)

它是如何工作的...

第一步执行加载和氨基酸序列对齐，就像我们在之前的食谱中使用msa包时看到的那样，返回一个MsaAAMultipleAlignment对象。

第二步使用as.phyDat()函数将对齐转换为可以由phangorn函数使用的phyDat对象。

在步骤 3中，我们实际上生成树。树是由距离矩阵构建的，我们可以使用dist.ml()和我们的对齐来计算距离矩阵（这是一个最大似然距离度量；如果需要，也可以使用其他函数）。然后将dist_mat传递给upgma()和NJ()函数，分别生成 UPGMA 和邻接法树。这些函数返回标准的phylo对象，可以在许多其他函数中使用。这里，我们直接绘制：

在最后一步，我们使用bootstraps.phyDat()函数来计算树中各个分支的自助法支持。第一个参数是phyDat对象aln，而第二个参数FUN需要一个函数来计算树。这里，我们使用一个匿名函数，包装了我们最初用于生成nj_tree的NJ()方法。bs参数告诉函数计算多少次自助法。最后，我们可以使用plotBS()函数将结果图绘制到树上。

第五章：宏基因组学

高通量测序技术已经大大推动了宏基因组学的发展，从一个专注于研究单一序列变异的领域（例如 16S 核糖体 RNA (rRNA) 序列）发展到研究样本中可能存在的多个物种的整个基因组。识别物种或分类单元及其在样本中的丰度是一项计算挑战，要求生物信息学家处理序列准备、分类分配、分类比较以及定量分析。为此，许多专业实验室已经开发出相关的包，这些包创造了特定于宏基因组学中序列处理的新工具和新可视化。

在本章中，我们将查看一些在 R 中进行宏基因组学复杂分析的配方：

  • 使用phyloseq加载层级分类数据

  • 使用metacoder进行计数稀疏化以纠正样本差异

  • 使用dada2从原始读取中读取扩增子数据

  • 使用metacoder中的热图树可视化分类丰度

  • 使用vegan计算样本多样性

  • 将序列文件拆分为操作性分类单元

技术要求

你需要的示例数据可以从本书的 GitHub 库获取：github.com/PacktPublishing/R-Bioinformatics-Cookbook. 如果你想按原样使用代码示例，你需要确保这些数据位于你的工作目录的子目录中。

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

  • ape

  • Bioconductor

    • dada2

    • phyloseq

  • corrplot

  • cowplot

  • dplyr

  • kmer

  • magrittr

  • metacoder

  • RColorBrewer

  • vegan

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

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

进一步的信息可以通过以下链接获取：www.bioconductor.org/install/。

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

有时候，在一个配方的中间，我会中断代码，方便你查看一些中间输出或对象结构，这对于理解非常重要。每当这种情况发生时，你会看到一个代码块，其中每一行都以##（双井号）符号开头。请看下面的命令：

letters[1:5]

这将给我们以下输出：

## a b c d e

注意，输出行前面有##的前缀。

使用 phyloseq 加载层级分类数据

元基因组学管道通常从大型测序数据集开始，这些数据集会在强大且功能齐全的程序中处理，如 QIIME 和mothur。在这些情况下，我们希望将这些工具的结果整理成报告或进一步进行特定分析。在这个步骤中，我们将展示如何将 QIIME 和mothur的输出导入 R 中。

准备工作

对于这个简短的步骤，我们需要从Bioconductor安装phyloseq包，并从本书数据仓库的datasets/ch5文件夹中获取文件。

如何操作…

使用phyloseq加载层次分类数据可以通过以下步骤完成：

  1. 加载库：
library(phyloseq)
  1. 导入 QIIME 的.biom文件：
biom_file <- file.path(getwd(), "datasets", "ch5", "rich_sparse_otu_table.biom")
qiime <- import_biom(biom_file)
  1. 访问phyloseq对象的不同部分：
tax_table(qiime)
## Taxonomy Table:     [5 taxa by 7 taxonomic ranks]:
##          Rank1         Rank2               Rank3                   
## GG_OTU_1 "k__Bacteria" "p__Proteobacteria" "c__Gammaproteobacteria"
## GG_OTU_2 "k__Bacteria" "p__Cyanobacteria"  "c__Nostocophycideae"

otu_table(qiime)
## OTU Table:          [5 taxa and 6 samples]
##                      taxa are rows
##          Sample1 Sample2 Sample3 Sample4 Sample5 Sample6
## GG_OTU_1       0       0       1       0       0       0
## GG_OTU_2       5       1       0       2       3       1

sample_data(qiime)
##         BarcodeSequence  LinkerPrimerSequence BODY_SITE Description
## Sample1    CGCTTATCGAGA CATGCTGCCTCCCGTAGGAGT       gut   human gut
  1. 导入mothur数据文件：
mothur <- import_mothur(
  mothur_list_file = file.path(getwd(), "datasets", "ch5", "esophagus.fn.list"),
  mothur_group_file = file.path(getwd(), "datasets", "ch5", "esophagus.good.groups"),
  mothur_tree_file = file.path(getwd(), "datasets", "ch5", "esophagus.tree")
)
  1. 访问phyloseq对象中的otu对象：
otu_table(mothur)
## OTU Table:          [591 taxa and 3 samples]
##                      taxa are rows
##          B C  D
## 9_6_14   2 0  0
## 9_6_25   1 0  0

它是如何工作的…

在这个简单的步骤中，我们创建对象并使用访问器函数。

在第 1 步中，我们按惯例加载phyloseq库。

然后，在第 2 步中，我们定义一个文件，并将其作为import_biom()函数的第一个参数。该函数可以读取来自 QIIME 的现代biom格式输出，支持未压缩的 JSON 和压缩的 HDF5 格式。类型会自动检测。我们将得到一个完全填充的phyloseq对象。

在第 3 步中，我们使用访问器函数获取对象的子部分，使用tax_table()获取分类数据，使用otu_table()获取 OTU，使用sample_data()获取样本数据；这些都可以作为类似矩阵的对象在后续处理中轻松使用。

在第 4 步中，我们改变方向，使用mothur的输出数据。我们至少需要一个列表文件和一个分组文件，这些文件路径通过mothur_list_file和mothur_group_file参数指定。这里，我们还通过mothur_tree_file参数指定一个 Newick 格式的树。

同样，我们可以使用phyloseq的访问器函数otu_table()来获取 OTU。对于最小的mothur数据，我们在此指定不能获取样本数据或分类学表。

还有更多…

如果你有来自旧版 QIIME 的专有格式数据，可以使用import_qiime()函数。若你附加了树对象，也有一个访问器函数——phy_tree()。

另请参见

QIIME 和mothur程序的官网和 wiki 页面非常出色地展示了如何在 R 中处理它们管道输出的数据，尤其是mothur。如果你想要一些数据分析的思路，可以试试它们。

使用 metacoder 稀释计数并校正样本差异。

在宏基因组学中，一个常见的问题是询问某个样本中存在哪些物种，以及两个或多个样本之间的差异。由于样本可能由不同数量的观察值组成——在宏基因组学中，这意味着生成的不同数量的读取——因此样本的分类丰富度会随着测序深度的增加而增加。为了公平地评估样本中不同类群的多样性，宏基因组学家通常会对计数进行稀释，以确保样本的测序深度相同。本质上，这意味着将样本深度减少到最小的样本深度。在这个步骤中，我们将对来自 biom 文件的 OTU 计数进行稀释。

准备工作

对于这个步骤，你将需要 metacoder 包和 datasets/ch5/rich_high_count_otu.biom，这是一个包含六个样本（标记为 Sample1–Sample6）和五个 OTU 的示例 biom 文件。当然，这是一个非常小的文件，只用于学习代码的工作原理。真实的宏基因组数据集要大得多。

如何操作...

使用 metacoder 对计数进行稀释并修正样本差异可以通过以下步骤完成：

  1. 加载库和文件：
library(metacoder)
biom_file <- file.path(getwd(), "datasets", "ch5", "rich_high_count_otu.biom")
taxdata <- parse_qiime_biom(biom_file)
  1. 创建样本中计数的直方图：
sample_ids <- paste0("Sample", 1:6)
hist_data <- colSums(taxdata$data$otu_table[, sample_ids])
hist(hist_data, prob= TRUE, breaks=3)
lines(density(hist_data, adjust = 2), col="red")
  1. 调用稀释函数并过滤掉可能生成的低 OTU：
taxdata$data$rarefied_otus <- rarefy_obs(taxdata, "otu_table", other_cols = TRUE)
low_otu_index <- rowSums(taxdata$data$rarefied_otus[, sample_ids]) <=20
taxdata <- filter_obs(taxdata, "rarefied_otus", ! low_otu_index)
taxdata$data$rarefied_otus

工作原理...

这里的总体模式是加载文件，检查样本 OTU 计数的分布，并应用稀释法。

第一步是加载库并导入示例文件。我们通过准备 rich_high_count_otu.biom 文件来完成此操作，并将其传递给 parse_qiime() 函数。这个 metacoder 函数只是读取生物群落文件并返回一个 taxmap 对象（另一种用于保存分类数据的对象），我们可以在 metacoder 函数中使用它。

接下来，我们希望检查样本 OTU 计数的分布，这可以通过准备一个直方图来完成。我们使用 paste() 函数创建一个包含样本名称的字符向量，并用它来通过命名索引提取 otu_table 中的计数列。这些列的子集被传递给 colSums() 函数，后者会获取 hist_data 向量中每个样本的总计数。现在，我们可以使用 hist() 创建这些计数的直方图，并通过 lines() 和 density() 函数在 hist_data 上添加密度曲线。请注意，结果图（在以下直方图中）看起来较为稀疏，因为示例文件中的样本数量较少。这里的最低计数为我们提供了最低测序样本的概念。如果存在明显较低的样本，最好先去除这些列：

现在，我们可以进行稀释分析。我们在 taxdata 上使用 rarefy_obs() 函数；第二个参数（值为 "otu_table"）是 taxdata 对象中包含 OTU 计数的槽位名称。由于稀释会减少计数，因此我们现在需要去除在所有样本中已经减少过多的计数。因此，我们使用 rowSums() 函数，并通过样本名称索引 taxdata$data$rarefied_otus 对象，得到一个逻辑向量，表示哪些 OTU 的总计数低于 20。最后，我们在 taxdata 上使用 filter_obs() 函数；第二个参数（值为 "rarefied_otus"）是 taxdata 对象中包含稀释 OTU 计数的槽位名称。! 字符用于反转低计数 OTU 的逻辑向量，因为 filter_obs() 会保留通过的行，而我们希望移除这些行。最终输出的是一个稀释后的 OTU 计数集。

请注意，在以下输出中，OTU 行 3 已通过低计数被移除：

## # A tibble: 4 x 8
##   taxon_id otu_id   Sample1 Sample2 Sample3 Sample4 Sample5 Sample6
##   <chr>    <chr>      <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
## 1 ax       GG_OTU_1      24    1004     847    1979    1070    1170
## 2 ay       GG_OTU_2     872       0     704     500    1013     689
## 3 ba       GG_OTU_4     875    1144    1211     217       0    1180
## 4 ax       GG_OTU_5    1270     893     276     338     953       0

还有更多……

我们可以通过稀释曲线估算一个有用的计数水平。通过这些曲线，计数会在不同的样本大小下随机抽取，并计算 OTU 中的物种数。当物种数量停止增加时，我们就知道我们已经有足够的读取数据，并且不会从处理额外的计数中获得更多价值。rarecurve() 函数在 vegan 包中可以帮助我们完成这一任务。我们提供一个 OTU 表（请注意，该函数需要样本按行排列，因此我们必须使用 t() 函数旋转我们的 taxdata OTU 表）。然后，我们为 sample 参数传入最小样本大小。我们使用 colSums() 和 min() 函数来获取该最小样本 OTU 计数。输出结果如下图所示：

在这里，我们可以清楚地看到，超过 20,000 的样本并没有增加物种的丰富度。

使用 dada2 从原始读取数据中读取扩增子数据

作为元基因组学中的一项长期技术，特别是对于有兴趣进行细菌微生物组研究的人，使用 16S 或 18S rRNA 基因的克隆拷贝（扩增子）测序来创建物种谱图。这些方法可以利用低通量测序以及目标序列的知识来对每个克隆序列进行分类，从而简化了将分类单元分配到读取序列的复杂任务。在这个方案中，我们将使用 dada2 包从原始 fastq 序列读取中运行扩增子分析。我们将执行质量控制和 OTU 分配步骤，并使用一种有趣的机器学习方法来对序列进行分类。

准备工作

对于这个食谱，我们需要 Bioconductor 的dada2包和 CRAN 的cowplot包。我们将使用来自短读档案实验SRR9040914的某些宏基因组序列数据，其中检测了一个旅游中心潮汐海湖水中的物种组成，因为游客向湖里投掷硬币并许愿。我们将使用二十个fastq文件，每个文件包含 2500 条数据，每个文件经过压缩，且可以在本书的仓库datasets/ch5/fq中找到。这是 Illumina 读取的一个小子集。我们还需要datasets/ch5/rdp_train_set_14.fa文件，它是由dada团队维护的 16S 序列训练集之一。更多训练集请参见benjjneb.github.io/dada2/training.html。

如何进行...

使用dada2从原始读取数据中读取扩增子数据，可以按照以下步骤进行：

  1. 加载库并为每个fastq文件准备绘图：
library(dada2)
library(cowplot)

fq_dir <- file.path(getwd(), "datasets", "ch5", "fq")
read_files <- list.files(fq_dir, full.names = TRUE,  pattern = "fq.gz")

quality_plots <- lapply(read_files, plotQualityProfile)
plot_grid(plotlist = quality_plots)
  1. 对文件进行质量修剪和去重复：
for (fq in read_files ){
  out_fq <- paste0(fq, ".trimmed.filtered")
  fastqFilter(fq, out_fq, trimLeft=10, truncLen=250,
                      maxN=0, maxEE=2, truncQ=2,
                      compress=TRUE)
}

trimmed_files <-  list.files(fq_dir, full.names = TRUE, pattern = "trimmed.filtered")
derep_reads <- derepFastq(trimmed_files)
  1. 从样本子集估算dada2模型：
trimmed_files <-  list.files(fq_dir, full.names = TRUE, pattern = "trimmed.filtered")
derep_reads <- derepFastq(trimmed_files)

dd_model <- dada(derep_reads[1:5], err=NULL, selfConsist=TRUE)
  1. 使用第 3 步中估算的参数推测样本的序列组成：
dada_all <- dada(derep_reads, err=dd_model[[1]]$err_out, pool=TRUE)
  1. 给序列分配分类：
sequence_tb <-makeSequenceTable( dada_all )
taxonomy_tb <- assignTaxonomy(sequence_tb, refFasta = file.path(getwd(), "datasets", "ch5", "rdp_train_set_14.fa")) 
taxonomy_tb[1, 1:6]

它是如何工作的...

我们首先通过将包含fastq目录的fq_dir变量传递给list.files()函数，创建一个所有fastq文件路径的向量。然后，我们使用循环函数lapply()，遍历每个fastq文件路径，并依次运行dada函数plotQualityProfile()。每个结果绘图对象都会保存到列表对象quality_plots中。当将绘图列表传递给plotlist参数时，cowplot函数plot_grid()将把这些绘图以网格形式显示出来。

我们得到下图所示的绘图。请注意，fastq质量分数在前 10 个核苷酸左右较差，且在大约 260 个核苷酸后出现问题。这些将是下一步的修剪点：

为了进行修剪，我们对 read_files 中的 fastq 文件运行一个循环。在每次循环迭代中，我们通过将文本 "trimmed.filtered" 附加到文件名上来创建一个输出的 fastq 文件名 out_fq（因为我们将修剪后的读取结果保存到新文件中，而不是内存中），然后运行 fastqFilter() 修剪函数，传递输入文件名 fq、输出文件名 out_fq 和修剪参数。循环结束后，我们将得到一个包含修剪过的读取文件的文件夹。通过再次使用 list.files() 函数加载这些文件的名称，这次只匹配文件名中带有 "trimmed.filtered" 模式的文件。所有这些文件都被加载到内存中，并使用 derepFaistq() 函数进行去重。接着，我们通过对部分文件使用 dada() 函数，计算组成推断步骤的参数。我们通过索引 derep_reads 对象，传递前五组去重后的文件。通过将 err 设置为 NULL 和 selfConsist 设置为 TRUE，我们强制 dada() 从数据中估算参数，并将结果保存在 dd_model 变量中。

接下来，我们对所有数据运行 dada() 函数，将 err 参数设置为先前估算并存储在 dd_model 中的值。此步骤计算整个数据集的最终序列组成。

最后，我们可以使用 dada() 函数的结果创建序列表，并利用该表通过 assignTaxonomy() 查找 OTU。此函数使用朴素贝叶斯分类器将序列分配到分类群中，基于提供的 rdp_train_set_14.fa 文件中的训练集分类。该函数的输出是每个序列的分类。结果表格 taxonomy_tb 中的一行如下所示：

## Kingdom Phylum 
## "Bacteria" "Cyanobacteria/Chloroplast" 
## Class Order 
## "Chloroplast" "Chloroplast" 
## Family Genus 
## "Bacillariophyta" NA

另见

本配方中使用的函数 fastqFilter() 和 derepFastQ() 也有用于配对读取的变体。

使用 metacoder 可视化分类丰度的热树

无论我们如何估算分类丰度，通常都需要创建一种可视化方法，以便在单一图形中总结数据的广泛趋势。一个表达力强且易于解读的可视化方式是热树。这些是对感兴趣的分类群的系统发育树的呈现，其中数据被映射到图形的视觉元素上。例如，一个分类群被观察到的次数可能通过改变树枝的颜色或粗细来表示。不同的数据集可以通过比较每个数据集的树形图来轻松发现差异。在本配方中，我们将构建一棵热树并进行自定义。

准备工作

我们需要输入 .biom 文件 datasets/ch5/rich_high_count_otu.biom 和 metacoder、RColorBrewer 包。

如何实现...

使用 metacoder 可视化分类丰度的热树可以通过以下步骤完成：

  1. 加载库和输入文件：
library(metacoder)
library(RColorBrewer)
biom_file <- file.path(getwd(), "datasets", "ch5", "rich_high_count_otu.biom")
taxdata <- parse_qiime_biom(biom_file)
  1. 将自定义选项传递给树绘制函数：
heat_tree(taxdata, 
    node_label = taxon_names, 
    node_size = n_obs, 
    node_color = n_supertaxa, 
    layout = "gem", 
    title = "sample heat tree", 
    node_color_axis_label = "Number of Supertaxa", 
    node_size_axis_label = "Number of OTUs", 
    node_color_range = RColorBrewer::brewer.pal(5, "Greens")
)

工作原理...

首先，我们加载库并使用parse_qiime_biom()函数从biom文件获取一个metacoder taxmap对象。

然后我们使用heat_tree()函数来渲染树。只需要传递taxdata taxmap对象——这将生成默认的树——其他的参数则是用来定制树的。node_label指定在taxdata对象中用于节点标签的列；在这里，我们使用taxon_names，特别注意这里没有加引号，因为该函数使用非标准评估方式，类似于你在tidyverse和ggplot函数中可能已经熟悉的方式。node_size根据给定的列控制节点大小。在这里，n_obs和node_color提供了影响节点颜色变化的参数（注意，这不是颜色的集合——而是应该被相同/不同颜色标记的内容）。接下来，layout参数告诉函数如何在渲染中展开树的分支。在接下来的三个参数标题中，node_color_axis和node_size_axis_label仅仅是图形的标签。最后，node_color_range获取一个颜色标识符向量，用来绘制图形。这里，我们使用RColorBrewer包的函数brewer.pal()，它返回这样的内容。它的第一个参数是要返回的颜色数量，第二个参数是要选择的调色板名称。设置好所有这些后，我们从我们的输入小文件中得到如下图：

使用 vegan 计算样本多样性

在生态学和宏基因组学研究中，一个常见的任务是估算样本内或样本之间的物种（或分类学）多样性，以查看哪些样本多样性较高或较低。有多种度量方法可以衡量样本内外的多样性，包括辛普森指数和布雷指数。在这个示例中，我们将查看能够从常见的 OTU 表格中返回多样性度量的函数。

准备工作

我们需要样本的.biom输入文件，datasets/ch5/rich_high_count_otu.biom，以及vegan包。

如何实现……

使用vegan计算样本多样性可以通过以下步骤完成：

  1. 加载库并从样本文件准备 OTU 表格：
library(vegan)
biom_file <- file.path(getwd(), "datasets", "ch5", "rich_high_count_otu.biom")
taxdata <- metacoder::parse_qiime_biom(biom_file)
otu_table <- taxdata$data$otu_table[, paste0("Sample", 1:6)]
  1. 计算α多样性：
alpha_diversity <- diversity(otu_table, MARGIN=2, index="simpson")
barplot(alpha_diversity)
  1. 计算β多样性：
between_sample <- vegdist(t(otu_table), index = "bray")

between_sample_m <- as.matrix(between_sample, ncol = 6)
corrplot::corrplot(between_sample_m, method="circle", type = "upper", diag = FALSE )

它是如何工作的……

第一步非常简单。在这里，我们使用metacoder parse_qiime_biom()函数加载我们的biom文件，然后对生成的taxdata$data$otu_table插槽进行子集化，提取一个简单的 OTU 表格到otu_table中。

现在我们可以调用vegan包中的diversity()函数。index参数设置为"simpson"，所以函数将使用辛普森指数来计算样本内的多样性。MARGIN参数告诉函数样本是按行还是按列排列：1 表示按行，2 表示按列。diversity()函数返回一个命名的向量，便于使用barplot()函数进行可视化，生成如下图：

现在我们可以使用vegdist()函数运行样本间的多样性度量；同样，index参数设置要使用的指数，这里我们选择 Bray 指数。vegdist()期望样本数据是按行排列的，所以我们使用t()函数来旋转otu_table。生成的对象存储在between_sample中——它是一个成对相关性表格，我们可以在corrplot中可视化它。为了做到这一点，我们需要通过as.matrix()将其转换为矩阵；ncol参数应与样本的数量匹配，以便为每个样本生成一列。返回的矩阵between_sample_m可以传递给corrplot()函数进行渲染。通过将method设置为circle，type设置为upper，并将diag设置为false，我们可以得到一个只显示矩阵上三角部分的图，而没有自我与自我的比较，从而减少图中的冗余。

输出结果如下：

另见...

本示例中的相关性图明确显示了几个样本的相关性，但在非常大的实验中可能变得不易处理。在此阶段，您可能需要考虑 PCA 或其他某种多维尺度方法。

将序列文件拆分为 OTU

对于经过清理和修剪的测序数据，最常见的任务之一是将序列划分为 OTU。在这方面有很多方法；在本示例中，我们将探讨一种将序列拆分为指定长度的子序列并对其执行某种层次聚类的方式，以创建组。

准备工作

这里的关键包是kmer包，我们将使用datasets/ch5/fq文件夹中的一个示例fastq序列文件。我们还将使用dplyr和magrittr包以提高便利性。

如何做...

将序列文件拆分为 OTU 可以通过以下步骤完成：

  1. 加载数据并计算 OTU：
library(kmer)
library(magrittr)
library(ape)
seqs <- ape::read.fastq(file.path(getwd(), "datasets", "ch5","fq", "SRR9040914ab.fq.gz")
otu_vec <- otu(seqs, k = 6, threshold = 0.99 )
  1. 计算每个 OTU 簇的频率：
data.frame(
    seqid = names(otu_vec), 
    cluster = otu_vec, 
    row.names = NULL) %>% 
dplyr::group_by(cluster) %>% 
dplyr::summarize(count = dplyr::n() )

它是如何工作的...

加载库后，我们使用ape包中的read.fastq()函数获取表示序列的DNAbin对象。kmer包中的关键函数otu()可以直接使用DNAbin seqs对象创建长度为k的 k-mer，并对其执行层次聚类。threshold参数设置 OTU 的身份阈值。该函数返回一个命名向量，其中名称为序列 ID，每个值为该序列所属簇的 ID。

然后我们可以使用otu_vec通过data.frame构建一个中间数据框，使用names属性设置一个seqid列，并将聚类成员信息放入一个名为cluster的列中。通过将row.names设置为NULL，我们去掉行名称。接着，我们使用magrittr管道符号%>%，通过dplyr::group()按聚类对数据框进行分组，并通过dplyr::summarize()创建一个汇总数据框。通过将计数设置为dplyr::n()函数的结果，我们可以得到每个聚类在命名向量中出现的次数——即每个聚类中分配到的读取次数。

第六章：从光谱到注释的蛋白质组学

质谱（MS）数据通常包含必须经过生物信息学处理的光谱，以识别候选肽。这些肽包括分配，计数可以使用各种技术和包进行分析。用于蛋白质组学的各种图形用户界面驱动工具意味着出现了多种文件格式，最初可能很难处理。这些配方将探索如何利用RforProteomics项目中的优秀解析器和格式转换工具来分析和验证光谱，甚至向你展示如何在基因组浏览器中查看你的肽以及基因模型等其他基因组信息。

本章将涵盖以下配方：

  • 以视觉方式表示原始 MS 数据

  • 在基因组浏览器中查看蛋白质组学数据

  • 可视化肽命中次数的分布以寻找阈值

  • 转换 MS 格式以在工具之间传输数据

  • 使用 protViz 将光谱与肽匹配进行验证

  • 应用质量控制滤波器到光谱

  • 识别与肽匹配的基因组位置

技术要求

你将需要的样本数据可以从本书的 GitHub 仓库中获取，地址是github.com/danmaclean/R_Bioinformatics_Cookbook[.] 如果你想按原样使用代码示例，那么你需要确保该数据位于工作目录的子目录中。

以下是你将需要的 R 包。通常，你可以通过install.packages("package_name")来安装这些包。列在Bioconductor下的包需要通过专用的安装程序进行安装，如此处所述。如果你需要做其他事情，安装方法将在这些包所使用的配方中描述：

  • Bioconductor

    • EnsDb.Hsapiens.v86

    • MSnID

    • MSnbase

    • mzR

    • proteoQC

    • rtracklayer

  • data.table

  • dplyr

  • ggplot2

  • protViz

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

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

然后，你可以使用此代码安装软件包：

BiocManager::install("package_name")

更多信息可以在www.bioconductor.org/install/找到。

在 R 中，用户通常会加载一个库并直接通过名称使用其中的函数。这在交互式会话中非常方便，但在加载多个包时可能会造成困惑。为了明确在某一时刻我正在使用哪个包和函数，我偶尔会使用packageName::functionName()这种惯例。

有时，在配方中间，我会中断代码，以便你查看一些中间输出或对你理解非常重要的对象结构。每当发生这种情况时，你会看到一个代码块，其中每行以##（即双哈希符号）开头。考虑以下命令：

letters[1:5]

这将产生以下输出：

## a b c d e

请注意，输出行以##为前缀。

以可视化方式表示原始质谱数据

蛋白质组学分析的原始数据就是由质谱仪生成的光谱。每种类型的质谱仪都有不同的本地文件格式来编码光谱。检查和分析光谱从加载文件并将其强制转换为通用的对象类型开始。在这个示例中，我们将学习如何加载不同类型的文件，查看元数据，并绘制光谱本身。

准备就绪

对于这个示例，我们需要使用Bioconductor包、mzR包以及本书数据仓库中的一些文件，这些文件位于datasets/ch6文件夹中。我们将使用三个不同的文件，它们的选择并不是基于文件中的数据，而是因为它们分别代表了最常见的三种质谱文件类型：mzXML、mzdata和mzML。这些示例文件都来自mzdata包。由于它们是提取出来的，你不需要安装该包，但如果你需要更多的示例文件，这个包是一个不错的选择。

如何做到...

原始质谱数据可以通过以下步骤以可视化方式表示：

  1. 加载库：
library(mzR)
library(MSnbase)
  1. 将文件加载到对象中：
mzxml_file <- file.path(getwd(), "datasets", "ch6", "threonine_i2_e35_pH_tree.mzXML")
ms1 <- openMSfile(mzxml_file)

mzdata_file <- file.path(getwd(), "datasets", "ch6", "HAM004_641fE_14-11-07--Exp1.extracted.mzdata")
ms2 <- openMSfile(mzdata_file)

mzml_file <-  file.path(getwd(), "datasets", "ch6", "MM8.mzML")
ms3 <- openMSfile(mzml_file)
  1. 查看可用的元数据：
runInfo(ms3)

## $scanCount
## [1] 198
## 
## $lowMz
## [1] 95.51765
## 
## $highMz
## [1] 1005.043
## 
## $dStartTime
## [1] 0.486
## 
## $dEndTime
## [1] 66.7818
## 
## $msLevels
## [1] 1
## 
## $startTimeStamp
## [1] "2008-09-01T09:48:37.296+01:00"

sampleInfo(ms1)

## [1] ""
  1. 绘制光谱：
msn_exp <- MSnbase::readMSData(mzxml_file)
MSnbase::plot(msn_exp, full = TRUE)
MSnbase::plot(msn_exp[5], full = TRUE)

它是如何工作的...

在步骤 1中，我们加载了所需的库。主要的库是mzR。

在步骤 2中，我们使用系统无关的file.path()函数定义了我们将要加载的文件路径，该函数返回一个包含文件名的字符向量。然后，我们使用该文件名通过mzR中的openMSfile()函数来创建一个代表相应文件中数据的mzR对象。请注意，我们本质上运行了相同的代码三次，只是每次更改文件和输入文件类型。openMSfile()函数会自动检测文件的格式。

在步骤 3中，我们使用mzR包的访问器函数runInfo()和sampleInfo()来提取输入文件中的一些元数据。请注意，sampleInfo()与ms1一起使用时没有返回任何内容——这是因为该特定文件中没有包含这些数据。可以返回的元数据取决于文件和文件类型。

在步骤 4中，我们使用MSnbase包的readMSData()函数加载文件。该函数在后台使用mzR，因此可以执行相同的操作，但它返回一个修改后的MSnbase类对象。这意味着一些通用的绘图函数将会生效。接着，我们使用plot()函数创建文件中所有光谱的图像：

然后，通过使用索引，我们创建了仅包含文件中第五个光谱的图像：

在基因组浏览器中查看蛋白质组学数据

一旦我们获取了质谱数据，并使用如 Xtandem、MSGF+或 Mascot 等搜索引擎软件识别了光谱描述的肽段和蛋白质，我们可能希望在其基因组上下文中查看这些数据，并与其他重要数据一起展示。在本食谱中，我们将展示如何从搜索文件中提取肽段和 Uniprot ID，找到这些 Uniprot ID 映射的基因，并创建一个显示这些基因的基因组浏览器轨道。可以将这些发送到 UCSC 人类基因组浏览器，互动网页将会自动在本地浏览器中加载。

准备工作：

对于这个食谱，您需要安装 Bioconductor 包MSnID、EnsDB.Hsapiens.v86和rtracklayer，以及本书仓库中datasets/ch6文件夹下的HeLa_180123_m43_r2_CAM.mzid.gz文件。为了使这个食谱生效，您还需要连接到互联网，并拥有一个能够运行 UCSC 基因组浏览器的现代浏览器，网址是genome.ucsc.edu。

如何操作...

通过以下步骤可以在基因组浏览器中查看蛋白质组学数据：

  1. 加载库：
library(MSnID)
library(EnsDb.Hsapiens.v86)
library(rtracklayer)
  1. 创建并填充搜索文件对象：
msnid <- MSnID() 
msnid <- read_mzIDs(msnid, file.path(getwd(), "datasets", "ch6", "HeLa_180123_m43_r2_CAM.mzid.gz"))
  1. 提取包含有用命中的行和包含有用信息的列：
real_hits <- msnid@psms[! msnid@psms$isDecoy, ]
required_info <- real_hits[, c("spectrumID", "pepSeq", "accession", "start", "end")]
  1. 从accession列提取 Uniprot ID：
uniprot_ids <- unlist(lapply(strsplit(required_info$accession, "\\|"), function(x){x[2]}) )
uniprot_ids <- uniprot_ids[!is.na(uniprot_ids)]
  1. 创建数据库连接并获取与我们的 Uniprot ID 匹配的基因：
edb <- EnsDb.Hsapiens.v86
genes_for_prots <- genes(edb, 
    filter = UniprotFilter(uniprot_ids), 
    columns = c("gene_name", "gene_seq_start", "gene_seq_end", "seq_name"))
  1. 设置基因组浏览器轨道：
track <- GRangesForUCSCGenome("hg38", 
    paste0("chr",seqnames(genes_for_prots)), 
    ranges(genes_for_prots), 
    strand(genes_for_prots), 
    genes_for_prots$gene_name, 
    genes_for_prots$uniprot_id )
  1. 设置浏览器会话并查看：
session <- browserSession("UCSC")
track(session, "my_peptides") <- track

first_peptide <- track[1] 
view <- browserView(session, first_peptide * -5, pack = "my_peptides")

工作原理...

步骤 1是我们的标准库加载步骤。

步骤 2是数据加载步骤。这一步有点特殊。我们不仅仅调用一个文件读取函数，而是首先创建并清空MSnID对象，并将数据加载到其中。我们使用MSnID()函数创建msnid，然后将其传递给read_mzid()函数，实际上将数据加载到其中。

步骤 3关注于从msnid对象中提取我们关心的信息。我们需要的是匹配实际命中的行，而不是诱饵行，因此我们直接访问msnid@psms槽，该槽包含有用数据，并筛选出isDecoy值为FALSE的行。这将给我们一个对象，我们将其保存在real_hits变量中。接下来，我们使用real_hits从原始对象中的许多列中选择一些有用的列。

步骤 4 帮助我们提取嵌入在 accession 列字段中的 Uniprot ID。需要注意的是，这些值来自于搜索引擎数据库中使用的名称。自然地，这一步会根据数据库的精确格式有所不同，但一般的模式适用。我们有一组嵌套得相当密集的函数，其解析过程如下：内层的匿名函数function(x){x[2]}返回它所传入的任何向量的第二个元素。我们使用lapply()将这个函数应用于从 strsplit() 函数返回的 accession 列中的每个元素。最后，由于lapply()返回的是列表，我们使用unlist()将其展开成我们所需的向量。有时，这会生成 NAs，因为某些 Uniprot ID 可能不存在，所以我们通过子集和 is.na() 将其从向量中删除。

在 步骤 5 中，我们连接到 Ensembl 数据库包，并使用 genes() 函数获取与我们的 Uniprot ID 匹配的 Ensembl 基因。Uniprot ID 的向量被传递到 UniprotFilter() 函数中，并且通过 columns 参数，我们选择从数据库中返回的数据。这为我们提供了一个 GRanges 对象，包含了构建浏览器轨道所需的所有信息。

在 步骤 6 中，我们使用辅助函数GRangesForUCSCGenome()，并传入我们希望查看的基因组版本——hg38，然后是基本的染色体名称、坐标和链信息，这是创建 GRanges 对象所需的数据。我们可以使用 seqnames()、ranges() 和 strand() 访问函数，从我们之前创建的 genes_for_prots 对象中提取这些信息。UCSC 中的序列名称以 chr 为前缀，因此我们使用 paste 将其添加到我们的序列名称数据中。我们还为基因名称和基因 ID 创建了列，以便在最终的视图中保留这些信息。我们将结果保存在 track 变量中。

最后，在 步骤 7 中，我们可以渲染我们创建的轨道。首先，我们创建一个表示 UCSC 会话的会话对象，并分别使用 session() 和 track() 函数将轨道添加到其中。我们通过将第一个肽传递给 view() 函数来选择关注的肽，view() 函数会实际打开一个新的网页浏览器窗口，显示请求的数据。view() 的第二个参数指定缩放级别，通过将参数公式化为 first_peptide * -5，我们可以获得一个可以容纳五个请求特征的缩放视图。

在写作时，这个配方生成了如下视图。请注意，最顶部的轨道是我们的 my_peptides 轨道：

还有更多内容……

你可能注意到，这个示例实际上绘制的是整个基因，而不是我们最初开始时的肽段命中数据。绘制基因是最简单的情况，但要绘制肽段只需要稍作修改。在步骤 5中，我们创建了一个对象genes_for_prots，它给出了基因的起始和结束位置。早期的msnid@psms对象包含这些基因内肽段的起始和结束位置，这些位置是从命中开始处索引的，因此通过将两个数据合并在一起，就可以创建一个代表肽段而非基因的对象。

对于那些没有使用 UCSC 浏览器中的生物体的用户，仍然可以生成命中的 GFF 文件，并上传到其他基因组浏览器——许多浏览器都提供这种功能。只需在步骤 5结束时停止该示例，并使用rtracklayer::export()函数创建一个 GFF 文件。

可视化肽段命中数分布以找到阈值

每个质谱实验都需要了解哪些肽段命中数代表噪声或异常特征，例如在蛋白质组中出现过度表达的肽段。在这个示例中，我们将使用tidyverse工具，如dplyr和ggplot，结合一些巧妙的可视化技巧，创建图形，帮助你了解质谱实验中肽段命中的分布和限制。

准备工作

对于这个示例，你将需要MSnId、data.table、dplyr和ggplot包。我们将使用来自本书仓库datasets/ch6文件夹中的mzid文件HeLa_180123_m43_r2_CAM.mzid.gz。

如何实现...

可视化肽段命中数分布以找到阈值，可以使用以下步骤完成：

  1. 加载库和数据：
library(MSnID)
library(data.table)
library(dplyr)
library(ggplot2)
msnid <- MSnID() 
msnid <- read_mzIDs(msnid, file.path(getwd(), "datasets", "ch6", "HeLa_180123_m43_r2_CAM.mzid.gz")) 
peptide_info <- as(msnid, "data.table")
  1. 过滤掉虚假数据行，并统计每个肽段出现的次数：
per_peptide_counts <- peptide_info %>% 
 filter(isDecoy == FALSE) %>%
 group_by(pepSeq) %>%
 summarise(count = n() ) %>% 
 mutate(sample = rep("peptide_counts", length(counts) ) )
  1. 创建一个小提琴图和抖动图，显示命中计数：
per_peptide_counts %>% 
 ggplot() + aes( sample, count) + geom_jitter() + geom_violin() + scale_y_log10()
  1. 创建一个肽段命中数的累计图，并按命中数排序：
per_peptide_counts %>%
 arrange(count) %>%
 mutate(cumulative_hits = cumsum(count), peptide = 1:length(count)) %>%
 ggplot() + aes(peptide, cumulative_hits) + geom_line()
  1. 过滤掉非常低和非常高的肽段命中数，然后重新绘制它们：
filtered_per_peptide_counts <- per_peptide_counts %>%
  filter(count >= 5, count <= 2500) 

filtered_per_peptide_counts %>% 
 ggplot() + aes( sample, count) + geom_jitter() + geom_violin() + scale_y_log10()

它是如何工作的...

在步骤 1中，我们进行了一些库加载，并添加了数据加载步骤。如前所述，使用MSnID时，情况有些不同。我们不是直接调用文件读取函数，而是首先创建并清空MSnID对象，然后将数据加载到其中。我们通过MSnID()函数创建msnid，然后将其传递给read_mzid()函数，以实际将数据放入其中。接下来，我们使用as()函数将msnid转换为data.table对象——一个像数据框的对象，专门为大数据集优化。

在步骤 2中，我们使用tidyverse包中的dplyr和ggplot准备图形。tidyverse包相互协作得非常好，因为它们专注于处理数据框。通常的工作方式是使用管道操作符%>%将数据从一个函数传递到另一个函数，而无需保存中间对象。按照惯例，上游函数的结果作为下游函数的第一个参数传递，因此我们不需要显式指定它。这就形成了我们所见的结构。我们将peptide_info对象通过%>%操作符传递给dplyr filter()函数，后者执行工作并将结果传递给group_by()函数，依此类推。每个函数执行完后将数据传递给下一个函数。因此，在这个管道中，我们使用filter()保留所有不是诱饵的行，然后使用group_by(pepSeq)将长data.table根据pepSeq行的值分组为子表——实际上就是按肽序列获取一个表格。接下来的步骤使用summarise()，它生成一个包含count列的汇总表，count列的内容是n()函数的结果，n()函数统计表格中的行数，得出的表格每行代表一个肽，告诉我们这个肽在表中出现的次数。如果不清楚这些对象是如何逐步构建起来的，逐个函数地调试代码是个好主意。最后，我们使用mutate()添加一个名为sample的新列到表中，这列的长度与当前表相同，并填充为peptide_counts，然后将其添加到表中。该表被保存在名为per_peptide_counts的变量中。

在步骤 3中，我们将per_peptide_counts数据传递给ggplot()函数，该函数会设置一个ggplot对象。这些是内置层，因此我们使用+操作符通过aes()函数添加一个美学层。这个层通常包含要绘制在 x 轴和 y 轴上的变量——在这里，它们是sample和count。然后，我们再次使用+来添加一个geom——一个定义图形外观的层。首先，我们添加geom_jitter()，它绘制数据点，并在 x 轴和 y 轴上添加一些随机噪声，使点略微分散开。接着，我们添加另一个 geom，geom_violin()，它生成一个小提琴密度图。最后，我们添加一个尺度层，将尺度转换为以 10 为底的对数尺度。最终的图形如下所示：

在第 4 步中，我们通过将per_peptide_counts数据传递给arrange()函数，创建了一个累计命中图表。该函数按照指定的变量（在此案例中为计数）对数据框进行升序排序。结果被传递给mutate()，以添加一个新列cumulative_hits，它的值是cumsum()函数对计数列的计算结果。我们还添加了一个名为peptide的列，它获取表格的行号，同时也提供了一个方便的变量，使我们可以在图表中按顺序排列肽段。我们可以通过将排序后的数据直接传递给ggplot()并添加aes()函数来生成图表，使得peptide显示在 x 轴，cumulative_hits显示在 y 轴。然后，通过添加geom_line()，生成的图表如下所示：

从两个图表中，我们可以看到数据点的分布，并评估我们希望应用的阈值。

在第 5 步中，我们再次使用filter()函数，保留计数值大于 5 且小于 2500 的行，并将这些新数据放入我们在第 3 步中创建的同一绘图配方中。这将生成如下图表，显示了阈值外的点被移除：

转换 MS 格式以在工具之间移动数据

在生物信息学中，我们花费大量时间在不同文件格式之间转换，这是不可避免的事实。在这个简短的配方中，我们将介绍一些 R 语言中的方便方法，帮助我们在 MS 数据格式之间进行转换。

准备工作

对于这个配方，我们需要使用mzR包和来自本书仓库datasets/ch6文件夹中的threonine_i2_e35_pH_tree.mzXML文件。某些依赖项依赖于封装的 Java 代码，因此你需要为你的系统安装Java 运行时环境（JRE）；有关安装说明，请参考docs.oracle.com/goldengate/1212/gg-winux/GDRAD/java.htm。请在安装 R 包之前安装 JRE。

如何操作...

使用以下步骤可以在工具之间转换 MS 格式：

  1. 加载库并导入源数据文件：
library(mzR)
mzxml_file <- file.path(getwd(), "datasets", "ch6", "threonine_i2_e35_pH_tree.mzXML")
mzdata <- openMSfile(mzxml_file)
  1. 提取标题和峰数据：
header_info <- header(mzdata)
peak_data_list <- spectra(mzdata)
  1. 将数据写入新格式的文件：
writeMSData(peak_data_list, 
 file.path(getwd(), "datasets", "ch6", "out.mz"), 
 header = header_info, 
 outformat = "mzml", 
 rtime_seconds = TRUE 
)

它是如何工作的...

第一步是一个简单的数据加载步骤，我们在之前的配方中见过。我们使用openMSfile()函数，它会自动检测输入文件类型。

第 2 步是关键步骤；为了创建输出，我们需要制作一个标题对象和一个峰值列表。因此，我们使用header()和spectra()访问器函数从我们的mzdata对象中提取它们。输出函数将需要一个列表，因此如果文件中只有一个光谱，使用list()函数将spectra()函数包装起来。

最后的步骤是写入文件；在这里，第一个参数是峰值列表，第二个参数是要创建的文件名称，第三个参数是你选择的输出格式—可以选择mzml、mzxml和mzdata。最后一个参数表示保留时间是否以秒为单位编码；选择FALSE会将输出设置为以分钟为单位。

使用 protViz 进行光谱与肽段的匹配验证

尽管大多数光谱/肽段匹配是在高通量搜索引擎中完成的，但有时你可能希望检查竞争的模糊匹配的质量，或者与一个完全任意的感兴趣序列进行比较。运行整个搜索引擎管道可能是过于复杂的操作，因此在本教程中，我们将介绍一种便捷的方法，用于将单个光谱与单个肽段序列匹配，并生成理论离子大小与光谱中存在的离子之间的吻合图。

准备工作

对于本教程，我们只需要protViz包、mzR包以及本书仓库中datasets/ch6文件夹下的MM8.mzml文件。

如何实现...

使用protViz将光谱与肽段匹配可以通过以下步骤完成：

  1. 加载库和 MS 数据：
library(mzR)
library(protViz)
mzml_file <- file.path(getwd(), "datasets", "ch6", "MM8.mzML")
ms <- openMSfile(mzml_file)
  1. 从光谱中提取峰值和保留时间：
p <- peaks(ms,2)
spec <- list(mZ = p[,1], intensity = p[,2])
  1. 创建理论与观测离子质量的图表：
m <- psm("PEPTIDESEQ", spec, plot=TRUE)

它是如何工作的...

在步骤 1中，我们加载所需的库，并使用mzR函数openMSFile()来创建代表质谱数据的对象。

在步骤 2中，我们使用peaks()函数，该函数会将保留时间和峰值强度提取为矩阵对象。注意，第一列包含保留时间，而第二列包含强度。peaks()的第二个参数是我们需要的光谱的索引，因此我们正在获取该文件中的第二个光谱。如果省略此参数，则会返回所有光谱的列表。接下来的步骤中，我们需要将保留时间和强度数据封装在一个列表中，我们使用list()函数完成此操作，列表的成员命名为mZ和intensity。

最后，我们可以使用psm()函数生成图表。该函数的第一个参数是一个序列（这里是一个无意义的序列，保证匹配不佳），第二个参数是我们之前创建的光谱数据列表。通过将图表参数设置为TRUE，我们可以得到如下结果图表：

在图表中，每个点表示预测离子质量与光谱中最接近的观测质量之间的差异。我们可以看到，离子 b8、b7 和 c1 的质量偏差约为 1 Da，或者与任何预测质量相比更为分散，表明该肽段序列与光谱的拟合效果较差。

对光谱应用质量控制过滤器

原始蛋白质组学数据的质量控制是确保管道和分析能够给出可信和有用结果的关键步骤。需要大量的指标和数据图表来评估特定实验是否成功，这意味着在我们开始从数据中提取任何新知识之前，必须进行大量的分析。在这个配方中，我们将查看一个集成的管道，该管道执行一系列相关且有用的质量控制步骤，并将结果呈现为一个单一的、易于阅读的报告。

准备工作

在这个配方中，我们将研究大肠杆菌细胞膜的蛋白质组学实验。由于该文件太大，无法在本书的代码库中托管，因此我们将通过代码直接下载它。因此，在进行此实验时，您需要保持在线。我们还需要目标有机体的肽段文件，即Escherichia_coli.pep.all.fa文件，可以在本书代码库的datasets/ch6文件夹中找到。我们的主要功能将来自proteoQC库。

如何操作...

可以使用以下步骤将质量控制过滤器应用于光谱：

  1. 加载库并下载源数据：
library(proteoQC) 
online_file <- "ftp://ftp.pride.ebi.ac.uk/pride/data/archive/2017/11/PXD006247/CS_130530_ORBI_EMCP2156_b2469_narQ_DDM_AmH_X_5.mzXML"

mzxml_file <- file.path(getwd(), "datasets", "ch6", "PXD006247_mz.xml.gz" )

download.file(online_file, mzxml_file, "internal")
  1. 创建设计文件：
design_df <- data.frame(
 file = c(mzxml_file),
 sample = c(1),
 bioRep = c(1),
 techRep = c(1),
 fraction = c(1)
 )
design_file <- file.path(getwd(), "datasets", "ch6", "design_file.txt")
write.table(design_df, file = design_file, quote = FALSE, row.names = FALSE)
  1. 设置质量控制管道并运行以下命令：
qc <- msQCpipe(
 spectralist = design_file,
 fasta = file.path(getwd(), "datasets", "ch6", "Escherichia_coli.pep.all.fa"),
 outdir = file.path(getwd(), "qc_result"),
 enzyme = 1, varmod = 2, fixmod =1,
 tol = 10, itol = 0.6, cpu = 2,
 mode = "identification"
)

它是如何工作的...

在步骤 1中加载库后，我们设置要从www.proteomexchange.org/获取的文件的 URL；我们只需要获取PXD006247这个文件，并将 URL 保存在online_file变量中。我们还创建了一个指向本地文件系统中不存在的文件PXD006247_mz.xml.gzX的mzmxl_file变量——这将是下载文件保存的名称。download.file()函数实际上执行下载；第一个参数是在线源，第二个参数是文件下载到本地机器的存放位置。最后一个参数internal是下载方法。我们选择的设置应该使用一个系统无关的下载器，可以在任何地方使用，但如果你喜欢，也可以将其更改为其他更快或更符合系统的设置。文档将解释这些选项。

在步骤 2中，我们创建一个设计文件，描述实验。在我们的小示例中，我们只有一个文件，但你可以在此处指定更多文件。在第一部分，我们创建一个数据框，包含文件、样本、生物重复、技术重复和分级列。我们只有一个文件，所以表格只有一行。它看起来是这样的：

文件样本生物重复技术重复分级
PXD006247_mz.xml.gz1111

如果你的实验更为复杂，那么每个文件会有更多的行来描述样本和生物重复（bioRep）。接下来，我们使用write.table()和适当的选项将这个文件保存到磁盘，以便在下一步中使用。需要注意的是，虽然为了演示目的，我们是通过编程方式创建了这个文件，但如果我们通过电子表格程序或文本编辑器手动创建它，这个文件同样有效。

最后，我们在第 3 步中设置并运行 QC 流程。主要函数msQCpipe()是核心部分，需要一些选项设置。spectralist选项需要指向我们创建的设计文件路径，以便知道打开哪些文件以及如何处理它们。fasta选项需要目标生物体蛋白质序列的fasta格式文件。这使得 QC 流程能够使用rtandem包中的XTandem进行谱肽识别。outdir参数指定了一个新文件夹的路径，用于保存将要创建的众多报告文件。在这里，我们的文件夹将被命名为qc_result，并且它将是当前工作目录的一个子目录。enzyme、varmod和fixmod参数分别描述了用于消化的酶（1 = 胰蛋白酶）、可能存在的可变修饰以及所有残基上将存在的固定修饰。tol和itol参数指定了肽质量值的容差和误差窗口。cpu参数指定了源机器上要使用的计算核心数，而mode参数指定了运行的类型。

当 QC 流程完成后，我们会在qc_result文件夹中得到一系列报告。qc_report.html文件包含了可以浏览的 QC 结果。多个描述结果的页面应该能够让你了解实验的成功程度。

还有更多…

要找到合适的enzyme、varmod和fixmod变量值，你可以使用showMods()和showEnzymes()函数查看列表及其关键数字。

确定与肽匹配的基因组位点

找到基因组中肽匹配的确切位置可能是一个具有挑战性的任务，尤其是当基因组没有出现在原始搜索文件中时。在这个方法中，我们将结合使用经典的命令行 BLAST 方法，在翻译后的基因组序列中寻找短的、几乎精确的肽匹配，并通过针对 BLAST 命中的GRanges对象，将其与各种 R 基因组学管道相结合。

准备就绪

对于这个食谱，我们将使用MSnID、dplyr、withR、GenomicRanges和Biostrings包，以及来自大肠杆菌谱图的搜索引擎输出文件，该文件可以在本书的datasets/ch6文件夹中的PXD006247.mzXML.mzid文件中找到。你还需要本地安装 BLAST+版本。可以通过 conda 包管理器使用conda install -c bioconda blast命令安装。你还需要知道 BLAST+中的 tblastn 程序安装的位置。可以在 macOS 和 Linux 系统上使用终端命令which tblastn找到，Windows 系统亦然。

如何做...

可以使用以下步骤来识别与肽段匹配的基因组位点：

  1. 加载库和数据：
library(MSnID)
library(dplyr)
library(Biostrings)

msnid <- MSnID() # create object
msnid <- read_mzIDs(msnid, file.path(getwd(), "datasets", "ch6", "PXD006247.mzXML.mzid")) 

peptide_info <- as(msnid, "data.table") %>%
 filter(isDecoy == FALSE) %>%
 select(spectrumID, pepSeq, ) %>%
 mutate(fasta_id = paste0( spectrumID, ":", 1:length(spectrumID)) )
  1. 提取肽段序列并将其保存为 fasta 文件：
string_set <- AAStringSet(peptide_info$pepSeq )
names(string_set) <- peptide_info$fasta_id
writeXStringSet(string_set[1], file.path(getwd(), "datasets", "ch6", "peptides.fa"))
  1. 准备 BLAST 运行的文件名：
input_seqs <- file.path(getwd(), "datasets", "ch6", "peptides.fa")
genome_seqs <- file.path(getwd(), "datasets", "ch6", "ecoli_genome.fasta")
output_blast <- file.path(getwd(), "datasets", "ch6", "out.blast")
  1. 准备BLAST命令：
command <- paste0(
 "tblastn", 
 " -query ", input_seqs ,
 " -subject ", genome_seqs, 
 " -out ", output_blast,
 " -word_size 2 -evalue 20000 -seg no -matrix PAM30 -comp_based_stats F -outfmt 6 -max_hsps 1"
 )
  1. 将 BLAST 作为后台进程运行：
library(withr)
with_path("/Users/macleand/miniconda2/bin", system(command, wait = TRUE) )
  1. 将 BLAST 转换为GFF和GRanges：
results <- read.table(output_blast)

blast_to_gff <- function(blst_res){
 blst_res %>% 
 mutate(
 seqid = V2,
 source = rep("tblastn", length(V1)),
 type = rep(".", length(V1)),
 start = V9,
 end = V10,
 score = V3,
 strand = rep(".", length(V1)),
 phase = rep(".", length(V1)),
 attributes = paste("Name=",V1)
 ) %>%
 select( - starts_with("V") )
}

gff_df <- blast_to_gff(results)

library(GenomicRanges)
granges<-makeGRangesFromDataFrame(gff_df)

它是如何工作的...

步骤 1加载库，并使用MSnID包将数据加载到一个对象中，然后使用dplyr管道进行处理，正如本章中食谱 3的步骤 2所描述的。如果你不熟悉这种语法，可以查看那里进行深入解释。简而言之，虽然管道移除了假设行，但它仅保留spectrumID和pepSeq列，并添加了一个名为fasta_id的新列，该列将谱图 ID 粘贴为唯一编号。结果数据框被保存到peptide_info变量中。

步骤 2使用peptide_info$pepSeq列和peptide_info$fasta_id列的名称，通过names()函数创建一个Biostrings对象。然后，使用writeXStringSet()函数将结果BioStrings字符串集对象以 fasta 格式写入名为peptides.fa的文件。请注意，string_set末尾的索引[1]；这是一个小技巧，确保只写入第一个肽段。我们之所以这样做，仅仅是因为这是一个演示，且我们希望代码在短时间内完成。对于真实分析，你可以完全去掉索引，写入所有的序列。

在步骤 3中，我们只是设置了 BLAST 运行的输入和输出文件名。请注意，我们映射到的参考基因组ecoli_genome.fasta将位于本书的datasets/ch6文件夹中的仓库内。

在步骤 4中，我们指定了BLAST命令，这里的代码只是简单地将变量和文本粘贴在一起，形成一个长字符串并保存到命令中。值得详细查看的是，第一行指定了要运行的 BLAST+程序；这里是tblastn，它使用蛋白质输入和翻译后的核苷酸数据库。接下来的三行指定了输入的肽序列、用于 BLAST 的参考基因组，以及保存结果的输出文件。最后的长行指定了 BLAST+选项，以允许进行短小且几乎精确的匹配。设置了这些特定选项后，BLAST 运行可能需要一些时间，因此在开发过程中建议只运行一个序列。

在步骤 5中，指定了BLAST命令后，我们可以执行实际的 BLAST。我们这里的主要函数是基本的 R 函数system()，它将在后台运行系统命令。然而，为了帮助这个函数在不同系统之间移植，我们使用了withR库中的with_path()函数，它暂时将一个特定文件夹添加到系统的 PATH 中——这个 PATH 是包含程序的文件夹列表。这个步骤是必要的，因为有时 R 和 RStudio 不会识别非标准的安装位置，比如 conda 包管理器使用的那些位置。因此，这里第一个参数是tblastn文件夹的路径。注意，/Users/macleand/miniconda2/bin是我机器上的路径；你需要使用类似which tblastn的命令在终端或命令行中获取你机器上的路径，并进行替换。添加路径后，with_path()将运行其第二个参数，我们的system()函数，进而运行 BLAST。实际运行 BLAST 程序会花费一些时间。

一旦命令完成，在步骤 6中，我们首先通过read.table()函数将 BLAST 生成的输出文件加载到结果变量中。然后，我们创建一个自定义函数，将结果的行转换为 GFF 兼容的表格。blast_to_gff()函数使用dplyr mutate()函数添加相关列，然后使用select()函数与-选项选择不以字母 V 开头的列，因为所有原始列的名称都是以 V 开头的。现在我们可以使用GenomicRanges函数makeGRangesFromDataFrame()，将我们的 GFF 风格数据框转换为GRanges对象。这是最后一步，我们现在拥有一个匹配肽的基因组位置对象，可以在 R 的所有标准基因组学管道中使用，并且可以在本书中的基因组学配方中使用。

第七章：生成出版物和 Web-ready 可视化图形

设计和制作出版质量的可视化图形是生物信息学家的核心任务之一，也是他们在数据处理过程中最具成就感的部分。R 提供了许多优秀的图形制作包，超越了强大的基础图形系统和ggplot2。在本章的配方中，我们将探讨如何为许多不同类型的数据制作图表，这些数据通常不是典型的条形图/散点图类型。我们还将讨论网络、交互式图形、3D 图形和圆形基因组图。

本章将介绍以下配方：

  • 使用 ridgeplots 可视化多个分布

  • 为双变量数据创建色图

  • 将关系数据表示为网络

  • 使用 plotly 创建交互式 Web 图形

  • 使用 plotly 构建三维图形

  • 构建多组学数据的圆形基因组图

技术要求

本章所需的示例数据可以在本书的 GitHub 仓库中找到：github.com/danmaclean/R_Bioinformatics_Cookbook. 如果你想按照书中的代码示例运行，你需要确保数据位于工作目录的子目录中。

以下是你需要的 R 包。你可以使用install.packages("package_name")来安装它们。列在Bioconductor下的包需要使用专用的安装器进行安装，安装步骤也在本节中描述。如果你需要做其他操作，安装步骤将在包使用的配方中描述：

  • circlize

  • dplyr

  • ggplot2

  • ggridges

  • gplots

  • plotly

  • RColorBrewer

  • readr

  • magrittr

  • tidyr

  • viridis

Bioconductor非常庞大，且有自己的安装管理器。你可以使用以下代码安装该管理器：

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

然后，你可以使用以下代码安装这些包：

BiocManager::install("package_name")

更多信息可以在www.bioconductor.org/install/中找到。

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

有时，在配方的中间，我会打断代码流程，让你看到一些中间输出或对理解非常重要的对象结构。每当发生这种情况，你会看到一个代码块，其中每行以##双哈希符号开头。请考虑以下命令：

letters[1:5]

这将产生以下输出：

## a b c d e

请注意，输出行以##为前缀。

使用 ridgeplots 可视化多个分布

可视化某些测量数量的分布是生物信息学中的一项常见任务，R 的基础功能通过hist()和density()函数以及通用的plot()方法，能很好地处理这种任务，后者可以创建对象的图形。ggplot图形系统以一种整洁的方式绘制多个密度图，每个因子水平对应一个图形，从而生成一个紧凑且非常易读的图形——即所谓的山脊图。在这个示例中，我们将看看如何创建一个山脊图。

准备工作

在这个示例中，我们将使用ggplot和ggridges包。数据集方面，我们将使用datasets包中的一个数据集，该包通常是 R 的预安装包之一。我们将使用airquality数据集。你可以直接在 R 控制台中输入airquality来查看它。

如何实现...

使用山脊图可视化多个分布可以通过以下步骤完成：

  1. 加载库：
library(ggplot2)
library(ggridges)
  1. 构建一个ggplot描述：
ggplot(airquality) + aes(Temp, Month, group = Month) + geom_density_ridges()
  1. 显式将Month转换为因子：
ggplot(airquality) + aes(Temp, as.factor(Month) ) + geom_density_ridges()
  1. 给山脊上色：
ggplot(airquality) + aes(Temp, Month, group = Month,  fill = ..x..) + 
  geom_density_ridges_gradient() + 
  scale_fill_viridis(option = "C", name = "Temp")
  1. 重塑数据框并添加面板：
library(tidyr)
airquality %>%
  gather(key = "Measurement", value = "value", Ozone, Solar.R, Wind, Temp ) %>%
  ggplot() + aes(value, Month, group = Month) + 
  geom_density_ridges_gradient() +
  facet_wrap( ~ Measurement, scales = "free")

工作原理...

在步骤 1中加载所需的库后，步骤 2中我们使用ggridges包的geom_ridges()函数创建了一个标准的ggplot描述。如果你以前没见过ggplot图，它们非常简单。一个ggplot图有三个层次，至少需要通过三个函数构建——ggplot()函数始终是第一个，它允许我们指定数据集。接下来，使用+运算符添加的函数是aes()函数或美学函数，我们可以将其视为我们希望在图表中看到的内容。第一个参数表示X轴上的元素，而第二个参数表示Y轴上的元素。group = Month参数是特定于山脊图的，它告诉绘图函数如何对数据点进行分组。由于Month数据是数字型，而不是因子型，因此这里需要此参数。最后，我们添加geom_density_ridges()来创建正确类型的图。

在步骤 3中，我们遵循了与步骤 2类似的程序，但这次，我们作为替代方法使用了as.factor(Month)，它在处理和渲染分组之前显式地将Month数据转换为因子类型。因此，Month的步骤变得不再必要。这些步骤生成的图形如下，步骤 2在左，步骤 3在右：

在第 4 步中，我们为山脊添加了颜色。基本上，ggplot 构建过程没有变化，只是在 aes() 函数中增加了 fill = ..x..，它告诉绘图填充 x 轴方向的颜色。然后我们使用了一个稍微不同的 geom 函数，geom_density_ridges_gradient()，它能够为山脊着色。在最后一层，我们使用 scale_fill_viridis()，选择了来自 viridis 颜色库的颜色尺度（该库在顶部加载）。"C" 选项指定了特定的颜色尺度，而 name 用来指定尺度的名称。结果图像如下所示：

最后，在第 5 步中，我们通过进一步的维度对数据进行了拆分，并在相同样式的图中添加了包含同一数据集不同方面的面板。为了实现这一点，airquality 数据需要进行一些预处理。我们加载了 tidyr 包，并使用 gather() 函数将命名列的值（具体是 Ozone、Solar.R、Wind 和 Temp）合并到一个名为 value 的列中，并添加一个名为 Measurement 的新列，记录原始列中观察值的来源。然后，我们将结果传递到 ggplot()。构建过程几乎与之前相同（请注意，我们的 x 轴现在是 value，而不是 Temp，因为温度数据已经存储在重塑后的数据框中），最后添加了 facet_wrap() 函数，它使用公式符号选择数据的子集并显示在各个面板中。选项 scales 为 "free"，允许每个结果面板有独立的坐标轴刻度。结果如下所示：

为双变量数据创建颜色映射

颜色映射，也称为热图，是二维矩阵的绘图，其中数值在特定尺度上转换为颜色。我们可以通过多种方式在 R 中绘制这些图表；大多数图形包都提供了这类功能。在本教程中，我们将使用基础包的 heatmap() 函数来可视化一些矩阵。

准备工作

我们只需要 ggplot 包和内置的 WorldPhones 数据集。

操作步骤...

为双变量数据创建颜色映射可以通过以下步骤完成：

  1. 创建基本的热图：
heatmap(WorldPhones)
  1. 移除树状图：
heatmap(WorldPhones, Rowv = NA, Colv = NA)
  1. 为分组添加颜色尺度：
cc <- rainbow(ncol(WorldPhones), start = 0, end = .3)

heatmap(WorldPhones, ColSideColors = cc)
  1. 更改调色板：
library(RColorBrewer)
heatmap(WorldPhones, ColSideColors = cc, 
        col = colorRampPalette(brewer.pal(8, "PiYG"))(25))
  1. 重新缩放数据：
heatmap(WorldPhones, ColSideColors = cc, scale = "column")

它是如何工作的...

在第 1 步中，我们将基础的 heatmap() 函数应用于一个矩阵，返回一个如下所示的图：

在第 2 步中，我们使用了 Rowv 和 Colv 参数来去除树状图。请注意，在结果图中，列的顺序与矩阵中的顺序一致。通过使用树状图，我们可以重新排列列和行。没有树状图的图表如下所示：

在第 3 步中，我们使用rainbow()函数创建了一个调色板对象，该函数返回绘图所需的颜色。rainbow()的第一个参数是颜色的数量。这里，我们使用ncol(WorldPhones)来为数据集的每一列获取一种颜色。start和end参数指定了彩虹中颜色选择的起始和结束位置。然后，我们可以在ColSideColors参数中使用CC调色板对象，以获取列的颜色键。我们可以使用更多相似的列来获得更多相似的颜色，具体如下：

在第 4 步中，我们为col参数提供了一个调色板对象，以更改热图的整体调色板。我们使用了colorRampPalette()函数，从一小部分特定颜色列表中生成顺序调色板。这个函数会插值生成完整的调色板。我们将RColorBrewer包中的brewer.pal()函数作为参数传给colorRampPalette()，该函数根据提供的选项返回粉色-黄色-绿色（PiYG）调色板中的八种颜色。生成的热图颜色如下所示：

最后，在第 5 步中，我们在可视化步骤中对数据进行了数值转换。我们使用heatmap()的scale选项来规范化图表中的数据。请注意，设置为column时按列进行规范化，而设置为row时按行进行规范化。默认的基础包scale()函数用于此操作。图表中的数字重新缩放是颜色变化的原因，而不是直接从调色板中选择的结果。图表如下所示：

另见

heatmap()函数随后有了其他包，这些包遵循类似的语法，但扩展了其功能。尝试使用gplots包中的heatmap.2()和heatmap.3()。以下直方图展示了heatmap.2()的绘图效果。它与heatmap()非常相似，但默认添加了颜色键和直方图：

表示关系数据为网络

网络，或图形，是表示实体之间关系的强大数据表示方式，对于许多生物学研究至关重要。网络分析可以揭示生态学研究中的社区结构，揭示蛋白质-蛋白质相互作用中的潜在药物靶点，并帮助我们理解复杂代谢反应中的相互作用。表示网络的基础数据结构可能是复杂的。幸运的是，R 有一些非常强大的包，特别是igraph和ggraph，我们可以使用它们来访问网络信息并进行绘图。在本教程中，我们将探索一些生成适当大小网络图的方式。

准备工作

在这个示例中，我们需要ggraph和igraph包及其依赖项，包括magrittr、readr和dplyr。我们还需要来自本书代码库中datasets/ch7文件夹的bio-DM-LC.edges文件。这个文件包含来自 WormNet 的一些基因功能关联数据。网络包含大约 1,100 条边和 650 个节点。你可以在这里阅读更多关于数据的信息：networkrepository.com/bio-DM-LC.php。

如何操作...

将关系数据表示为网络可以通过以下步骤完成：

  1. 加载包并准备数据框：
library(ggraph)
library(magrittr)

df <- readr::read_delim(file.path(getwd(), "datasets", "ch7", "bio-DM-LC.edges"), 
                        delim = " ", 
                        col_names = c("nodeA", "nodeB", "weight")) %>% 
      dplyr::mutate(edge_type = c("A","B"), length(weight), replace = TRUE))
  1. 创建一个igraph对象并在基础图形中使用它：
graph <- igraph::graph_from_data_frame(df)

ggraph(graph, layout = "kk") + 
  geom_edge_link() + 
  geom_node_point() + 
  theme_void()
  1. 根据边的值或类型给边着色：
ggraph(graph, layout = "fr") + 
  geom_edge_link(aes(colour = edge_type)) + 
  geom_node_point() + 
  theme_void()
  1. 添加节点属性并相应地为节点着色：
igraph::V(graph)$category <- sample(c("Nucleus", "Mitochondrion", "Cytoplasm"), length(igraph::V(graph)), replace = TRUE )
igraph::V(graph)$degree <- igraph::degree(graph)

ggraph(graph, layout = "fr") + 
  geom_edge_link(aes(colour = edge_type)) + 
  geom_node_point(aes(size = degree, colour = category)) + 
  theme_void()

它是如何工作的...

在第 1 步中，我们加载了所需的库，并从边列表文件中准备了数据框。输入文件基本上是一个边列表。每一行描述了一个连接，第一列是目标节点之一，第二列是另一个目标节点。第三列包含一个值，表示这两个节点之间交互的强度，我们可以将其视为边的权重。字段由空格分隔，文件没有列名标题。因此，我们适当地设置了delim和col_names参数的值。然后，我们将数据框传递给dplyr::mutate()函数，添加一个名为edge_type的额外列。在该列中，我们使用sample()函数随机为每一行分配"A"或"B"。最终的对象保存在df变量中。

在第 2 步中，我们使用igraph::graph_from_data_frame()函数从df创建了igraph对象，并将其保存在graph变量中。我们将igraph图形对象作为第一个对象传递给ggraph()函数，后者与ggplot()类似。它接收graph对象和一个布局参数。（在这里，我们使用"kk"，但具体使用哪个布局将高度依赖于数据本身。）然后，我们使用+运算符添加层。首先，我们添加了geom_edge_link()层，它绘制了边，然后是geom_node_point()，它绘制了节点，最后，我们添加了theme_void()，它去除了背景的灰色面板和白色线条，为网络留下了一个清晰的背景。初始图形如下所示：

在第 3 步中，我们添加了一些基于数据的自定义设置。我们首先将布局算法更改为"fr"，这种算法提供了更美观且更分散的视图。然后，我们在geom_edge_link()中使用aes()函数将边的颜色映射到edge_type值。剩下的层像之前一样添加。通过这样做，我们得到了以下图形：

在 第 4 步 中，我们为节点设置了一些属性。这比看起来要简单。igraph 中的 V() 函数返回一个简单的节点 ID 向量（在 igraph 行话中，节点被称为顶点），因此我们计算该向量的长度，并用它生成一个包含 Nucleus、Mitochondrion 和 Cytoplasm 值的随机向量。然后我们可以通过使用带有 $ 索引的 V() 函数将这些新值分配给节点。我们可以创建任何我们喜欢的属性，因此 igraph::V(graph)$category 创建了一个名为 category 的新属性。我们可以使用标准的 *<-* 赋值运算符将新值直接分配给该属性。接下来的步骤类似；igraph::V(graph)$degree 创建了一个名为 degree 的属性。在我们的案例中，我们将 igraph::degree() 函数的结果分配给它。Degree 是图论中的术语，表示与一个节点相连的边的数量。现在我们拥有了新的属性，可以根据这些属性对图形进行上色。ggraph() 构建过程与之前相同，但在 geom_node_point() 层中，我们使用 aes() 将颜色映射到我们新的 category 属性，并将大小映射到我们新的 degree 属性。最终的图形如下所示：

还有更多内容...

蜂窝图是绘制网络的好方法，尤其是在你有三种节点类型或某种方向性结构时。你可以使用我们已经拥有的相同类型的数据创建一个蜂窝图，如下所示：

ggraph(graph, 'hive', axis = 'category') +
  geom_edge_hive(aes(colour = edge_type, alpha = ..index..)) + 
  geom_axis_hive(aes(colour = category)) + 
  theme_void()

在这里，我们设置布局类型为 hive，并指定用于制作轴的属性为 category。在 geom_edge_hive() 中的边描述与之前差不多，alpha 参数为 ..index..，它根据边的绘制顺序为边缘添加透明度元素。geom 节点被替换为 geom_axis_hive()，我们在其中使用 aes() 将颜色映射到 category 属性。生成的图形如下所示：

使用 plotly 创建交互式网页图形

通过图形用户界面交互式地探索数据集是一种有益且启发性的方式，能够分析和审视数据。动态地在图形中添加和删除数据、放大或缩小特定部分，或让图形随着时间变化（依赖于底层数据）可以让我们看到一些在静态图形中看不见的趋势和特征。在这个实例中，我们将使用 plotly 库在 R 中创建交互式图形，从一个基本图形开始，逐步构建一个更复杂的图形。

准备工作

在这个实例中，我们将使用内置的 Orange 数据集，它描述了橙树树干围长随时间的变化。该数据集是 datasets 包的一部分（通常已经预装），因此你应该能够立即访问它。

如何操作...

使用 plotly 创建交互式网页图形可以通过以下步骤完成：

  1. 加载库并创建一个基本图形：
library(plotly)
plot_ly(data = Orange, x = ~age, y = ~circumference)
  1. 映射标记的颜色和大小以及悬停文本到数据：
plot_ly(data = Orange, x = ~age, y = ~ circumference,
        color = ~Tree, size = ~age,
        text = ~paste("Tree ID: ", Tree, "<br>Age: ", age, "Circ: ", circumference)
        )
  1. 添加第二个系列/迹线：
trace_1 <- rnorm(35, mean = 120, sd = 10)
new_data <- data.frame(Orange, trace_1)

plot_ly(data = new_data, x = ~age, y = ~ circumference,
        color = ~Tree, size = ~age,
        text = ~paste("Tree ID: ", Tree, "<br>Age: ", age, "Circ: ", circumference)

)  %>% add_trace(y = ~trace_1, mode = 'lines' ) %>%
  add_trace(y = ~circumference, mode = 'markers' )
  1. 添加一个下拉菜单，以便您可以选择图表类型：
plot_ly(data = Orange, x = ~age, y = ~ circumference,
        color = ~Tree, size = ~age,
        text = ~paste("Tree ID: ", Tree, "<br>Age: ", age, "Circ: ", circumference)

)  %>%
  add_trace(y = ~circumference, mode = 'marker' ) %>%
  layout(
    title = "Plot with switchable trace",
      updatemenus = list(
        list(
          type = "dropdown",
            y = 0.8,
            buttons = list(
              list(method = "restyle",
                   args = list("mode", "markers"),
                   label = "Marker"
              ),
              list(method = "restyle",
                   args = list("mode", "lines"),
                   label = "Lines"
                   )

            )
          )
      )
  )

它是如何工作的...

在第 1 步加载库后，我们使用核心的plot_ly()函数创建了可能最简单的图。我们向plot_ly()传递了数据框的名称，以及作为公式的x和y轴的列，因此使用了~符号。在这一点上，我们还没有明确指定迹线类型，即plotly称其为系列或数据轨道，所以它猜测并创建了一个散点图，如下图所示：

注意图表顶部的菜单图标以及当鼠标悬停在数据点上时显示的悬停文本。这些图形在交互式 R 会话中可以很好地进行交互，但更适合于 HTML-based 文档，例如编译的 R markdown。

在第 2 步中，我们将图中的特征映射到数据的各个方面。我们将大小和颜色设置为与树木 ID 和年龄列相对应的公式，再次使用~语法。我们还为每个数据点设置了悬停文本，并使用paste()函数编译了确切的格式。请注意，悬停文本是基于 HTML 的，我们可以使用诸如<br>这样的标签来按照选择的方式格式化悬停效果。我们的图现在改进为如下所示：

在第 3 步中，我们的主要变化是明确指定了迹线数据。为了突出显示迹线可以携带超出原始数据框之外的数据，我们使用rnorm()创建了一个名为trace_1的新数据向量，其中包含 35 个均值为 120，标准差为 1 的随机数。我们与在第 2 步创建图表的方式相同，但这次我们使用了magrittr管道将图形对象发送到add_trace()函数。在这里，我们将新的trace_1对象作为我们的y值传递，并将mode设置为"lines"以获得折线图。再次，我们将其传递到另一个add_trace()函数（我们可以通过多个迹线系列构建一个图），但这次使用原始数据框列周长，并将mode设置为"markers"。生成的图表看起来像这样：

在第 4 步中，我们将菜单引入了我们的图表。我们实现的菜单将允许我们在不同的追踪类型之间切换——从线条到标记，再到线条。此步骤以相同的基本plot_ly()调用开始，然后这次只传入一个追踪。接下来，我们传入了layout函数，该函数接受一个图表标题作为title参数，并为updatemenus参数提供一个复杂的选项列表。你必须将一个包含三个成员的列表传递给updatemenus——type、y和buttons。type设置菜单类型——在这种情况下，我们想要一个下拉菜单；y设置菜单在y轴上的位置，值范围为 0 到 1，而buttons需要另一个列表，其中每个子列表描述一个菜单选项。每个子列表都有members方法，以及args和labels。setting方法用于"restyle"，意味着图表将在选择菜单时更新。args成员需要另一个列表，指定当前菜单选项的"mode"和"type"。最后，label指定此菜单选项在菜单中显示的文本。当我们在下拉菜单中选择Marker时，图表如下所示，渲染在左侧：

使用 plotly 构建三维图表

我们在生物信息学中生成的大多数图表都是静态的，并且局限于二维平面，但现代网络技术使我们能够与三维对象进行互动并动态渲染。plotly库提供了渲染不同类型 3D 图表的工具，在这个教程中，我们将展示如何构建一个 3D 表面图和一个带有x、y、z轴的散点图。

准备就绪

在这个教程中，我们将再次使用plotly库，并且使用内置的longley经济数据集。

如何做...

使用plotly构建三维图表可以通过以下步骤完成：

  1. 设置数据对象：
library(plotly)

d <- data.frame(
  x <- seq(1,10, by = 0.5),
  y <- seq(1,10, by = 0.5)
)

z <- matrix(rnorm(length(d$x) * length(d$y) ), nrow = length(d$x), ncol = length(d$y))
  1. 创建基本的表面图：
plot_ly(d, x = ~x , y = ~y, z = ~z) %>% 
  add_surface()
  1. 添加一个响应式等高线图层：
plot_ly(d, x = ~x , y = ~y, z = ~z) %>% 
  add_surface(
    contours = list(
    z = list(
      show=TRUE,
      usecolormap=TRUE,
      highlightcolor="#ff0000",
      project=list(z=TRUE)
      )
    )
  )

它是如何工作的...

在第 1 步中，我们首先构建一个适合图表类型的数据集。对于表面三维图，我们需要一个包含x和y坐标的数据框，我们直接使用data.frame()函数创建，并将其保存在一个名为d的变量中。数据框d包含了 x 和 y 列之间从 0 到 10 的 20 个值（每列 10 个值）。你可以将这个数据框视为指定三维场的宽度和长度，而不是实际的数据值。数据值是通过与数据框指定的维度相对应的矩阵对象获得的。我们使用matrix()函数创建了一个适当维度的矩阵，矩阵中的值是通过rnorm()函数生成的随机正态值。一旦我们有了这两个结构，就可以像之前为二维图表做的那样，使用它们在plot_ly()函数中指定d、x和y，并添加新的z轴来获取我们的矩阵。结果被传递给add_surface()函数，将数据渲染为三维表面。该图表的效果类似如下：

请注意，点击并拖动图表区域，你可以调整相机视角。

在第 2 步中，我们通过在三维表面下方（或上方）添加一个反应式等高线图来进一步完善图表。我们在add_surface()函数中使用了 contours 选项。它接受一个选项列表。第一个z指定如何处理等高线。它还接受一个包含多个成员的列表，用于控制等高线图的外观，最重要的选项是highlightcolor，它指定了在等高线图上绘制的颜色，用于显示鼠标悬停的当前三维图层。渲染后的图像类似于以下效果：

在第 3 步中，我们改变了方向，绘制了一个三维散点图。这更加直接。我们将 Longley 数据传递给plot_ly()函数，并指定了维度和数据列进行映射。我们还添加了一个标记选项，将颜色映射到 GNP 列。最后，我们将基本图对象传递给add_markers()函数，得到最终图表，效果如下所示：

构建多组学数据的圆形基因组图

对多条数据序列的整体基因组分析通常以圆形方式呈现，采用同心圆，每个圆显示不同种类的数据，每个数据在圆中的表示方式不同。这些图表被称为 Circos 图，具有极大的表现力，可以在紧凑的形式中显示大量的密集信息。在本教程中，我们将学习如何在 R 中构建这样的图表，使用常见的基因组数据文件。

准备工作

为了制作 Circos 图，我们将使用circlize包以及本书仓库中datasets/ch7/文件夹下的四个以arabidopsis为前缀的文件。

如何操作...

构建多组学数据的圆形基因组图可以通过以下步骤完成：

  1. 加载库并读取染色体长度信息：
library(circlize)
df <- readr::read_tsv(file.path( getwd(), "datasets", "ch7", "arabidopsis.gff"), col_names = FALSE) %>% 
  dplyr::select(X1, X4, X5)
  1. 初始化图表和染色体轨迹，然后添加链接：
circos.genomicInitialize(df)
circos.link("Chr4", c(9000000, 1200000), 
            "Chr5", c(12000000,15000000),
            col = "red")
  1. 从文件加载链接信息并绘制：
circos.clear()

source_links <- read.delim(file.path(getwd(), "datasets", "ch7", "arabidopsis_out_links.bed"), header = FALSE) 
target_links <- read.delim(file.path(getwd(), "datasets", "ch7", "arabidopsis_in_links.bed"), header = FALSE)

circos.genomicInitialize(df)
circos.genomicLink(source_links, target_links, col = "blue")
  1. 加载基因位置并添加密度轨迹：
circos.clear()

gene_positions <- read.delim(file.path(getwd(), "datasets", "ch7", "arabidopsis_genes.bed"), header = FALSE)
circos.genomicInitialize(df)
circos.genomicDensity(gene_positions, window.size = 1e6, col = "#0000FF80", track.height = 0.1)
  1. 加载热图数据。然后，添加一个热图轨迹：
circos.clear()

heatmap_data <- read.delim(file.path(getwd(), "datasets", "ch7", "arabidopsis_quant_data.bed"), header = FALSE)

col_fun = colorRamp2(c(10, 12, 15), c("green", "black", "red"))
circos.genomicInitialize(df)
circos.genomicHeatmap(heatmap_data, col = col_fun, side = "inside", border = "white")
  1. 合并这些轨迹：
circos.clear()

circos.genomicInitialize(df)
circos.genomicHeatmap(heatmap_data, col = col_fun, side = "inside", border = "white")
circos.genomicDensity(gene_positions, window.size = 1e6, col = "#0000FF80", track.height = 0.1)
circos.genomicLink(source_links, target_links, col = "blue")

它是如何工作的...

在步骤 1中，我们开始通过读取arabidopsis.gff文件来构建图表，该文件描述了我们希望在图表中使用的染色体长度。我们只需要名称、开始和结束列，因此我们将数据传递给dplyr::select()函数，保留适当的列，即X1、X4和X5。由于.gff文件没有列标题，read_tsv()函数会给出列名 X1 ... Xn。我们将结果保存在df对象中。

在步骤 2中，我们开始构建图表。我们使用circos.genomicInitialize()函数和df来创建图表的骨架和坐标系统，然后手动添加了一个单独的链接。circos.link()函数允许我们使用染色体名称和 c(start, end)格式创建单个来源和目标，从而将链接着色为请求的颜色。当前图表如下所示：

在步骤 3的开始，我们使用了circos.clear()来完全重置图表。重置仅在本教程中需要，因为我们想逐步构建内容；在你自己的编码中，你可能可以忽略它。下一步是加载一个基因组区域文件，代表一些链接的来源，以及一个独立的基因组区域文件，代表一些链接的目标。这两个文件应该是 BED 格式，并且源文件中的第 N 行必须与目标文件中的第 N 行对应。然后，我们使用circos.genomicInitialize()重新初始化图表，并使用circos.genomicLink()在一个命令中添加多个链接，传入源链接数据和目标数据的对象，然后将它们全部着色为蓝色。图表如下所示：

在步骤 4中，在清除图表后，我们从arabidopsis_genes.bed文件中读取另一个基因位置的 BED 文件。我们希望将这些信息添加为一个密度轨迹，计算用户指定长度的窗口中功能的数量，并将其绘制为密度曲线。为此，我们使用circos.genomicDensity()函数，传入gene_positions数据框，选择一个 1 百万的窗口大小，一种颜色（注意颜色使用八位十六进制格式，允许我们为颜色添加透明度），以及track.height，它指定了图表中此轨迹所占的比例。轨迹如下所示：

在第 5 步中，我们添加了一个更复杂的轨迹——一个热图，可以表示许多列的定量数据。这里的文件格式是扩展的 BED 格式，包含染色体名称、起始位置和结束位置，以及后续列中的数据。在我们的示例文件arabidopsis_quant_data.bed中有三个额外的数据列。我们通过read.delim()将 BED 文件加载到heatmap_data中。接下来，我们创建了一个颜色函数，并将其保存为col_fun，用于绘制热图。colorRamp2()函数接受数据的最小值、中值和最大值的向量作为其参数，然后使用第二个参数中指定的颜色。因此，使用10、12和15以及绿色、红色和黑色，我们分别将 10 绘制为绿色，12 绘制为黑色，15 绘制为红色。介于这些值之间的颜色会由colorRamp2()自动计算。为了绘制热图，我们使用了circos.genomicHeatmap()函数，并将col_fun传递给col参数。side参数指定是否在圆内或圆外绘制，而border参数指定热图元素之间线条的颜色。绘图效果如下：

最后，在第 6 步中，我们将所有内容结合在一起。通过清除并重新初始化绘图，我们通过按从外到内的顺序调用相关函数来指定轨迹的顺序：

最终的图形，如前面所见，先是circos.genomicHeatmap()，然后是circos.genomicDensity()，最后是circos.genomicLink()，以得到圆形基因组图。

第八章：与数据库和远程数据源的工作

大规模模式生物有机体测序项目，如人类基因组计划（HGP）或 1001 植物基因组测序项目，已经使大量基因组数据公开可用。同样，个别实验室的开放数据共享也使基因组和转录组的原始测序数据得到了广泛的共享。通过编程处理这些数据意味着可能需要解析或本地存储一些非常庞大或复杂的文件。因此，许多努力已投入使这些资源通过 API 和其他可查询接口（如 BioMart）尽可能易于访问。在本章中，我们将介绍一些食谱，帮助我们在不下载整个基因组文件的情况下搜索注释，并能够在数据库中找到相关信息。我们还将了解如何在代码中提取实验的原始读取数据，并借此机会研究如何对下载的数据进行质量控制。

本章将覆盖以下食谱：

  • 从 BioMart 检索基因和基因组注释

  • 检索和处理 SNP

  • 获取基因本体信息

  • 从 SRA/ENA 查找实验和读取数据

  • 对高通量测序读取进行质量控制和过滤

  • 完成使用外部程序的读到参考比对

  • 可视化读到参考比对的质量控制图

技术要求

你所需的示例数据可以在本书的 GitHub 仓库中找到，地址为 github.com/PacktPublishing/R-Bioinformatics-Cookbook. 如果你想按原样使用代码示例，你需要确保这些数据位于你的工作目录的子目录中。

以下是你需要的 R 包。一般来说，你可以使用install.packages("package_name")来安装这些包。在Bioconductor下列出的包需要使用专用的安装器进行安装。具体安装方法将在本节中描述。如果需要进一步操作，安装过程将在包使用的食谱中进行描述：

  • Bioconductor

    • biomaRt

    • ramwas

    • ShortRead

    • SRAdb

Bioconductor 非常庞大，并且拥有自己的安装管理器。你可以使用以下代码来安装管理器：

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

然后，你可以使用以下代码安装这些包：

BiocManager::install("package_name")

更多信息请参见 www.bioconductor.org/install/。

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

有时，在配方中间，我会中断代码，这样您就可以看到一些中间输出或者重要对象的结构。每当这种情况发生时，您将看到每行以##（双重井号）符号开头的代码块。考虑以下命令：

letters[1:5]

这将给我们以下输出：

## a b c d e

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

从BioMart检索基因和基因组注释

一旦准备好了某个基因组序列的草稿，就会进行大量的生物信息学工作，以找到基因和其他功能特征或重要的基因座。这些注释很多，执行和验证都很困难，通常需要大量的专业知识和时间，并且不希望重复。因此，基因组项目联合体通常会通过某种方式共享他们的注释，通常是通过一些公共数据库。BioMart是一种常见的数据结构和 API，通过它可以获取注释数据。在这个配方中，我们将看看如何以编程方式访问这些数据库，以获取我们感兴趣的基因的注释信息。

准备工作

对于这个配方，我们需要名为biomaRt的Bioconductor包以及一个可用的互联网连接。我们还需要知道要连接的BioMart服务器 —— 全球约有 40 个这样的服务器，提供各种信息。最常访问的是Ensembl数据库，这些是这些包的默认设置。您可以在这里查看所有BioMart的列表：www.biomart.org/notice.html。我们将开发的代码将适用于这些BioMart中的任何一个，只需稍微修改表名和 URL。

如何实现...

可以通过以下步骤从BioMart检索基因和基因组注释：

  1. 列出所选示例数据库gramene中的mart列表：
library(biomaRt)
listMarts(host = "ensembl.gramene.org")
  1. 创建到所选mart的连接：
gramene_connection <- useMart(biomart = "ENSEMBL_MART_PLANT", host = "ensembl.gramene.org")
  1. 列出该mart中的数据集：
data_sets <-  listDatasets(gramene_connection)
head(data_sets)

data_set_connection <- useMart("atrichopoda_eg_gene", biomart = "ENSEMBL_MART_PLANT", host = "ensembl.gramene.org")
  1. 列出我们实际可以检索的数据类型：
attributes <- listAttributes(data_set_connection)
head(attributes)
  1. 获取所有染色体名称的向量：
chrom_names <- getBM(attributes = c("chromosome_name"), mart = data_set_connection )
head(chrom_names)
  1. 创建一些用于查询数据的过滤器：
filters <- listFilters(data_set_connection)
head(filters)
  1. 获取第一个染色体上的基因 ID：
first_chr <- chrom_names$chromosome_name[1]
genes <- getBM(attributes = c("ensembl_gene_id", "description"), filters = c("chromosome_name"), values = c(first_chr), mart = data_set_connection )head(genes)
head(genes)

工作原理...

该配方围绕对数据库进行一系列不同的查找操作，每次获取一些更多信息来处理。

在第一步中，我们使用listMarts()函数获取指定主机 URL 上所有可用的BioMart列表。在需要连接到不同服务器时，请根据需要更改 URL。我们得到一个可用mart的数据框，并使用该信息。

在第二步中，我们使用useMart()函数创建一个名为gramene_connection的连接对象，传入服务器 URL 和第一步中的具体BioMart。

在步骤 3中，我们将gramene_connection传递给listDatasets()函数，以检索该biomart中的数据集。选择其中一个数据集（atrichopda_eg_gene）后，我们可以运行useMart()函数来创建一个到该biomart中数据集的连接，并将对象命名为data_set_connection。

在步骤 4中，我们几乎完成了确定可以使用哪些数据集的工作。在这里，我们使用在listAttributes()函数中创建的data_set_connection，获取我们可以从该数据集中检索的各种信息类型的列表。

在步骤 5中，我们最终通过主要函数getBM()获取一些实际信息。我们将attributes参数设置为我们希望返回的数据的名称；在这里，我们获取chromosome_name的所有值，并将它们保存到一个向量chrom_names中。

在步骤 6中，我们设置过滤器——即接收哪些值的限制。我们首先询问data_set_connection对象，使用listFilters()函数查看我们可以使用哪些过滤器。从返回的filters对象中可以看到，我们可以在chromosome_name上进行过滤，所以我们将使用这个。

在步骤 7中，我们设置一个完整的查询。在这里，我们打算获取第一个染色体上的所有基因。请注意，我们已经在步骤 5中获得了一个染色体列表，因此我们将使用chrom_names对象中的第一个元素作为过滤器，并将其保存在first_chr中。为了执行查询，我们使用getBM()函数，并指定ensembl_gene_id和description属性。我们将filter参数设置为我们希望过滤的数据类型，并将values参数设置为我们希望保留的过滤器值。我们还将data_set_connection对象作为要使用的 BioMart 传递。最终生成的genes对象包含了第一个染色体上的ensembl_gene_id和描述信息，如下所示：

## ensembl_gene_id           description
## 1 AMTR_s00001p00009420    hypothetical protein 
## 2 AMTR_s00001p00015790    hypothetical protein 
## 3 AMTR_s00001p00016330    hypothetical protein 
## 4 AMTR_s00001p00017690    hypothetical protein 
## 5 AMTR_s00001p00018090    hypothetical protein 
## 6 AMTR_s00001p00019800    hypothetical protein

获取和处理 SNP

SNP 和其他多态性是重要的基因组特征，我们通常希望在特定基因组区域内检索已知的 SNP。在这里，我们将介绍如何在两个不同的 BioMart 中执行此操作，这些 BioMart 存储了不同类型的 SNP 数据。在第一部分中，我们将再次使用 Gramene 来查看如何获取植物 SNP。在第二部分中，我们将了解如何在主要的 Ensembl 数据库中查找人类 SNP 的信息。

准备工作

如前所述，我们只需要biomaRt包，它来自Bioconductor，并且需要一个正常工作的互联网连接。

如何操作……

获取和处理 SNP 可以通过以下步骤完成：

  1. 从 Gramene 获取数据集、属性和过滤器列表：
library(biomaRt)
listMarts(host = "ensembl.gramene.org")
gramene_connection <- useMart(biomart = "ENSEMBL_MART_PLANT_SNP", host = "ensembl.gramene.org")
data_sets <- listDatasets(gramene_connection)
head(data_sets)
data_set_connection <- useMart("athaliana_eg_snp", biomart = "ENSEMBL_MART_PLANT_SNP", host = "ensembl.gramene.org")

listAttributes(data_set_connection)
listFilters(data_set_connection)
  1. 查询实际的 SNP 信息：
snps <- getBM(attributes = c("refsnp_id", "chr_name", "chrom_start", "chrom_end"), filters = c("chromosomal_region"), values = c("1:200:200000:1"), mart = data_set_connection )
head(snps)

它是如何工作的……

步骤 1将与之前的食谱中的步骤 1到6类似，我们将建立初始连接并让它列出我们可以在此 BioMart 中使用的数据集、属性和过滤器；这是相同的模式，每次使用 BioMart 时都会重复（直到我们能熟记它）。

在步骤 2中，我们利用收集到的信息提取目标区域的 SNP。同样，我们使用getBM()函数并设置chromosomal_region过滤器。这允许我们指定一个描述基因组中特定位点的值。value参数采用Chromosome:Start:Stop:Strand格式的字符串；具体而言，1:200:20000:1，这将返回染色体 1 上从第 200 到 20000 个碱基的所有 SNP，且位于正链上（注意，正链 DNA 的标识符为1，负链 DNA 的标识符为-1）。

还有更多...

从 Ensembl 中查找人类 SNP 的步骤基本相同。唯一的区别是，由于 Ensembl 是默认服务器，我们可以在useMart()函数中省略服务器信息。对于人类的类似查询会像这样：

data_set_connection <- useMart("hsapiens_snp", biomart = "ENSEMBL_MART_SNP")
human_snps <- getBM(attributes = c("refsnp_id", "allele", "minor_allele", "minor_allele_freq"), filters = c("chromosomal_region"), value = c("1:200:20000:1"), mart = data_set_connection)

另见

如果你拥有dbSNP refsnp ID编号，可以通过rnsps包和ncbi_snp_query()函数直接查询这些 ID。只需将有效的refsnp ID 向量传递给该函数。

获取基因本体论信息

基因本体论（GO）是一个非常有用的受限词汇，包含用于基因和基因产物的注释术语，描述了注释实体的生物过程、分子功能或细胞成分。因此，这些术语在基因集富集分析和其他功能-组学方法中非常有用。在本节中，我们将看看如何准备一个基因 ID 列表，并为它们获取 GO ID 和描述信息。

准备工作

由于我们仍在使用biomaRt包，所以只需要该包以及一个有效的互联网连接。

如何操作...

获取基因本体论信息的步骤如下：

  1. 连接到 Ensembl BioMart 并找到适当的属性和过滤器：
library(biomaRt)

ensembl_connection <- useMart(biomart = "ENSEMBL_MART_ENSEMBL")
 listDatasets(ensembl_connection)

data_set_connection <- useMart("hsapiens_gene_ensembl", biomart = "ENSEMBL_MART_ENSEMBL")

att <- listAttributes(data_set_connection)
fil <- listFilters(data_set_connection)
  1. 获取基因列表，并使用它们的 ID 获取 GO 注释：
genes <- getBM(attributes = c("ensembl_gene_id"), filters = c("chromosomal_region"), value = c("1:200:2000000:1"), mart = data_set_connection)

go_ids <- getBM(attributes = c("go_id", "goslim_goa_description"), filters = c("ensembl_gene_id"), values = genes$ensembl_gene_id, mart = data_set_connection )

工作原理...

如同前两个步骤，步骤 1 包括找到适合的 biomart、数据集、属性和过滤器的值。

在步骤 2中，我们使用getBM()函数获取特定染色体区域中的ensembl_gene_id属性，将结果保存在genes对象中。然后我们再次使用该函数，使用ensembl_gene_id作为过滤器，go_id和goslim_goa_description来获取仅选定基因的 GO 注释。

查找 SRA/ENA 中的实验和读取数据

短序列数据归档（SRA）和欧洲核苷酸库（ENA）是记录原始高通量 DNA 序列数据的数据库。每个数据库都是相同高通量序列数据集的镜像版本，这些数据集是来自世界各地各个生物学领域的科学家提交的。通过这些数据库免费获取高通量序列数据意味着我们可以设想并执行对现有数据集的新分析。通过对数据库进行搜索，我们可以识别出我们可能想要处理的序列数据。在本配方中，我们将研究如何使用SRAdb包查询 SRA/ENA 上的数据集，并以编程方式检索选定数据集的数据。

准备工作

这个配方的两个关键元素是来自Bioconductor的SRAdb包和一个有效的互联网连接。

如何操作...

查找来自 SRA/ENA 的实验和数据读取可以通过以下步骤完成：

  1. 下载 SQL 数据库并建立连接：
library(SRAdb)
sqlfile <- file.path(system.file('extdata', package='SRAdb'), 'SRAmetadb_demo.sqlite')
sra_con <- dbConnect(SQLite(),sqlfile)
  1. 获取研究信息：
dbGetQuery(sra_con, "select study_accession, study_description from study where study_description like '%coli%' ")
  1. 获取关于该研究包含内容的信息：
sraConvert( c('ERP000350'), sra_con = sra_con )
  1. 获取可用文件的列表：
listSRAfile( c("ERR019652","ERR019653"), sra_con, fileType = 'sra' )
  1. 下载序列文件：
getSRAfile( c("ERR019652","ERR019653"), sra_con, fileType = 'fastq', destDir = file.path(getwd(), "datasets", "ch8") )

它是如何工作的...

加载库后，第一步是设置一个本地的 SQL 文件，名为sqlfile。该文件包含了关于 SRA 上研究的所有信息。在我们的示例中，我们使用的是包内的一个小版本（因此，我们通过system.file()函数提取它）；真实的文件大小超过 50GB，因此我们暂时不会使用它，但可以通过以下替换代码获取：sqlfile <- getSRAdbfile()。一旦我们拥有一个sqlfile对象，就可以使用dbConnect()函数创建与数据库的连接。我们将连接保存在名为sra_con的对象中，以便重用。

接下来，我们使用dbGetQuery()函数对sqlfile数据库执行查询。该函数的第一个参数是数据库文件，第二个参数是一个完整的 SQL 格式的查询。写出的查询非常直观；我们希望返回包含coli一词的描述时的study_accession和study_description。更复杂的查询是可能的——如果你准备好用 SQL 编写它们。关于这方面的教程超出了本配方的范围，但有许多书籍专门介绍这个主题；你可以尝试阅读 Upom Malik、Matt Goldwasser 和 Benjamin Johnston 合著的《SQL for Data Analytics》，由 Packt Publishing 出版：www.packtpub.com/big-data-and-business-intelligence/sql-data-analysis。该查询返回一个类似于以下内容的 dataframe 对象：

## study_accession    study_description
## ERP000350    Transcriptome sequencing of E.coli K12 in LB media in early exponential phase and transition to stationary phase

步骤 3 使用我们提取的访问号，利用sraConvert()函数获取与该研究相关的所有提交、样本、实验和运行信息。这将返回如下表格，我们可以看到该研究的运行 ID，展示了包含序列的实际文件：

##    study submission    sample experiment       run
## 1 ERP000350 ERA014184 ERS016116 ERX007970 ERR019652 
## 2 ERP000350 ERA014184 ERS016115 ERX007969 ERR019653

在步骤 4中，我们使用 listSRAfile() 函数获取运行中特定序列的实际 FTP 地址。这将提供 SRA 格式文件的地址，这是一个压缩且方便的格式，如果你希望了解的话：

     run     study    sample experiment    ftp
## 1 ERR019652 ERP000350 ERS016116 ERX007970 ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/ERR/ERR019/ERR019652/ERR019652.sra 
## 2 ERR019653 ERP000350 ERS016115 ERX007969 ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/ERR/ERR019/ERR019653/ERR019653.sra

但是在步骤 5中，我们使用 getSRAfile() 函数，并将 fileType 参数设置为 fastq，以便以标准的 fastq 格式获取数据。文件将下载到 destDir 参数指定的文件夹中。

还有更多...

不要忘记定期刷新本地的 SQL 数据库，并使用以下代码获取完整版本：sqlfile <- getSRAdbfile()。

对高通量序列读取进行质量控制和过滤

当我们有一组新的序列读取需要处理时，无论是来自新实验还是数据库，我们都需要执行质量控制步骤，去除任何序列接头、去除质量差的读取，或者根据需要修剪质量差的序列。在这个配方中，我们将介绍如何在 R 中使用 Bioconductor ShortRead 包来实现这一点。

准备工作

你需要安装 ShortRead 包，并且需要运行本章中 查找 SRA/ENA 中的实验和读取 这一配方的代码。在该配方的最后一步，会创建两个文件，我们将使用其中一个。运行该代码后，文件应该位于本书的 datasets/ch8/ERRR019652.fastq.gz 目录中。

如何执行...

对高通量序列读取进行质量控制和过滤可以通过以下步骤完成：

  1. 加载库并连接到文件：
library(ShortRead)
fastq_file <- readFastq(file.path(getwd(), "datasets", "ch8", "ERR019652.fastq.gz") )
  1. 过滤掉质量低于 20 的任何核苷酸的读取：
qualities <- rowSums(as(quality(fastq_file), "matrix") <= 20) 
fastq_file <- fastq_file[qualities == 0]
  1. 修剪读取的右侧：
cut_off_txt <- rawToChar(as.raw(40))
trimmed <- trimTails(fastq_file, k =2, a= cut_off_txt)
  1. 设置自定义过滤器以移除 N 和同源重复序列：
custom_filter_1 <- nFilter(threshold=0)
custom_filter_2 <- polynFilter(threshold = 10, nuc = c("A", "T", "C", "G"))
custom_filter <- compose(custom_filter_1, custom_filter_2)
passing_rows <- custom_filter(trimmed)
trimmed <- trimmed[passing_rows]
  1. 写出保留的读取：
writeFastq(trimmed, file = file.path(getwd(), "datasets", "ch8", "ERR019652.trimmed.fastq.gzip"), compress = TRUE)

它是如何工作的...

第一步将读取加载到 ShortReadQ 对象中，该对象表示 DNA 读取及其相关的质量评分；这个特殊的对象使我们能够在一次操作中处理序列和质量。

第二步让我们找到所有质量分数都高于 20 的读取。这里的代码有点特殊，所以请花时间理解它。首先，我们在 fastq_file 上使用 quality() 函数提取质量值，然后将其传递给 as() 函数，要求返回一个矩阵。接着，我们对该矩阵使用 rowSums() 计算每行的总和，最终通过比较得到一个逻辑向量 qualities，它指示哪些 rowSums() 的值小于 20。在下一行中，我们使用 qualities 向量来对 fastq_file 进行子集筛选，去除低质量的读取。

在步骤 3中，我们修剪读取的右侧（以纠正读取质量低于阈值的地方）。这里的主要功能是trimTails()，它接受两个参数：k，开始修剪所需的失败字母数，以及a，开始修剪的字母。这当然意味着我们所认为的 Phred 数值质量分数（如步骤 2中，我们仅使用了 20）需要根据质量分数的文本编码转换为其 ASCII 等价物。这就是第一行发生的事情；数字 40 通过as.raw()转换为原始字节，然后通过rawToChar()转换为字符。生成的文本可以通过将其存储在cut_off_txt变量中来使用。

步骤 4应用了一些自定义过滤器。第一行，custom_filter_1，为包含名为N的碱基的序列创建过滤器，阈值参数允许序列包含零个N。第二行，custom_filter_2，为长度等于或大于阈值的同质聚合物读取创建过滤器。nuc参数指定要考虑的核苷酸。一旦指定了过滤器，我们必须使用compose()函数将它们合并成一个单一的过滤器，该函数返回一个我们称之为custom_filter()的过滤器函数，然后对修剪后的对象进行调用。它返回一个SRFFilterResult对象，可以用来对读取进行子集化。

最后，在步骤 5中，我们使用writeFastQ()函数将保留的读取写入文件。

使用外部程序完成读取到参考的比对

高通量读取的比对是本书中许多方法的一个重要前提，包括 RNAseq 和 SNP/INDEL 调用。在第一章，执行定量 RNAseq，以及第二章，使用 HTS 数据查找遗传变异中我们对其进行了深入讨论，但我们没有涉及如何实际执行比对。我们通常不会在 R 中执行此操作；进行这些比对所需的程序是强大的，并且作为独立进程从命令行运行。但 R 可以控制这些外部进程，因此我们将探讨如何运行外部进程，以便你可以从 R 包装脚本中控制它们，最终使你能够开发端到端的分析管道。

准备中...

本教程仅使用基础的 R 语言，因此你无需安装任何额外的包。你需要准备参考基因组 FASTA 文件datasets/ch8/ecoli_genome.fa，以及我们在寻找 SRA/ENA 实验和读取数据教程中创建的datasets/ch8/ERR019653.fastq,gz文件。本教程还需要系统中安装 BWA 和samtools的可执行文件。相关软件的网页可以在samtools.sourceforge.net/和bio-bwa.sourceforge.net/找到。如果你已经安装了conda，可以通过conda install -c bioconda bwa和conda install -c bioconda samtools来安装它们。

如何操作……

使用以下步骤完成读取到参考的比对，借助外部程序：

  1. 设置文件和可执行文件路径：
bwa <- "/Users/macleand/miniconda2/bin/bwa"
samtools <- "/Users/macleand/miniconda2/bin/samtools"
reference <- file.path(getwd(), "datasets", "ch8", "ecoli_genome.fa")
  1. 准备index命令并执行：
command <- paste(bwa, "index", reference)
system(command, wait = TRUE)
  1. 准备alignment命令并执行：
reads <- file.path(getwd(), "datasets", "ch8", "ERR019653.fastq.gz")
output <- file.path(getwd(), "datasets", "ch8", "aln.bam")
command <- paste(bwa, "mem", reference, reads, "|", samtools, "view -S -b >", output)
system(command, wait = TRUE)

它是如何工作的……

第一步非常简单：我们仅仅创建了几个变量来保存程序和文件所在目录的路径。bwa和samtools分别保存了这些程序在系统中的路径。请注意，你的系统上的路径可能与此不同。在 Linux 和 macOS 系统中，你可以通过which命令在终端中查找路径；在 Windows 机器上，你可以使用where命令或等效命令。

在步骤 2中，我们概述了运行系统命令的基本模式。首先，通过paste()函数，我们将命令作为字符串创建，并保存在一个名为command的变量中。在这里，我们正在准备一个命令行，在执行 BWA 进行读取比对之前创建所需的索引。然后，我们将该命令作为第一个参数传递给system()函数，后者会实际执行该命令。命令作为一个全新的进程在后台启动，默认情况下，一旦进程开始，控制权会立即返回给 R 脚本。如果你打算在后台进程输出后立即在 R 中继续工作，你需要将system()函数中的wait参数设置为TRUE，这样 R 进程将在后台进程完成后才继续。

在步骤 3中，我们扩展了模式，创建了读取和输出变量，并组合了一个更长的命令行，展示了如何构建一个有效的命令行。然后我们重复执行system命令。此过程将生成一个最终的 BAM 文件，位于datasets/ch8/aln.bam。

可视化读取到参考比对的质量控制

一旦完成读取的比对，通常建议检查比对的质量，确保读取的模式和期望的插入距离等没有异常。这在草拟参考基因组中尤其有用，因为高通量读取的异常比对可能会揭示参考基因组的拼接错误或其他结构重排。在本教程中，我们将使用一个名为ramwas的包，它有一些易于访问的图形，可以帮助我们评估比对的质量。

正在准备中...

对于这个配方，我们需要准备好的 bam_list.txt 和 sample_list.txt 信息文件，位于本书的 datasets/ch8 目录下。我们还需要来自同一位置的小文件 ERR019652.small.bam 和 ERR019653.small.bam。

如何执行...

可通过以下步骤可视化读取到参考基因组的比对质量控制：

  1. 设置运行的参数：
library(ramwas)
param <- ramwasParameters( dirbam = ".", filebamlist = "bam_list.txt", 
                            dirproject = file.path(getwd(), "datasets", "ch8"), 
                            filebam2sample = "sample_list.txt")
  1. 执行质量控制：
ramwas1scanBams(param)
qc <- readRDS(file.path(getwd(), "datasets", "ch8", "rds_qc", "ERR019652.small.qc.rds")$qc
  1. 查看图表：
plot(qc$hist.score1)
plot(qc$bf.hist.score1)
plot(qc$hist.length.matched)
plot(qc$bf.hist.length.matched)

它是如何工作的...

步骤 1 使用 ramwasParameters() 函数设置一个包含参数的对象。我们只需提供信息文件（bam_list.txt 和 sample_list.txt），分别指定要使用的 BAM 文件及其包含的样本位置。dirproject 参数指定结果应该写入系统的路径。请注意，这些结果会被写入磁盘，而不是直接返回内存。

步骤 2 使用参数通过 ramwas1scanBams() 函数运行质量控制（QC）。结果被写入磁盘，因此我们使用基本的 R readRDS() 函数加载返回的 RDS 文件。qc 对象包含许多成员，代表了不同的比对质量控制方面。

步骤 3 使用通用的 plot 函数来创建一些 qc 对象中的质量控制统计图。

avatar
文章
阅读
粉丝