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

R 生物信息学秘籍（三）

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

译者：飞龙

协议：CC BY-NC-SA 4.0

第九章：有用的统计和机器学习方法

在生物信息学中，对不同大小和组成的数据集进行统计分析是常见的任务。R 无疑是一个功能强大的统计语言，拥有丰富的选项来处理各种任务。在本章中，我们将重点关注一些有用但不常讨论的方法，尽管这些方法本身并不构成完整的分析，但它们可以成为你经常进行的分析的有力补充。我们将查看模拟数据集的食谱，以及用于类别预测和降维的机器学习方法。

本章将涵盖以下食谱：

• 校正 p 值以考虑多重假设

• 生成一个模拟数据集来表示背景

• 在数据中学习分组并使用 kNN 进行分类

• 使用随机森林预测类别

• 使用 SVM 预测类别

• 在没有先验信息的情况下对数据进行分组学习

• 使用随机森林识别数据中最重要的变量

• 使用 PCA 识别数据中最重要的变量

技术要求

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

这里是你需要的 R 包。通常，你可以使用install.packages("package_name")来安装这些包。列在Bioconductor下的包需要使用专门的安装器安装。如果需要进一步的操作，安装过程将在使用这些包的食谱中进行描述：

• Bioconductor

  • Biobase

• caret

• class

• dplyr

• e1071

• factoextra

• fakeR

• magrittR

• randomForest

• RColorBrewer

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

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

校正 p 值以考虑多重假设

在生物信息学，特别是在基因组学项目中，我们常常在分析中执行数千次统计检验。但这可能会导致我们结果中的显著误差。考虑一个基因表达实验，每个处理的测量数很少（通常只有三次），但基因数量有成千上万。一个执行统计检验的用户，在p <= 0.05的情况下，会错误地拒绝零假设 5%的时间。对进行多个假设检验的校正可以帮助我们减少此类分析中的错误率。我们将查看一种简单的校正方法。

准备开始

所有需要的函数都是 R 的基础函数，我们将通过代码创建自己的数据。

如何实现...

校正 p 值以考虑多个假设的步骤如下：

1. 运行 10,000 次 t 检验：

set.seed(1)
random_number_t_test <- function(n){
  x <- rnorm(10)
  y <- rnorm(10)
  return(t.test(x,y)$p.value)
}

p_values <- sapply(1:10000, random_number_t_test )

1. 评估 p*-*值数量，<= 0.05：

sum(p_values <= 0.05)

1. 调整 p 值：

adj_p_values <- p.adjust(p_values, method = "holm")

1. 重新评估<= 0.05的 p 值数量：

sum(adj_p_values <= 0.05)

它是如何工作的...

步骤 1中的第一行代码简单地固定了随机数生成器，以便我们能够在不同计算机之间获得一致的结果；除了为了对比本书中的结果，你不需要这部分代码。接下来是创建一个自定义函数，生成两组（x和y）10 个随机数，然后执行 t 检验并返回 p 值。由于这些只是来自相同分布的随机数，因此没有实际差异。最后一行使用 sapply()函数运行我们自定义的函数并创建一个包含 10,000 个 p 值的向量。

在步骤 2中，我们仅仅统计 p 值小于 0.05 的数量。我们得到如下结果：

## [1] 506

这表示我们有 506 个错误判定为显著的结果。

在步骤 3中，我们使用p.adjust()函数应用一种校正方法。argument方法可以是几种可用方法之一。实践中，最好尝试holm或BH（Benjamini Hochberg），因为这些方法提供了准确的假阳性率。一种广泛使用但效果不太好的方法是Bonferroni；大多数情况下应避免使用它。

在步骤 4中，我们重新评估小于 0.05 的 p 值数量。这次，结果正如我们预期的那样：

## [1] 0

生成一个代表背景的模拟数据集

构建模拟数据集以作为合理的对照，进行与预期背景分布的适当比较，并拥有一个合适的背景人群来抽取样本，这是许多研究的重要方面。在本配方中，我们将探讨从头开始或通过混合现有数据框来生成这些数据的各种方式。

准备开始

我们将使用fakeR包和内置的iris数据集。

如何实现...

生成一个代表背景的模拟数据集可以按照以下步骤进行：

1. 创建一个具有与给定数据集相同特征的随机数据集：

library(fakeR)
fake_iris <- simulate_dataset(iris)

1. 创建一个均值和标准差与给定向量相同的正态分布随机数向量：

sample_mean <- mean(iris$Sepal.Length)
sample_sd <- sd(iris$Sepal.Length)
random_sepal_lengths <- rnorm(iris$Sepal.Length, mean = sample_mean, sd = sample_sd)
hist( random_sepal_lengths)

1. 在一个范围内创建一个均匀分布的随机整数向量：

low_num <- 1
high_num <- 6
hist(runif(1500, low_num, high_num))

1. 创建一个二项分布成功次数的向量：

number_of_coins <- 1
p_heads <- 0.5
hist(rbinom(1500, number_of_coins, p_heads ))
number_of_coins <- 5
hist(rbinom(1500, number_of_coins, p_heads ))

1. 从列表中随机选择元素，且每个元素的选择概率不同：

random_from_list <- sample(c("Low", "Medium", "High"), 100, replace = TRUE, prob = c(0.2, 0.6, 0.2))
table(random_from_list)

它是如何工作的……

步骤 1 使用了fakeR包中的simulate_dataset()函数，生成一个新的数据集，该数据集与源数据集（iris）具有相同数量的值、相同的列名、相同数量的因子水平和水平名称，以及相同数量的行。值是随机化的，但数据框架是完全相同的。注意，使用str()函数报告iris和新的fake_iris对象的结构是完全相同的：

str(iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...

str(fake_iris)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.26 6.69 5.63 5.21 5.28 6.45 6.8 5.71 6.01 6.44 ...
##  $ Sepal.Width : num  2.84 2.78 2.83 2.44 2.19 3.87 3.14 2.58 2.78 3.25 ...
##  $ Petal.Length: num  4.03 4.84 2.64 2.83 5.37 3.63 5.54 4.74 4.63 4.29 ...
##  $ Petal.Width : num  1.63 1.33 0.7 0.61 2.03 1.17 2.05 1.6 1.57 1.32 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 3 2 2 3 1 2 1 3 3 1 ...

在步骤 2中，我们的目标是生成一个随机数向量，它的均值和标准差与iris数据集中Sepal.Length列的均值和标准差相同。为此，我们首先使用mean()和sd()函数计算这些数值。然后，我们将这些值作为参数传递给rnorm()函数的mean和sd参数。通过运行hist()绘制生成的random_sepal_lengths向量，我们可以确认其分布和参数。

在步骤 3中，我们希望创建一个数值（浮动点）向量，这些数值能够以相等的概率出现——这类似于反复掷骰子：每个选项的可能性是一样的。事实上，在这个步骤中，我们将范围的低值（low_num）设置为 1，将高值（high_num）设置为 6，以模拟这个过程。我们通过runif()函数生成 1,500 个符合这些低值和高值的随机数，并通过再次使用hist()绘制结果，我们可以看到每个区间内的频率相对均匀，从而确认这些数值的均匀性。

在步骤 4中，我们希望模拟一个掷硬币式的概率实验——所谓的二项成功概率分布。我们首先必须决定每次实验的试验次数——在掷硬币实验中，这指的是我们掷的硬币数量。在这里，我们将number_of_coins变量设置为 1。我们还必须决定成功的概率。同样，模拟掷硬币实验时，我们将p_heads变量设置为 0.5。为了进行模拟，我们将这些值传递给rbinom()函数，请求进行 1,500 次独立的实验重复。hist()函数显示，0 次成功（掷出反面）和 1 次成功（掷出正面）的频率在所有 1,500 次重复中大致相等。接下来，我们通过改变number_of_coins变量的值，将试验次数改为 5，模拟每次重复实验时使用五枚硬币。我们再次使用rbinom()函数并通过hist()绘制结果，这时我们可以观察到，两个和三个成功（正面）是每次五枚硬币实验中最常见的结果。

最后，在步骤 5中，我们通过sample()函数来选择向量中的项目。sample的第一个参数是要从中采样的向量——这里是整数 1 到 10。第二个参数是要选择的项目数量——这里我们选择 10 个。请注意，默认情况下，sample()会进行不放回抽样，因此每个项目不会被选择两次，尽管向量中的每个项目每次都有相等的被选中概率。sample()的第二种用法将replacement参数设置为TRUE，意味着可以重复选择项目。此用法还设置了prob参数——一个包含选择初始向量中每个值的概率的向量。运行此采样并将结果传递给table()函数确认我们获得的选择符合预期的近似概率。

数据中的学习分组与 kNN 分类

**k-最近邻（kNN）**算法是一种监督学习算法，它会根据给定的数据点，尝试根据其与一组已知类别的训练样本的相似性对其进行分类。在这个食谱中，我们将学习如何使用数据集，将其划分为测试集和训练集，并根据训练集构建的模型预测测试集的类别。这种方法广泛应用于生物信息学，并且在聚类分析中非常有价值，尤其是当我们有一些已知的目标类别示例时。

准备工作

对于这个食谱，我们需要几个新的包：caret、class、dplyr和magrittr。作为数据集，我们将使用内置的iris数据集。

如何进行...

使用以下步骤可以在数据中进行学习分组并进行 kNN 分类：

1. 标准化数据并删除非数值列：

set.seed(123)
scaled_iris <- iris %>% mutate_if( is.numeric, .funs = scale)
labels <- scaled_iris$Species
scaled_iris$Species <- NULL

1. 提取训练集和测试集：

train_rows <- sample(nrow(scaled_iris), 0.8 * nrow(scaled_iris), replace = FALSE)
train_set <- scaled_iris[train_rows, ]
test_set <- scaled_iris[-train_rows, ]
train_labels <- labels[train_rows]
test_set_labels <- labels[-train_rows]

1. 构建模型并对测试集进行预测：

test_set_predictions <- knn(train = train_set, test = test_set, cl = train_labels, k = 10)

1. 比较预测结果与实际类别：

caret::confusionMatrix(test_set_predictions,  test_set_labels)

它是如何工作的...

在步骤 1中，我们首先使用set.seed()确保随机数的可复现性，然后使用dplyr mutate_if()函数对数据集的每一列进行标准化。mutate_if()的第一个参数是要测试的条件；.funs参数是如果条件为真则应用的函数。在这里，我们对iris数据框的一列应用scale()函数，如果该列是数值型，则返回我们称之为scaled_iris的数据框。列间的标准化在 kNN 中非常重要，因为实际值的大小可能会有很大影响，所以我们需要确保它们在各列之间的尺度相似。接下来，我们从数据中复制Species列，因为它包含类标签，并通过将列赋值为NULL将其从数据框中删除——接下来的步骤中，数据框应仅包含数值数据。

在步骤 2中，我们决定哪些行应该包含在训练集和测试集中。我们使用sample()函数从 1 到iris行数的向量中选择；我们无替代地选择 80%的行号，因此train_rows是一个整数向量，表示scaled_iris中的行，我们将在训练集中使用它。在这一步的其余部分，我们使用子集和负子集来准备我们需要的scaled_iris子集。

在步骤 3中，我们应用 kNN 算法，通过knn()函数构建模型，并在一次操作中对测试集进行分类。train参数表示我们为训练预留的数据部分，test参数表示我们为测试预留的数据部分，cl（类别）参数表示训练集的标签。k参数是用于分类每个未知测试点的邻居数。该函数返回一个包含测试数据中每行预测类别的向量，我们将其保存在test_set_predictions中。

在步骤 4中，我们使用caret包的confusionMatrix()函数评估预测结果。该函数接受预测类别和真实类别，并生成一组统计数据，包括以下表格，其中行表示Real标签，列表示Predicted标签。该模型错误地将一个versicolor行预测为virginica，但其余所有预测都正确：

##             Reference
## Prediction   setosa versicolor virginica
##   setosa          8          0         0
##   versicolor      0          9         1
##   virginica       0          0        12

使用随机森林进行类别预测

随机森林是另一种监督学习算法，它使用决策树的集合来进行多次类别预测，最终预测最常出现的类别作为模型的最终预测。随机森林通常是有用的，因为它可以同时处理分类和数值数据，并且可以应用于分类和回归任务。在本章的使用随机森林识别数据中最重要的变量配方中，我们将再次使用它来预测数据中最重要的变量。在这个配方中，我们将使用随机森林来预测数据的类别。

做好准备

对于这个配方，我们需要caret和randomForest包以及内置的iris数据集。

如何执行此操作...

使用随机森林进行类别预测可以通过以下步骤完成：

1. 从iris数据集中准备训练集：

library(randomForest)

train_rows <- sample(nrow(iris), 0.8 * nrow(iris), replace = FALSE)
train_set <- iris[train_rows, ]
test_set <- iris[-train_rows, ]

1. 在训练数据上构建模型：

model <- randomForest(Species ~ . , data = train_set, mtry = 2)

1. 使用模型对测试数据进行预测：

test_set_predictions <- predict(model, test_set, type = "class")
caret::confusionMatrix(test_set_predictions,  test_set$Species)

它是如何工作的...

步骤 1的整个过程是准备训练集和测试集。我们使用sample()函数从 1 到iris行数的向量中选择；我们无替代地选择 80%的行号，因此train_rows是一个整数向量，表示iris中的行，我们将在训练集中使用它。在这一步的其余部分，我们使用子集和负子集来准备我们需要的iris子集。

在步骤 2 中，我们直接进行模型构建和预测。randomForest()函数以其第一个参数为 R 公式，命名要预测的列（即Species，响应变量），并且是训练数据的数据框列——在这里，我们使用所有列，表示为一个.字符。data参数是源数据框的名称，mtry参数是一个可调参数，告诉算法使用多少分裂。这个最佳值通常是列数的平方根，但优化它可能会有帮助。生成的模型保存在一个名为model的变量中，可以打印进行检查。

在步骤 3 中，我们使用predict()函数与model、test_set数据和设置为class的type参数来预测测试集的类别。然后，我们使用caret::confusionMatrix()来评估它们，得到以下结果：

##             Reference
## Prediction   setosa versicolor virginica
##   setosa         13          0         0
##   versicolor      0          8         0
##   virginica       0          0        9
## 

结果表明测试集完美分类。

还有更多内容

可以使用非常相似的方法进行回归（预测数值）预测。看看以下建立回归模型和进行评估的代码的相似性。在这里，我们基于其他列预测萼片长度。模型构建后，我们像以前一样运行预测；注意我们如何删除type参数（因为实际上回归是默认的）。最后，我们通过计算均方误差（MSE）来评估，在其中我们平方预测值和萼片长度的实际值之间的差异，然后取均值：

model <- randomForest(Sepal.Length ~ . , data = train_set, mtry = 2)
test_set_predictions <- predict(model, test_set)

mean( (test_set$Sepal.Length - test_set_predictions )² ) 

使用 SVM 预测类别

支持向量机（SVM）算法是一种分类器，通过在数据的多个维度中找到类之间的最大距离——有效地是类之间最大的间隙——并使用该间隙的中点作为分类的边界。在这个配方中，我们将使用 SVM 来执行监督类别预测，并通过图形化地说明边界。

准备就绪

我们将继续使用内置的iris数据集和e1071包。

如何做...

使用 SVM 预测类别可以通过以下步骤完成：

1. 构建训练集和测试集：

library(e1071)
train_rows <- sample(nrow(iris), 0.8 * nrow(iris), replace = FALSE)

train_set <- iris[train_rows, ]
test_set <- iris[-train_rows, ]

1. 构建模型：

model <- svm(Species~., data=train_set, type="C-classification", kernel="radial", gamma=0.25)

1. 绘制模型边界：

cols_to_hold <- c("Sepal.Length", "Sepal.Width")
held_constant <- lapply(cols_to_hold, function(x){mean(train_set[[x]])})
names(held_constant) <- cols_to_hold

plot(model, train_set, Petal.Width ~ Petal.Length, slice = held_constant)

1. 在测试集上进行预测：

test_set_predictions <- predict(model, test_set, type = "class")
caret::confusionMatrix(test_set_predictions,  test_set$Species)

工作原理...

在步骤 1 中，我们有可能熟悉的训练集和测试集生成步骤，我们在前面的配方中讨论过。简而言之，在这里，我们创建一个行号向量作为训练集，并使用子集和负子集来提取新的子数据集。

在第 2 步中，我们使用svm()函数创建模型。第一个参数是一个 R 公式，指定了要用作类别的列（响应变量Species），在~后，我们使用.字符表示所有其他列将作为构建模型的数据。我们将data参数设置为train_set数据框，并为kernel和gamma类型选择合适的值。type可以是基于分类或回归的；kernel是为不同数据和问题设计的多种函数之一；而gamma是核函数的一个参数。你可能希望查看函数文档以获取详细信息，这些值也可以通过经验进行优化。

在第 3 步中，我们创建了一些对象，用来在二维空间中呈现四维边界。首先，我们选择不需要绘制的列（这些列将保持不变），然后使用lapply()函数遍历这些列名的字符向量，应用一个函数计算命名列的均值。我们将列名添加到cols_to_hold变量中的结果列表中。接着，我们使用通用的plot()函数，传入模型、训练数据、绘图的两个维度（通过公式Petal.Width ~ Petal.Length指定），以及一个slice参数，从held_constant列表中提取其他列的均值。

结果如下所示，显示了每个类别的边界颜色：

在第 4 步中，我们使用predict()对测试集进行预测，并通过caret::confusionMatrix()生成混淆矩阵以查看准确性。

在没有先验信息的情况下学习数据中的分组

在生物信息学中，通常需要在没有预先知道分组信息或数量的情况下将事物分类。这一过程通常被称为聚类，是一种无监督的机器学习方法。这种方法常见于基因组学实验，特别是在 RNA 测序及相关表达技术中。在本教程中，我们将从一个包含约 150 个样本的大型基因表达数据集开始，学习如何估算样本的分组数量，并应用一种基于主成分分析（PCA）降维的方法进行聚类，接着使用 k 均值聚类。

准备工作

对于本教程，我们需要factoextra和biobase库（后者来自Bioconductor），以及本书仓库中的datasets/ch1文件夹中的modencodefly_eset.RData文件。

如何操作……

在没有先验信息的情况下，了解数据中的分组可以通过以下步骤完成：

1. 加载数据并运行 PCA：

library(factoextra)
library(Biobase)

load(file.path(getwd(), "datasets", "ch1", "modencodefly_eset.RData") ) 
expr_pca <- prcomp(exprs(modencodefly.eset), scale=TRUE, center=TRUE ) fviz_screeplot(expr_pca)

1. 提取主成分并估算最佳聚类：

main_components <- expr_pca$rotation[, 1:3]
fviz_nbclust(main_components, kmeans, method = "wss")

1. 执行 k 均值聚类并进行可视化：

kmean_clus <- kmeans(main_components, 5, nstart=25, iter.max=1000)

fviz_cluster(kmean_clus, data = main_components,
             palette = RColorBrewer::brewer.pal(5, "Set2"),
             ggtheme = theme_minimal(),
             main = "k-Means Sample Clustering"
             )

它是如何工作的……

在步骤 1中，我们使用load()函数将modencodefly.eset对象导入内存；这是一个基因表达数据集。然后，我们使用Biobase函数exprs()提取表达测量值，并将其传递给prcomp()函数，后者执行 PCA 并返回一个 PCA 对象，我们将其存储在expr_pca变量中。

然后我们使用factoextra函数fviz_screeplot()绘制 PCA 图，并看到以下图示：

这显示了每个主成分所捕获的数据变异度。前三个成分捕获了超过 70%的变异度。因此，我们可以使用这三个成分，而不是整个 150 列的数据集，从而大大简化过程并加速分析。

在步骤 2中，我们通过对子集化expr_pca对象的旋转槽提取主要成分，提取前三列——这些列对应于前三个成分。我们将它们保存在一个名为main_components的变量中，并使用fviz_nbclust()函数对main_components和kmeans函数进行处理，生成以下图示：

在这个函数中，数据被划分为越来越多的聚类，同时计算wss (组内平方和)，这是衡量聚类内变异性的一个指标。图示表明，组内平方和在约 5 个聚类之前大幅下降，之后没有明显改善，表明数据大约包含 5 个聚类。

在步骤 3中，我们使用kmeans()函数执行 k-means 聚类，将main_components作为第一个参数的数据，将5作为聚类数的第二个参数。nstart和iter.max参数的值对于大多数算法运行来说是合理的选择。最后，我们将kmeans_clust对象传递给fviz_cluster()函数，并设置一些显示选项，生成以下图示：

还有更多

我们已经对这个数据集的样本或列进行了 k-means 聚类。如果你希望对基因或行做同样的操作，可以从步骤 2中的未旋转数据的x槽中提取主要成分：

main_components <- expr_pca$x[, 1:3]

如果你希望获取每个样本的实际聚类 ID，可以在kmeans_clus对象的cluster槽中找到：

kmean_clus$cluster[1:5]

## SRX007811 SRX008180 SRX008227 SRX008238 SRX008258 
##         2         2         2         2         2

使用随机森林识别数据中最重要的变量

在本章中，我们已经看到过随机森林算法的应用，在使用随机森林预测类别的示例中，我们用它进行类别预测和回归。这里，我们将用它来做不同的事情——尝试找出数据集中哪些变量对训练模型的分类或回归准确度贡献最大。这只需要对已有代码做一个简单的修改，并使用一个或两个新函数。

准备开始

我们将需要randomForest包和内置的iris数据集。

如何操作...

使用随机森林识别数据中最重要的变量可以通过以下步骤完成：

1. 准备训练数据和测试数据：

library(randomForest)

train_rows <- sample(nrow(iris), 0.8 * nrow(iris), replace = FALSE)
train_set <- iris[train_rows, ]
test_set <- iris[-train_rows, ]

1. 训练模型并创建importance图：

model <- randomForest(Species ~ . , data = train_set, mtry = 2, importance = TRUE)
varImpPlot(model)

它是如何工作的...

在步骤 1中，我们进行类似于之前几个配方中的数据集拆分。使用sample()函数，我们创建了一个包含原始iris数据 80%行号的列表，然后，使用子集和负子集提取这些行。

在步骤 2中，我们使用randomForest()函数训练模型。这里的第一个参数是一个公式；我们指定Species是我们希望预测的值，基于所有其他变量，这些变量由.描述。data是我们的train_set对象。此配方的关键是确保将importance变量设置为TRUE，这意味着模型将测试哪些变量在从模型构建中省略时，会导致准确度的最大下降。一旦模型构建并测试完毕，我们可以使用varImpPlot()函数可视化每个变量的重要性。这样，我们得到以下图表：

我们可以看到，当省略Petal.Width和Petal.Length变量时，会导致模型准确度的最大下降，因此，根据这一标准，它们是最重要的。

使用 PCA 识别数据中最重要的变量

我们已经在无先验信息的数据分组学习配方中看到过 PCA 的应用，它是一种降维技术——一种在保留重要信息的同时减少数据集大小的方法。正如你所想，这意味着我们可以了解哪些原始变量对降维后的表示贡献最大，因此可以确定哪些是最重要的。我们将在本配方中看到这一点。

准备就绪

对于此配方，我们将使用factoextra包和内置的iris数据集。

如何操作...

使用 PCA 识别数据中最重要的变量可以通过以下步骤完成：

1. 执行 PCA：

library(factoextra)
pca_result <- prcomp(iris[,-5], scale=TRUE, center=TRUE )

1. 创建变量图：

fviz_pca_var(pca_result, col.var="cos2")

它是如何工作的...

这个简短的配方在步骤 1开始时简单构建了从prcomp()函数获得的pca_result。我们将iris数据作为第一个参数（不包括第五个分类列），并对数据进行缩放和中心化——这可以避免不同量纲的测量差异占用不当的权重。

构建好pca_result后，我们可以使用fviz_pca_var()函数绘制变量，得到以下图表：

在其中，我们可以看到箭头表示每个变量。箭头从中心移动的角度表示变量的一个特征；箭头之间的距离越近，变量越相似——因此，Petal.Length和Petal.Width是高度相关的变量。箭头的颜色表示一个复杂的量（称为cos2），它代表变量贡献的质量。变量的贡献越高，cos2就越高。在这里，我们可以看到Sepal.Width和Petal.Length对主成分分析（PCA）的贡献较大。Petal.Width由于过于相似，无法被考虑。这与通过随机森林识别数据中最重要的变量的结果不同，因为这两种技术提出的问题不同。

第十章：使用 Tidyverse 和 Bioconductor 编程

R 是一种适合交互式使用的优秀语言；然而，这也意味着许多用户没有体验过将其作为编程语言使用——也就是说，在自动化分析和节省重复工作所需时间和精力方面的应用。在本章中，我们将介绍一些实现这一目标的技术——特别是，我们将探讨如何将基础 R 对象集成到 tidyverse 工作流中，如何扩展 Bioconductor 类以满足我们自己的需求，以及如何使用文献编程和笔记本风格的编码来保持我们的工作记录具有表达性和可读性。

本章将介绍以下食谱：

• 使基础 R 对象整洁

• 使用嵌套数据框

• 为 mutate 编写函数

• 使用编程方式操作 Bioconductor 类

• 开发可重用的工作流和报告

• 使用 apply 系列函数

技术要求

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

以下是你需要的 R 包。通常，你可以使用 install.packages("package_name") 来安装这些包。列在 Bioconductor 下的包需要使用专用的安装程序进行安装。如果需要做其他操作，安装过程将在使用这些包的食谱中描述：

• Bioconductor：

  • Biobase

  • biobroom

  • SummarizedExperiment

• broom

• dplyr

• ggplot2

• knitr

• magrittr

• purrr

• rmarkdown

• tidyr

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

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

使基础 R 对象整洁

tidyverse软件包集（包括dplyr、tidyr和magrittr）通过应用整洁的工作方式，在数据处理和分析中对 R 语言产生了巨大影响。本质上，这意味着数据保持在一种特定的整洁格式中，每一行包含一个单独的观测值，每一列包含一个变量的所有观测值。这样的结构意味着分析步骤有可预测的输入和输出，可以构建成流畅且富有表现力的管道。然而，大多数基础 R 对象并不整洁，往往需要大量的编程工作来提取所需的部分，以便下游使用。在这个教程中，我们将介绍一些函数，用于自动将一些常见的基础 R 对象转换为整洁的数据框。

准备工作

我们将需要tidyr、broom和biobroom包。我们将使用内置的mtcars数据和来自本书仓库datasets/ch1文件夹中的modencodefly_eset.RData。

如何操作...

使基础 R 对象整洁可以通过以下步骤完成：

1. 整理一个lm对象：

library(broom)
model <- lm(mpg ~ cyl + qsec, data = mtcars)
tidy(model) 
augment(model)
glance(model)

1. 整理一个t_test对象：

t_test_result <- t.test(x = rnorm(20), y = rnorm(20) )
tidy(t_test_result)

1. 整理一个 ANOVA 对象：

anova_result <- aov(Petal.Length ~ Species, data = iris)
tidy(anova_result)
post_hoc <- TukeyHSD(anova_result)
tidy(post_hoc)

1. 整理一个Bioconductor ExpressionSet对象：

library(biobroom)
library(Biobase)
load(file.path(getwd(), "datasets", "ch1", "modencodefly_eset.RData") ) 
tidy(modencodefly.eset, addPheno = TRUE)

它是如何工作的...

步骤 1 展示了使用lm()函数整理lm对象的一些函数。第一步是创建对象。在这里，我们使用mtcars数据执行一个多元回归模型。然后，我们对该模型使用tidy()函数，返回模型组件的对象摘要，例如系数，以整洁的数据框形式。augment()函数返回额外的逐观测数据，如果需要的话——同样是整洁格式。glance()函数检查模型本身，并返回关于模型的摘要——自然地，以整洁格式显示。glance()对于比较模型非常有用。

步骤 2 展示了相同的过程，适用于t.test对象。首先，我们对两个随机数向量进行 t 检验。tidy()函数会将所有细节以整洁的数据框形式返回。

在步骤 3中，我们对iris数据运行方差分析（ANOVA）。我们使用aov()函数查看Species对Petal.Length的影响。我们可以再次对结果使用tidy()，但它会给出模型组件的摘要。实际上，我们更感兴趣的可能是来自后验检验的比较结果，该检验使用TukeyHSD()函数进行，它也可以在tidy()中使用。

在步骤 4中，我们使用biobroom版本的tidy()对ExpressionSet对象进行处理。这将表达值的方阵转化为整洁的数据框，并包含样本和其他数据类型的列。额外的参数addPheno是特定于此类对象的，并将ExpressionSet元数据容器中的表型元数据插入其中。请注意，结果数据框超过 200 万行——生物学数据集可能很大，并且生成非常大的数据框。

使用嵌套数据框

数据框（dataframe）是整洁工作方式的核心，我们通常将其视为一个类似电子表格的矩形数据容器，每个单元格中仅包含一个值。实际上，数据框可以嵌套——即它们可以在特定的单个单元格中包含其他数据框。这是通过将数据框的向量列替换为列表列来实现的。每个单元格实际上是列表的一个成员，因此任何类型的对象都可以保存在外部数据框的概念性单元格中。在本教程中，我们将介绍创建嵌套数据框的方式以及与之合作的不同方法。

准备工作

我们将需要tidyr、dplyr、purrr和magrittr库。我们还将使用来自ggplot2包的diamonds数据，但我们不会使用任何函数。

工作原理...

使用嵌套数据框可以通过以下步骤完成：

1. 创建一个嵌套数据框：

library(tidyr)
library(dplyr)
library(purrr)
library(magrittr)
library(ggplot2)

nested_mt <- nest(mtcars, -cyl)

1. 添加一个新列表列，包含lm()的结果：

nested_mt_list_cols <- nested_mt %>% mutate(
 model = map(data, ~ lm(mpg ~ wt, data = .x))
)

1. 添加一个新列表列，包含tidy()的结果：

nested_mt_list_cols <- nested_mt_list_cols %>% mutate(
 tidy_model = map(model, tidy)
)

1. 解开整个数据框：

unnest(nested_mt_list_cols, tidy_model)

1. 在单个步骤中运行管道：

models_df <- nest(mtcars, -cyl) %>%
 mutate(
 model = map(data, ~ lm(mpg ~ wt, data = .x)),
 tidy_model = map(model, tidy)
 ) %>%
 unnest(tidy_model)

工作原理...

在步骤 1中，我们使用nest()函数将mtcars数据框进行嵌套。-选项告诉函数哪些列在嵌套时要排除；将cyl列转换为因子，用于创建不同的子集。从概念上讲，这类似于dplyr::group_by()函数。检查该对象会得到以下内容：

 A tibble: 3 x 2
## cyl   data 
## <dbl> <list> 
## 1 6   <tibble [7 × 10]> 
## 2 4   <tibble [11 × 10]>
## 3 8   <tibble [14 × 10]>

嵌套数据框包含一个名为data的新数据框列，以及被简化的cyl列。

在步骤 2中，我们通过使用mutate()在数据框中创建一个新列。在此过程中，我们使用purrr包中的map()函数，它遍历作为第一个参数提供的列表项（即我们的数据框列），并在作为第二个参数提供的代码中使用它们。在这里，我们对嵌套数据使用lm()函数，一次处理一个元素——注意，.x变量表示我们当前正在处理的内容——也就是列表中的当前项。运行后，列表看起来是这样的：

##  cyl   data              model 
##  <dbl> <list> <list>
## 1 6   <tibble [7 × 10]>  <lm> 
## 2 4   <tibble [11 × 10]> <lm> 
## 3 8   <tibble [14 × 10]> <lm>

新的model列表列包含我们的lm对象。

在确认添加新列表列的模式是使用mutate()并在其中嵌入map()后，我们可以同样整理lm对象。这就是在步骤 3中发生的情况。结果给我们带来了以下嵌套数据框：

##  cyl   data              model  tidy_model 
## <dbl> <list>            <list>  <list> 
## 1 6   <tibble [7 × 10]>  <lm>   <tibble [2 × 5]>
## 2 4   <tibble [11 × 10]> <lm>   <tibble [2 × 5]>
## 3 8   <tibble [14 × 10]> <lm>   <tibble [2 × 5]>

步骤 4使用unnest()函数将所有内容恢复为单个数据框；第二个参数tidy_model是需要解包的列。

步骤 5将步骤 1到步骤 4的整个过程合并为一个管道，突出显示这些只是常规的tidyverse函数，并且可以链式调用，无需保存中间步骤。

还有更多内容...

unnest()函数只有在嵌套列表列成员兼容且可以根据正常规则合理对齐和回收时才有效。在许多情况下，这种情况并不成立，因此你需要手动操作输出。以下示例展示了我们如何做到这一点。工作流程基本上与之前的示例相同，唯一的变化是在早期步骤中，我们使用dplyr::group_by()来创建nest()的分组。在mutate()中，我们传递自定义函数来分析数据，但其他步骤保持不变。最后一步是最大的变化，利用transmute()来删除不需要的列，并创建一个新的列，它是map_dbl()和自定义汇总函数的结果。map_dbl()类似于map()，但只返回双精度数值向量。其他的map_**函数也存在。

为dplyr::mutate()编写函数

dplyr的mutate()函数非常有用，可以根据现有列的计算结果向数据框中添加新列。它是一个矢量化函数，但通常被误解为按行工作，实际上它是按列工作的，也就是说，它对整个向量进行操作，并利用 R 内置的回收机制。这种行为往往会让那些想在复杂例子或自定义函数中使用mutate()的人感到困惑，因此，在本教程中，我们将探讨mutate()在某些情况下的实际行为，希望能够带来启示。

准备工作

为此，我们需要dplyr包和内置的iris数据。

如何做……

为dplyr::mutate()编写函数可以通过以下步骤完成：

1. 使用返回单个值的函数：

return_single_value <- function(x){
 sum(x) 
}
iris %>% mutate(
 result = return_single_value(Petal.Length)
)

1. 使用返回与给定值相同数量值的函数：

return_length_values <- function(x){
 paste0("result_", 1:length(x))
}
iris %>% mutate(
 result = return_length_values(Petal.Length)
)

1. 使用返回既不是单个值也不是与给定值相同数量值的函数：

return_three_values <- function(x){
 c("A", "b", "C")
}
iris %>% mutate(
 result = return_three_values(Petal.Length)
)

1. 强制函数的重复以适应向量的长度：

rep_until <- function(x){
 rep(c("A", "b", "C"), length.out = length(x))
}
iris %>% mutate(
 result = rep_until(Petal.Length)
)

它是如何工作的……

在步骤 1中，我们创建一个函数，给定一个向量，返回一个单一的值（长度为 1 的向量）。然后我们在mutate()中使用它，添加一个名为result的列，并得到如下结果：

## Sepal.Length Sepal.Width Petal.Length Petal.Width Species result 
## 1 5.1 3.5 1.4 0.2 setosa 563.7 
## 2 4.9 3.0 1.4 0.2 setosa 563.7 
## 3 4.7 3.2 1.3 0.2 setosa 563.7 
## 4 4.6 3.1 1.5 0.2 setosa 563.7

请注意，函数在result列中返回的单个值是反复出现的。对于length == 1的向量，R 会回收结果并将其放置到每个位置。

在步骤 2中，我们走到对立面，创建一个函数，给定一个向量，返回一个相同长度的向量（具体来说，它返回一个将result_与向量中的位置数字拼接在一起的向量）。当我们运行它时，我们得到如下结果：

## Sepal.Length Sepal.Width Petal.Length Petal.Width Species result 
## 1 5.1 3.5 1.4 0.2 setosa result_1 
## 2 4.9 3.0 1.4 0.2 setosa result_2 
## 3 4.7 3.2 1.3 0.2 setosa result_3 
## 4 4.6 3.1 1.5 0.2 setosa result_4

由于它的长度与数据框中其余列完全相同，R 会接受它并将其作为新列应用。

在步骤 3中，我们创建一个返回三个元素的向量的函数。由于长度既不是 1，也不是数据框中其他列的长度，代码会失败。

在第 4 步中，我们将探讨如何重复一个不兼容长度的向量，以便在需要时使其适应。rep_until()函数与length.out参数会重复其输入，直到向量的长度为length.out。通过这种方式，我们可以得到以下列，这就是我们在第 3 步中使用该函数时所期望看到的结果：

## Sepal.Length Sepal.Width Petal.Length Petal.Width Species result 
## 1 5.1 3.5 1.4 0.2 setosa A 
## 2 4.9 3.0 1.4 0.2 setosa b 
## 3 4.7 3.2 1.3 0.2 setosa C 
## 4 4.6 3.1 1.5 0.2 setosa A 
## 5 5.0 3.6 1.4 0.2 setosa b

使用 Bioconductor 类进行编程操作

Bioconductor的广泛应用意味着有大量的类和方法可以完成几乎任何你想要的生物信息学工作流程。不过，有时候，额外的数据槽或一些工具上的调整会帮助简化我们的工作。在本教程中，我们将探讨如何扩展一个现有的类，以包含一些特定于我们数据的额外信息。我们将扩展SummarizedExperiment类，以添加假设的条形码信息——一种元数据，表示一些核苷酸标签，这些标签可以标识包含在序列读取中的样本。

准备工作

对于本教程，我们只需要Bioconductor的SummarizedExperiment包。

如何操作...

使用以下步骤可以通过编程方式与Bioconductor类进行交互：

1. 创建一个继承自SummarizedExperiment的新类：

setClass("BarcodedSummarizedExperiment",
   contains = "SummarizedExperiment",
   slots = c(barcode_id = "character", barcode_sequence = "character")
 )

1. 创建构造函数：

BarcodedSummarizedExperiment <- function(assays, rowRanges, colData, barcode_id, barcode_sequence){
   new("BarcodedSummarizedExperiment", 
       SummarizedExperiment(assays=assays, rowRanges=rowRanges, colData=colData),
       barcode_id = barcode_id,
       barcode_sequence = barcode_sequence
   )
}

1. 向类中添加必需的方法：

setGeneric("barcode_id", function(x) standardGeneric("barcode_id"))
setMethod("barcode_id", "BarcodedSummarizedExperiment", function(x) x@barcode_id )

1. 构建新类的实例：

nrows <- 200
ncols <- 6
counts <- matrix(runif(nrows * ncols, 1, 1e4), nrows)
assays <- list(counts = counts)
rowRanges <- GRanges(    rep(c("chr1", "chr2"), c(50, 150)),
                         IRanges(floor(runif(200, 1e5, 1e6)), width=100),
                         strand=sample(c("+", "-"), 200, TRUE),
                         feature_id=sprintf("ID%03d", 1:200)
)
colData <- DataFrame(
                Treatment=rep(c("ChIP", "Input"), 3),
                row.names=LETTERS[1:6]
)

my_new_barcoded_experiment <- BarcodedSummarizedExperiment(
        assays = assays, 
        rowRanges = rowRanges, 
        colData = colData, 
        barcode_id = letters[1:6], 
        barcode_sequence = c("AT", "GC", "TA", "CG","GA", "TC") 
)

1. 调用新的方法：

barcode_id(my_new_barcoded_experiment)

它是如何工作的...

在第 1 步中，我们使用setClass()函数创建一个新的 S4 类。这个函数的第一个参数是新类的名称。contains参数指定我们希望继承的现有类（这样我们的新类就会包含这个类的所有功能以及我们创建的任何新功能）。slots参数指定我们希望添加的新数据槽，并要求我们为其指定类型。在这里，我们添加了文本数据槽，用于新的barcode_id和barcode_sequence槽，因此我们为这两个槽都使用character类型。

在第 2 步中，我们创建一个构造函数。该函数的名称必须与类名相同，并且我们在调用function()时指定创建新对象所需的参数。在函数体内，我们使用new()函数，其第一个参数是要实例化的类的名称。其余的部分用于填充实例数据；我们调用继承自SummarizedExperiment的构造函数，以填充新对象的那一部分，然后手动填充新的条形码槽。每次运行BarcodedSummarizedExperiment时，我们都会获得该类的一个新对象。

在步骤 3中，我们添加了一个新函数（严格来说，在 R 中，它被称为方法）。如果我们选择一个在 R 中尚不存在的Generic函数名称，我们必须使用setGeneric()注册该函数名，setGeneric()的第一个参数是函数名，第二个参数是一个模板函数。Generic函数设置完毕后，我们可以使用setMethod()函数添加实际的函数。新函数的名称是第一个参数，它将附加到的类是第二个参数，而代码本身是第三个参数。请注意，我们这里只是在创建一个访问器（getter）函数，它返回当前对象中barcode_id槽的数据。

在步骤 4中，我们的准备工作已经完成，可以构建类的实例。在这一步的前六行中，我们仅仅创建了构建对象所需的数据。这部分数据进入一个普通的SummarizedExperiment对象；你可以在文档中看到更多关于这里具体发生了什么的细节。然后，我们可以通过调用BarcodedSummarizedExperiment函数，并传入我们创建的数据以及新barcode_id和barcode_sequence槽的特定数据，实际创建my_new_barcoded_experiment。

现在，创建了对象后，在步骤 5中，我们可以像调用其他函数一样使用我们的方法，传入我们新创建的对象作为参数。

开发可重用的工作流程和报告

在生物信息学中，一个非常常见的任务是撰写结果，以便与同事交流，或者仅仅为了在实验记录中留有一份完整的记录（无论是电子版还是纸质版）。一个关键技能是尽可能让工作可重复，以便当我们需要回顾它时，或者有其他人对我们所做的工作感兴趣时，他们可以复制整个过程。一个日益流行的解决方案是使用文献编程技术和可执行笔记本，这些笔记本结合了可读的文本、分析代码和计算输出，打包成一个单一文档。在 R 中，rmarkdown包使我们能够以这种方式将代码和文本结合在一起，并生成多种格式的输出文档。

在这个食谱中，我们将查看一个可以通过rmarkdown编译的文档的大规模结构。RStudio 应用程序使得这个过程非常简便，所以我们将通过这个工具来查看编译过程。

准备工作

对于这个食谱，我们需要RStudio 应用程序和rmarkdown包。这个食谱的示例代码可以在本书的datasets/ch10/文件夹中的example_rmarkdown.rmd文件中找到。

如何做到这一点...

开发可重用的工作流程和报告可以通过以下步骤完成：

1. 在外部文件中，添加一个YAML头部：

---
title: "R Markdown Report"
author: "R Bioinformatics Cookbook"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
 html_document:
 df_print: paged
 bookdown::html_document2:
 fig_caption: yes
 keep_md: yes
 toc: yes
---

1. 然后，添加一些文本和代码进行解释：

We can include text and create code blocks, the code gets executed and the result passed in

```{r}

x <- iris$Sepal.Width

y <- iris$Sepal.Length

lm(y ~ x, data = iris)

```py

1. 文本可以使用最小化标记格式进行格式化：

## We can format text using Markdown
We can create many text formats including *italics* and **bold**,
We can make lists 
 1\. First item
 2\. Second item

1. 在块内应用更多选项并传递变量：

The whole document acts as a single R session - so variables created in earlier blocks can still be used later.
Plots are presented within the document. Options for blocks can be set in the header

```{r, fig.width=12 }

plot(x, y)

```py

它是如何工作的…

这里的代码很特殊，因为它必须从外部文档中运行；在 R 控制台中无法运行。运行文档的编译步骤有几种方法。在 RStudio 中，一旦安装了rmarkdown并且正在编辑一个.Rmd扩展名的文档，您会看到一个knit按钮。或者，您可以通过控制台使用rmarkdown::render()函数编译文档，尽管我建议使用 RStudio IDE 来完成此操作。

在步骤 1中，我们创建了一个YAML头部，描述了文档的渲染方式，包括输出格式、动态日期插入、作者和标题。这些内容将自动添加到您的文档中。

在步骤 2中，我们实际上创建了一些内容——第一行是纯文本，将作为段落文本传递到最终文档中，不做任何修改。块中的部分由py` ```py is code to be interpreted. Options for the block go inside the curly brackets—here, {r}` means this should be an R code block (some other languages are supported too). The code in this block is run in a new R session, its output captured; and inserted immediately after the code block.

In Step 3, we create some plaintext with the Markdown tags. ## gives us a line with a second-level heading, the **starred** text gives us different formatting options, and we can also create lists. Valid Markdown is interpreted and the reformatted text is passed into the eventual document.

In Step 4, we start with some more plaintext and follow with a new code block. The options for the code block are set in the curly brackets again—here, we set a width for figures in the plot. Note that the code in this block refers to variables created in an earlier block. Although the document creates a new R session without access to variables already in the usual console, the document itself is a single session so blocks can access earlier block's variables, allowing the code and text to be mixed up at whatever resolution the author requires. Finally, the resulting figure is inserted into the document just like code. 

Making use of the apply family of functions

Programming in R can sometimes seem a bit tricky; the control flow and looping structures it has, are a bit more basic than in other languages. As many R functions are vectorized, the language actually has some features and functions; that mean we don't need to take the same low-level approach we may have learned in Python or other places. Instead, base R provides the apply functions to do the job of common looping tasks. These functions all have a loop inside them, meaning we don't need to specify the loop manually. In this recipe, we'll look at using some apply family functions with common data structures to loop over them and get a result. The common thread in all of the apply functions is that we have an input data structure that we're going to iterate over and some code (often wrapped in a function definition) that we're going to apply to each item of the structure.

Getting ready

We will only need base R functions and data for this recipe, so you are good to go!

How to do it...

Making use of the apply family of functions can be done using the following steps:

1. Create a matrix and use apply to work on it:

m <- matrix(rep(1:10, 10, replace = TRUE), nrow = 10)

apply(m, 1, sum)

apply(m, 2, sum)

```py

2.  Use `lapply` over the vector:

numbers <- 1:3

number_of_numbers <- function(x){

rnorm(x)

}

my_list <- lapply(numbers, number_of_numbers)

3.  Use `lapply` and `sapply` over the list:

summary_function <- function(x){

mean(x)

}

lapply(my_list, summary_function)

sapply(my_list, summary_function)

4.  Use `lapply` over a dataframe:

list_from_data_frame <- lapply(iris, mean, trim = 0.1, na.rm = TRUE )

unlist(list_from_data_frame)

# How it works...

*Step 1* begins with the creation of a 10 x 10 matrix, with rows holding the same number and columns running from 1 to 10\. Inspecting it makes it clear, as partly shown in the following output:

> m

[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]

[1,] 1 1 1 1 1 1 1 1 1 1

[2,] 2 2 2 2 2 2 2 2 2 2

[3,] 3 3 3 3 3 3 3 3 3 3

We then use `apply()`: the first argument is the object to loop over, the second is the direction to loop in (or margin, 1 = rows, and 2 = columns), and the third is the code to apply. Here, it's the name of a built-in function, but it could be a custom one. Note it's the margin argument that affects the amount of data that is taken each time. Contrast the two `apply()` calls:

apply(m, 1, sum)

[1] 10 20 30 40 50 60 70 80 90 100

apply(m, 2, sum)

[1] 55 55 55 55 55 55 55 55 55 55

Clearly, `margin = 1` is taking each row at a time, whereas `margin = 2` is taking the columns. In any case, `apply()` returns a vector of results, meaning the results must be of the same type each time. It is not the same shape as the input data.

With *Step 2*, we move onto using `lapply()`, which can loop over many types of data structures, but always returns a list with one member for each iteration. Because it's a list, each member can be of a different type. We start by creating a simple vector containing the integers 1 to 3 and a custom function that just creates a vector of random numbers of a given length. Then, we use `lapply()` to apply that function over the vector; the first argument to `lapply()` is the thing to iterate over, and the second is the code to apply. Note that the current value of the vector we're looping over is passed automatically to the called function as the argument. Inspecting the resulting list, we see the following:

my_list

[[1]] [1] -0.3069078

[[2]] [1] 0.9207697 1.8198781

[[3]] [1] 0.3801964 -1.3022340 -0.8660626

We get a list of one random number, then two, then three, reflecting the change in the original vector.

In *Step 3*, we see the difference between `lapply()` and `sapply()` when running over the same object. Recall `lapply()` always returns a list but `sapply()` can return a vector (`s` can be thought of as standing for *simplify*). We create a simple summary function to ensure we only get a single value back and `sapply()` can be used. Inspecting the results, we see the following:

lapply(my_list, summary_function)

[[1]] [1] -0.3069078

[[2]] [1] 1.370324

[[3]] [1] -0.5960334

sapply(my_list, summary_function)

[1] -0.3069078 1.3703239 -0.5960334

Finally, in *Step 4*, we use `lapply()` over a dataframe, namely, the built-in `iris` data. By default, it applies to columns on a dataframe, applying the `mean()` function to each one in turn. Note the last two arguments (`trim` and `na.rm`) are not arguments for `lapply()`, though, it does look like it. In all of these functions, the arguments after the vector to iterate over and the code (in other words, argument positions 1 and 2) are all passed to the code being run—here, our `mean()` function. The column names of the dataframe are used as the member names for the list. You may recall that one of the columns in `iris` is categorical, so `mean()` doesn't make much sense. Inspect the result to see what `lapply()` has done in this case:

lapply(iris, mean, trim = 0.1, na.rm = TRUE )

$Sepal.Length [1] 5.808333

$Sepal.Width [1] 3.043333

$Petal.Length [1] 3.76

$Petal.Width [1] 1.184167

$Species [1] NA

It has returned `NA`. Also, it has generated a warning but not failed. This can be a source of bugs in later analyses.

With a simple list like this, we can also use `unlist()` to get a vector of the results:

unlist(list_from_data_frame)

Sepal.Length Sepal.Width Petal.Length Petal.Width Species

5.808333 3.043333 3.760000 1.184167 NA

如果存在名称，向量将被命名。


# 第十一章：为了代码重用构建对象和包

在最后一章中，我们将探讨如何将我们的代码从自己的机器中带出去，并与世界分享。我们最常分享的对象将是我们自己！因此，为了使我们的编程生活更加轻松和流畅，我们将学习如何创建对象和类来简化我们的工作流程，并如何将它们打包成可在其他项目中重用的包。我们将查看如何在 GitHub 等网站上共享代码的工具，以及如何检查代码是否按预期工作。

本章将涵盖以下食谱：

+   创建简单的 S3 对象来简化代码

+   利用 S3 类的通用对象函数

+   使用 S4 系统创建结构化和正式的对象

+   简单的代码打包方法，用于共享和重用

+   使用 `devtools` 从 GitHub 托管代码

+   构建单元测试套件，以确保函数按预期工作

+   使用 Travis 进行持续集成，以确保代码经过测试并保持最新

# 技术要求

您需要的示例数据可以从本书的 GitHub 仓库获取：[`github.com/PacktPublishing/R-Bioinformatics-Cookbook`](https://github.com/PacktPublishing/R-Bioinformatics-Cookbook)[。](https://github.com/danmaclean/R_Bioinformatics_Cookbook) 如果您想按照书中的代码示例使用，您需要确保这些数据位于您的工作目录的子目录中。

这里是您需要的 R 包。通常，您可以通过 `install.packages("package_name")` 来安装这些包。列在 `Bioconductor` 下的包需要通过专用的安装程序安装。如果您需要做其他安装操作，会在包使用的食谱中描述：

+   `devtools`

+   `usethis`

对于一些后续的食谱，我们还需要安装 `git` 版本控制系统。请访问官方网页以获取适用于您系统的最新版本：[`git-scm.com/downloads`](https://git-scm.com/downloads)。您还可以在 GitHub 网站上找到一个有用的 GitHub 账户。如果您还没有 GitHub 账户，请访问[`github.com/`](https://github.com/)。

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

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

`letters[1:5]` 这将给我们以下输出：

`## a b c d e` 请注意，输出行前面有`##`。

# 创建简单的 S3 对象以简化代码

创建你自己的对象可以大大简化代码和工作流，使它们更容易被你复现和重用，并且将程序的内部逻辑抽象掉，减少你作为程序员的认知负担，让你能够更多地集中精力在项目的生物信息学和分析方面。R 实际上有多种方法来创建对象和类。在这个教程中，我们将重点看它最简单、最临时的方法——S3。这是一种相当非正式的创建对象和类的方式，但在很多情况下都足够使用。

# 准备工作

在这个教程中，我们只需要基础的 R 函数，因此无需安装任何额外的工具。

# 如何实现...

创建简单的 S3 对象以简化代码可以通过以下步骤完成：

1.  创建一个构造函数：

```py
SimpleGenome <- function( nchr=NA, lengths = NA){

 genome <- list(
 chromosome_count = nchr,
 chromosome_lengths = lengths
 )
 class(genome) <- append(class(genome), "SimpleGenome")
 return(genome)
}

1. 调用构造函数来创建新对象：

ecoli <- SimpleGenome(nchr = 1, lengths = c(4600000) )
bakers_yeast <- SimpleGenome(nchr = 1, lengths=c(12100000))

它是如何工作的...

第 1 步 是所有工作的起点。这就是我们创建 S3 对象所需的全部内容。如你所见，这是非常简洁的代码。我们只是创建了一个生成并返回数据结构的函数。我们的类应该代表一个简化的基因组，并且我们希望它能保存一些关于基因组的基本信息。SimpleGenome() 函数是我们创建对象的构造函数。SimpleGenome 创建的基因组列表是构成最终对象主体的数据结构。这个列表的成员就是对象的槽，因此我们创建了名为 chromosome_count 和 chromosome_length 的成员，以表示基因组的一些特征。完成这一步后，我们进行一个重要的步骤——我们将类名（SimpleGenome）附加到基因组列表的类属性上。正是这一点使 R 能够识别这个对象属于 SimpleGenome 类。现在我们可以返回创建的 S3 对象。

在第 2 步中，我们只是简单地使用构造函数来创建类的实例。检查得到的对象如下所示：

> ecoli 
$chromosome_count 
[1] 1 
$chromosome_lengths 
[1] 4600000 
attr(,"class") 
[1] "list" "SimpleGenome"

> bakers_yeast 
$chromosome_count 
[1] 1 
$chromosome_lengths 
[1] 12100000 
attr(,"class") 
[1] "list" "SimpleGenome" 

我们可以看到对象的槽、对象之间的差异，以及包含新 SimpleGenome 对象的类。这就是我们创建 S3 对象的方式；它是一个简单但有效的做法。与仅仅创建一个普通数据结构（如列表）相比，优势并不是立刻显现出来，但当我们查看如何在下一个教程中创建方法时，原因将会更加明确。

利用 S3 类的通用对象函数

一旦我们有了一个 S3 对象，我们需要创建与之配套的函数。这些函数实际上是使得长期使用这些对象变得更容易的关键。正是在这些函数中，我们可以抽象化对象数据的处理，从而减少每次使用时所需的工作量。R 的对象系统基于通用函数。这些函数是具有相同基础名称但类特定名称扩展的分组函数。每个分组称为方法，R 会根据方法调用的对象的类来决定调用属于该方法的具体哪个函数。这意味着我们可以在A类的对象上调用plot()，并得到与在B类对象上调用时完全不同的图形。在本食谱中，我们将了解它是如何工作的。

准备工作

对于本食谱，我们将使用基础的 R 函数，因此无需安装任何包，但我们将使用内置的iris数据。

如何实现...

利用通用对象函数与 S3 类一起使用，可以通过以下步骤完成：

1. 在plot()方法中创建一个通用函数：

plot.SimpleGenome <- function(x){
 barplot(x$chromosome_lengths, main = "Chromosome Lengths")
}

1. 创建一个对象并在plot()中使用它：

athal <- SimpleGenome(nchr = 5, lengths = c(34964571, 22037565, 25499034, 20862711, 31270811 ) )
plot(athal)

1. 首先创建一个新方法：

genomeLength <- function(x){
 UseMethod("genomeLength", x)
}

genomeLength.SimpleGenome <- function(x){
 return(sum(x$chromosome_lengths))
}
genomeLength(athal)

1. 修改现有对象的类：

some_data <- iris
summary(some_data)
class(some_data) <- c("my_new_class", class(some_data) )
class(some_data)

1. 为新类创建通用函数：

summary.my_new_class <- function(x){
 col_types <- sapply(x, class)
 return(paste0("object contains ", length(col_types), " columns of classes:", paste (col_types, sep =",", collapse = "," )))
}
summary(some_data)

它是如何工作的...

在第 1 步中，我们创建了一个名为plot.SimpleGenome()的通用函数。这里的特殊命名约定标志着该函数是属于专门针对SimpleGenome类对象的通用绘图函数组。约定是method.class。这就是使通用绘图方法工作的所有必要步骤。

在第 2 步中，我们实际上创建一个SimpleGenome对象，就像在本章中的创建简单的 S3 对象以简化代码的食谱中所做的那样（你需要确保该食谱的第 1 步在当前会话中已经执行，否则这一步无法正常工作），然后调用plot()方法。plot方法查找适用于SimpleGenome对象的通用函数并运行该对象，从而生成我们预期的条形图，如下图所示：

在第 3 步中，我们深入了一些。在这一步中，我们想使用一个不存在的方法名（你可以使用methods()函数查看已经存在的方法），因此我们必须首先创建方法组。我们通过创建一个调用UseMethod()函数的函数来实现这一点，方法的名称作为封闭函数名称和第一个参数。完成后，我们就可以为SimpleGenome类创建通用函数，并通过简单地调用genomeLength()在对象上使用它。由于我们的通用函数只是将chromosome_lengths向量相加，我们得到如下结果：

> genomeLength(athal)
[1] 134634692

第 4 步展示了类查找系统的机制。我们首先复制iris数据，然后在其上使用summary()方法，得到数据框的标准结果：

> summary(some_data) 
Sepal.Length Sepal.Width Petal.Length Petal.Width Species 
Min. :4.300 Min. :2.000 Min. :1.000 Min. :0.100 setosa :50 
1st Qu.:5.100 1st Qu.:2.800 1st Qu.:1.600 1st Qu.:0.300 versicolor:50

接下来，在第 4 步中，我们使用class()函数为some_data对象添加了一个新类。注意，我们将其作为向量的第一个元素添加。我们可以看到，data.frame类仍然存在，但排在我们添加的类之后：

> class(some_data) 
[1] "my_new_class" "data.frame"

然后，在第 5 步中，我们为my_new_class创建了一个通用的summary()函数，使其返回一个完全不同类型的摘要。我们可以看到，当我们调用它时：

> summary(some_data) 
[1] "object contains 5 columns of classes:numeric,numeric,numeric,numeric,factor"

需要注意的一点是，尽管对象有多个类，但默认情况下会选择第一个与类匹配的通用函数。如果你想测试这一点，可以尝试交换class属性的顺序。

使用 S4 系统创建结构化和正式的对象

S4 是 S3 的更正式的对应版本，特别是因为它具有正式的类定义，所以不能随意使用，但它的工作方式与 S3 非常相似，因此我们已经学到的内容通常适用。在这个教程中，我们将快速介绍如何使用 S4 系统创建一个类似于本章前两节中SimpleGenome对象的类。了解 S4 会对你扩展Bioconductor（因为它是用 S4 编写的）有所帮助。

准备工作

再次，我们只使用基础 R，因此无需安装任何东西。

如何做到这一点...

使用 S4 系统创建结构化和正式的对象可以按照以下步骤进行：

1. 编写类定义：

S4genome <- setClass("S4genome", slots = list(chromosome_count = "numeric", chromosome_lengths = "numeric" ))

1. 创建通用函数：

setGeneric( "chromosome_count", 
 function(x){ standardGeneric("chromosome_count") }
)

1. 创建方法：

setMethod( "chromosome_count", "S4genome", function(x){ slot(x, "chromosome_count")} )

它是如何工作的

这里的大纲与前两个教程非常相似。在第 1 步中，我们使用setClass()函数创建了类定义；第一个参数是类的名称，slots参数是一个适当的插槽名称列表以及每个插槽的类型。S4 类需要定义类型。使用中的对象可以像 S3 那样实例化：

## > ecoli <- S4genome(chromosome_count = 1, chromosome_lengths = c(4600000) ) 
## > ecoli An object of class "S4genome" 
## Slot "chromosome_count": [1] 1 
## Slot "chromosome_lengths": [1] 4600000 

在第 2 步中，我们使用setGeneric()函数创建了一个通用函数chromosome_count，传入函数的名称和一个调用standardGeneric()函数的函数。这基本上是模板代码，因此请按照这个步骤进行，并在需要更多细节时查阅文档。

在第 3 步中，我们创建了方法。我们使用setMethod()函数创建了一个chromosome_count方法。第二个参数是这个方法将被调用的类，最后，我们传递了我们想要的代码。匿名函数仅仅调用了传递给它的对象上的slot()函数。slot()返回第二个参数中指定的槽的内容。

另请参见

如果你希望进一步了解 S4 来扩展 Bioconductor 类，请查看 Bioconductor 自己提供的教程：www.bioconductor.org/help/course-materials/2017/Zurich/S4-classes-and-methods.html。

分享和重用代码的简单方法

不可避免地，总会有这么一个时刻，你希望能够重用一些函数或类，而不必每次都重新输入（或者——更糟——复制粘贴）它们。将功能的唯一可靠版本放在一个地方，可以轻松管理，及时发现并解决代码中的错误和变化。因此，在这个教程中，我们将探讨两种简单的封装代码以便重用的方法。我们将讨论包创建的基础知识，尽管我们将创建的包将非常简陋，并且在考虑发布之前，需要进行大量的完善——特别是文档和测试方面。然而，以这种方式创建的包，将在你开发代码时提供帮助。

准备工作

为此，我们需要devtools和usethis包，以及本书仓库中datasets/ch11文件夹下的源代码文件my_source_file.R。

如何操作...

封装代码以便共享和重用可以通过以下步骤完成：

1. 加载现有的源代码文件：

source(file.path(getwd(), "datasets", "ch11", "my_source_file.R"))
my_sourced_function()

1. 创建一个包骨架：

usethis::create_package("newpackage")

1. 编写代码：

my_package_function <- function(x){
 return( c("I come from a package!") )
}

1. 将包代码加载到内存中：

devtools::load_all()

1. 将包安装到当前的 R 安装环境中：

devtools::install()
library(newpackage)

它是如何工作的...

这段代码的第一步展示了一种非常有效但非常基础的加载外部代码的方法。我们使用 source() 函数将 R 代码文件加载到当前的命名空间中。这里的特定文件包含普通的 R 函数，除此之外没有其他内容。source() 函数仅仅是读取外部文件中的代码并将其执行，就像它是直接在当前控制台中输入的一样。由于文件中仅包含函数，因此你必须将这些函数加载到内存中以便立即使用。

步骤 2 更进一步，使用 usethis::create_package() 函数创建了一个简陋的包。该函数会创建一个你提供名称的新文件夹（在这个例子中是 newpackage），并将包所需的所有基本文件和文件夹放入其中。现在你可以在包的 R/ 子文件夹中填充 R 代码，这些代码将在你加载包时加载。尝试使用 步骤 3 中的函数；将此函数添加到 R/ 文件夹中的 my_functions.R 文件中。R/ 文件夹中的文件名并不太重要，你可以有很多文件——但一定要确保它们以 .R 结尾。

步骤 4 将会使用 devtools::load_all() 函数加载你的源代码包到内存中。这大致模拟了我们调用 library() 函数时发生的情况，但并没有真正安装包。通过使用 devtools::load_all()，我们可以快速加载代码进行测试，而不必先安装它，这样如果我们需要更改代码，就不会有破损的版本被安装。我们没有提供任何参数，因此它会加载当前目录中的包（如果提供路径作为第一个参数，它会加载路径中找到的包）。

在 第 5 步 中，我们实际上将代码正确安装到 R 中。我们使用 devtools::install() 函数，它会构建包并将构建后的版本复制到 R 中的常规位置。现在，我们可以像加载任何其他包一样加载该构建版本 (newpackage)。请注意，这意味着我们有两个包的副本—一个是我们安装的版本，另一个是我们正在开发的版本。在开发更多代码并将其添加到包中时，你需要根据需要重复步骤四和五。

使用 devtools 从 GitHub 托管代码

良好的代码开发实践意味着将代码保存在某种版本控制系统中。Git 和 Git 共享网站 GitHub 是其中一种非常流行的系统。在本教程中，我们将使用 usethis 包添加一些有用的非代码文件，帮助描述其他用户如何重用我们的代码以及当前开发状态，并添加机制确保下游用户拥有与你的代码相关依赖的其他包。接下来，我们将看看如何将包上传到 GitHub 以及如何直接从 GitHub 安装。

准备工作

我们需要 usethis 和 devtools 包。

如何操作...

使用 devtools 从 GitHub 托管代码可以通过以下步骤完成：

1. 向包中添加一些有用的元数据和许可证文件：

usethis::use_mit_license(name = "Dan MacLean")
usethis::use_readme_rmd()
usethis::use_lifecycle_badge("Experimental")
usethis::use_version()

1. 向依赖项列表中添加将在安装包时自动安装的内容：

usethis::use_package("ggplot2")

1. 自动设置本地 Git 仓库并获取 GitHub 凭证：

usethis::use_git()
usethis::browse_github_token() 
usethis::use_github()

1. 从 GitHub 安装包：

devtools::install_github("user/repo")

工作原理...

第 1 步 中的代码非常简单，但它为包添加了很多内容。usethis::use_mit_license() 函数会添加一个名为 LICENSE 的文本文件，文件中包含 MIT 许可协议的内容。没有许可证文件，其他人很难知道在什么条件下可以使用该软件。MIT 许可协议是一种简单且宽松的许可，非常适合一般的开源软件，但也有其他替代方案；有关如何选择适合你的许可证，参见此网站：choosealicense.com/。查看 usethis 文档，了解有关许可证的相关函数，允许你添加其他许可证类型。所有这些函数的参数名称允许你指定软件的版权持有者—如果你为公司或机构工作，法律版权可能属于他们，值得检查这一点。

usethis::use_readme_rmd() 函数会添加一个空的 .Rmd 文件，你可以在其中添加代码和文本，该文件将被构建成一个常规的 markdown 文件，并作为 README 文件显示在你仓库的 GitHub 前端页面上。最少在此文件中添加描述你的包的目标、基本用法和安装说明。

向文档中添加一个有用的信息是指明开发的阶段。usethis::use_lifecycle_badge()函数可以让你创建一个漂亮的小图形徽章，显示你的包的开发进度。你可以作为第一个参数使用的术语可以在这里查看：www.tidyverse.org/lifecycle/。与此相关的是usethis::use_version()函数，它将帮助你递增软件的主要、次要或修复版本。

在步骤 2中，我们管理包所需的依赖项。当用户安装你的包时，包管理软件应会自动安装这些依赖项；R 要求它们放置在包元数据描述文件中的特定位置。usethis::use_package()函数会为你完成这项工作。

在步骤 3中，我们使用usethis::use_git()函数在当前目录创建一个本地的git仓库；它还会对当前代码进行初始提交。usethis::browse_github_token()函数会打开一个浏览器窗口，并导航到 GitHub 页面，允许你获取 GitHub 访问令牌，以便你的 R 会话能够与 GitHub 交互。获得令牌后，usethis::use_github()将把本地的git仓库上传至 GitHub，设置远程仓库并推送代码。你只需要执行一次。当git和 GitHub 仓库存在时，你需要手动管理版本，使用 RStudio 的git面板或git的命令行版本。

在步骤 4中，我们看到远程用户如何安装你的包，只需要使用devtools::install_github()函数，并提供适当的用户名和仓库名。

构建单元测试套件以确保函数按预期工作

大多数程序员会过度测试代码，单元测试的实践应运而生，让我们有了一种可以自动化并帮助减少构建中等复杂代码项目所需时间的正式测试方法。一个设计良好并维护的软件包，会为尽可能多的组件函数提供单元测试套件。在本食谱中，我们将学习如何使用usethis包添加组件文件和文件夹，以便创建一个使用testthat包的自动化测试套件。本书的范围不包括详细探讨为什么以及如何编写测试的哲学，但你可以查看testthat包的文档，网址是：testthat.r-lib.org/，作为一个很好的入门指南。

准备工作

我们需要usethis、testthat和devtools包。

如何操作...

使用以下步骤构建单元测试套件，以确保函数按预期工作：

1. 创建测试文件夹结构：

usethis::use_testthat()

1. 添加新的测试：

usethis::use_test("adds")

test_that("addition works", {
 expect_equal(1 + 1, 2)
})

1. 运行实际的测试：

devtools::test()

它是如何工作的...

第 1 步 是一个典型的 usethis 风格函数，创建一些你的包所需的常见文件系统组件——use_testthat() 只是构建了 testthat 测试引擎所需要的文件夹结构。

第 2 步 中，usethis::use_test() 函数开始工作，创建一个测试文件——它使用函数参数的值作为文件名的后缀，因此在此情况下，使用 adds 作为参数，会在 tests/testthat 文件夹中创建一个名为 test-adds.R 的文件。然后，我们可以将 tests 添加到该文件中。每个测试将遵循此步骤第二行所示的基本模式。调用 test_that() 函数；其第一个参数是打印到控制台的文本，显示当前正在进行的测试。第二个参数是来自 testthat 包的断言块，比较函数的输出与预期值。如果两者匹配，则测试通过；否则，测试失败。testthat 中有许多断言，可以测试多种类型的输出和对象。你可以在文档中查看：testthat.r-lib.org/。请注意，测试应该在测试文件中完成并保存，而不是直接在控制台中输入。

在 第 3 步 中，我们在控制台使用 devtools::test() 函数自动运行测试套件。测试结果会打印到控制台，你可以根据需要修改函数，然后重新运行此步骤。

使用 Travis 的持续集成来保持代码经过测试并保持最新

持续集成（CI）是一种团队编程实践，旨在帮助大团队在同一个项目中保持代码、依赖关系和测试的最佳协作。为此开发的工具同样可以帮助我们管理自己的软件项目，解决由于自己更新、所依赖包的更新，甚至在某些情况下，R 和操作系统更新引发的问题。Travis.CI 是 devtools 包支持的一个 CI 服务。通过将 Travis.CI 集成到你的项目中，Travis 服务器将构建一个新的虚拟计算机，安装操作系统，安装 R 及你的包所需的所有依赖包，然后安装你的包并运行测试套件。Travis 会将结果发送给你。这个过程会周期性地重复，尤其是每次你向 GitHub 推送代码时，以便你及时了解哪些代码出现了问题，并能及早解决问题。在这个食谱中，我们将介绍如何为你的包项目设置 Travis。

准备工作

对于这个食谱，你需要 usethis 包和一个托管在 GitHub 上的包项目。如果你还没有设置这些，前面章节的食谱将帮助你完成这些步骤。

如何实现...

为了使用 Travis 的持续集成来保持代码经过测试并保持最新，我们需要创建一个 .travis.yml 文件：

usethis::use_travis()

它是如何工作的...

这段代码中的唯一一行会在你的包目录的根目录中创建一个名为.travis.yml的文件。这个文件在 GitHub 上作为钩子使用，因此，一旦仓库被更新，Travis.CI服务器将会执行新的虚拟服务器构建和打包，并运行测试，然后将结果通过电子邮件发送到与你的 GitHub 账户关联的邮箱地址。尽管这只有一行代码，但这可能是本书中最具影响力的单行代码！.travis.yml文件包含了 Travis 构建的配置选项，并且可以添加很多内容来自定义输出。一个常见的添加内容如下：

warnings_are_errors: false

这将告诉 Travis，R 代码中的警告不应视为错误，也不会导致构建失败。

构建可能需要一些时间；即便是简单的代码也可能需要 15 分钟。更复杂的项目将需要更长的时间。