如果您已经知道一些R，则表格数据相对来说比较简单，但是分类数据是分层的，因此使用起来更加困难。当还有其他与分类法相关的数据时，它将变得更加复杂。
例如，如果我们要从数据中删除细菌门，是否还要删除该门内的所有分类单元？如果是这样，我们如何识别该门上的哪个分类单元？根据数据的格式，这可能很难。如果其他数据（例如OTU计数）与该门类中的分类单元相关联，我们是否也将其删除？毕竟，即使除去其中的门和分类单元，我们仍然知道所有相关数据都是细菌，所以我们是否将它们重新分配给细菌分类单元，还是将这些数据丢弃？这些问题的答案将取决于我们的数据以及我们要如何处理。
加载示例数据
如果您是从本节开始讲习班，或者在上一节中运行代码时遇到问题，请使用以下内容加载本节中使用的数据。您可以在此处下载“ parsed_data.Rdata”文件。如果您刚刚完成了上一节，则可以忽略并继续。
load("parsed_data.Rdata")
子集表格数据
taxa包提供的用于在taxmap对象中划分分类信息子集的函数(我们的示例数据采用的格式)是根据dplyr包中用于操作表格数据的函数建模的(Wickham和Francois(2015))。这些dplyr函数是标准R子集(例如[、[[和$)的替代品，它们被设计得更加直观和一致。由于理解这些函数将有助于处理分类数据的类似(和更复杂)函数，所以我们将从实践这些函数开始。幸运的是，我们有一个需要一些子集的表格数据集:我们的样本数据。
首先，让我们再看一下示例数据：
print(sample_data)
### A tibble: 1,698 x 16## SampleID Name Plant_ID Type Experiment Cohort Harvested Age Site Treatment Line Genotype
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 M1024P17… R_026… R_026 root fieldBCMA 2008 2011 3 LTM field 26 ril
## 2 M1024P17… R_073… R_073 root fieldBCMA 2008 2011 3 LTM field 73 ril
## 3 M1024P17… R_088… R_088 root fieldBCMA 2009 2011 2 LTM field 88 ril
## 4 M1024P18… R_156… R_156 root fieldBCMA 2009 2011 2 LTM field 156 ril
## 5 M1955P804 1_A_1… 1_A_1 root ecoGH NA 2011 NA Duke MAHsoil par1… PAR
## 6 M1956P837 1_A_1… 1_A_12 root ecoGH NA 2011 NA Duke MAHsoil mah3… MAH
## 7 M1957P983 1_A_3… 1_A_3 root ecoGH NA 2011 NA Duke MAHsoil silL… SIL
## 8 M1956P845 1_A_4… 1_A_4 root ecoGH NA 2011 NA Duke MAHsoil par1… PAR
## 9 M1957P987 1_A_7… 1_A_7 root ecoGH NA 2011 NA Duke MAHsoil par9… PAR
## 10 M1957P923 1_A_8… 1_A_8 root ecoGH NA 2011 NA Duke MAHsoil mil5… MIL
## # ... with 1,688 more rows, and 4 more variables: Block <chr>, oldPlate <chr>, newPlate <chr>,
## # Analysis <chr>
这是一个很大的数据集，并非Wagner等人的主要出版物中使用了全部1698个样本。（2016）。只有“生态型”实验中的样本才是主要研究的一部分。使用标准R子集，我们可以像这样删除它们：
sample_data <- sample_data[sample_data$Experiment == "ecotypes", ]
使用这些dplyr功能，相同的操作将是：
library(dplyr)
sample_data <- filter(sample_data, Experiment== "ecotypes")
注意我们如何Experiment在filter函数中单独使用该列。这是filter功能以及许多其他功能的特殊属性，dplyr称为非标准评估（NSE）。如果我们在函数之外尝试该操作，则会收到错误消息。例如，以下将不起作用：
is_ecotype <- Experiment== "ecotypes"
## Error in eval(expr, envir, enclos): object 'Experiment' not found
大多数dplyr功能都是以这种方式工作的，类似功能在taxa处理分类数据时也是如此。如果我们有多个过滤条件，我们可以简单地向函数添加更多参数。在下面的代码中，我们将表子集仅划分为三个站点之一且具有3年历史的植物的行。
sample_data <- filter(sample_data,Site %in% c("Mah", "Jam", "Sil"),Age == 3)
通常，这等效于&在常规子集中组合多个条件，但对于大型数据集则更快。NA与base R子集不同，它也不会返回条件求值的值。
到目前为止，我们使用filterR子集也可以轻松完成所有操作，但是dplyr引入了一些更高级的功能，可以节省时间/键入。我们不会详细介绍这些内容，但是很高兴知道它们的存在。一种是“分组”的概念。当按一列中的值对行进行“分组”时，在某些函数中，与该列中每个唯一值相对应的行将被视为一个单位。n()是一个这样的函数，它计算每个组中的项目数。我们可以通过计算每个“ Plant_ID”出现的次数来过滤掉没有根和叶样本的任何植物：
sample_data <- sample_data %>%
group_by(Plant_ID) %>%
filter(n() == 2)
剩下292个样本。注意使用%>%“管道”运算符。这将获取之前函数的结果，并将其用作之后函数的第一个参数。例如，seq(1:10) %>% length()与相同length(seq(1:10))。这允许以可读的方式将许多命令链接在一起。接下来，我们将把丰度数据子集化为这些样本，以将数据集缩减至标准个人计算机的可管理大小。
细分分类法
filter分类数据的类似物filter_taxa来自taxa包。这以相似的方式工作，但是对于分类关系以及使用与已删除的分类单元关联的数据进行了更多选择。让我们再看一下我们的数据对象：
print(obj)
## <Taxmap>
## 1558 taxa: aab. Unassigned, aac. Root ... chx. Methanospirillaceae, chy.
## 1558 edges: NA->aab, NA->aac, aac->aad, aac->aae ... bel->chw, ays->chx, bem->chy
## 1 data sets:
## otu_counts:
## # A tibble: 47,806 x 1,707
## taxon_id OTU_ID M1024P1833 M1551P81 M1551P57 M1551P85 M1551P28 M1551P29 M1551P38
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 aab 1 0. 0. 0. 0. 0. 0. 0.
## 2 ben 2 41. 22. 4. 726. 112492. 2. 413.
## 3 beo 3 67. 65. 12. 13514. 1. 4. 13314.
## # ... with 4.78e+04 more rows, and 1,698 more variables: M1551P90 <dbl>,
## # M1551P71 <dbl>, M1551P12 <dbl>, M1551P84 <dbl>, M1551P48 <dbl>, M1551P4 <dbl>,
## # M1551P52 <dbl>, M1551P3 <dbl>, M1551P15 <dbl>, M1551P31 <dbl>, …
## 0 functions:
如果我们看第二行，可以看到某些分类单元没有名称。这些是我们在上一节（例如，）中解析的分类中仅具有排名信息的分类。由于这些信息不会添加任何信息，因此请删除它们。o__"Root;k__Bacteria;p__Chlorobi;c__SJA-28;o__;f__"
library(taxa)
obj <- filter_taxa(obj, taxon_names != "")
print(obj)
## ## 959 taxa:aab. Unassigned, aac. Root ... chs. Alcanivoracaceae, chx. Methanospirillaceae
## 959 edges: NA->aab, NA->aac, aac->aad, aac->aae ... anq->cho, apt->chs, ays->chx## 1 data sets:
## otu_counts:
## # A tibble: 47,806 x 1,707## taxon_id OTU_ID M1024P1833 M1551P81 M1551P57 M1551P85 M1551P28 M1551P29 M1551P38
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>## 1 aab 1 0. 0. 0. 0. 0. 0. 0.
## 2 ben 2 41. 22. 4. 726. 112492. 2. 413.
## 3 beo 3 67. 65. 12. 13514. 1. 4. 13314.
## # ... with 4.78e+04 more rows, and 1,698 more variables: M1551P90 <dbl>,
## # M1551P71 <dbl>, M1551P12 <dbl>, M1551P84 <dbl>, M1551P48 <dbl>, M1551P4 <dbl>,
## # M1551P52 <dbl>, M1551P3 <dbl>, M1551P15 <dbl>, M1551P31 <dbl>, …
## 0 functions:
有与无名分类单元相关联的OTU，但是没有过滤掉OTU，那么它们发生了什么？默认情况下，它们被重新分配给通过过滤器的最接近的超分类单元。例如，对于分类 "Root;k__Bacteria;p__Chlorobi;c__SJA-28;o__;f__"，OTU将被重新分配给“ SJA-28”类。
像一样filter，filter_taxa可以使用输入对象中包含的变量，就好像它们是独立变量一样。在这种情况下，taxon_names实际上是这样引用的函数。
head(taxon_names(obj))
## aab aac aad aae aaf
## "Unassigned" "Root" "Bacteria" "Archaea" "Proteobacteria"
## aag
## "Actinobacteria"
可以以这种方式使用的其他值包括表中的列和存储在taxmap对象中obj$data列表中的列表/向量。由于taxmap对象可以存储映射到分类法的任意数量的列表、向量或表，因此可以通过名称引用许多变量。函数返回变量的名称和位置，可以这样使用:
head(all_names(obj),20)
## taxon_names taxon_ids taxon_indexes
## "taxon_names" "taxon_ids" "taxon_indexes"
## classifications n_supertaxa n_supertaxa_1
## "classifications" "n_supertaxa" "n_supertaxa_1"
## n_subtaxa n_subtaxa_1 n_leaves
## "n_subtaxa" "n_subtaxa_1" "n_leaves"
## n_leaves_1 taxon_ranks is_root
## "n_leaves_1" "taxon_ranks" "is_root"
## is_stem is_branch is_leaf
## "is_stem" "is_branch" "is_leaf"
## is_internode n_obs n_obs_1
## "is_internode" "n_obs" "n_obs_1"
## data$otu_counts$OTU_ID data$otu_counts$M1024P1833
## "OTU_ID" "M1024P1833"
length(all_names(obj))
## [1] 1724
接下来，让我们也过滤掉细菌中没有的任何东西，因为研究的重点是细菌。
obj <- filter_taxa(obj, taxon_names == "Bacteria", subtaxa = TRUE)
print(obj)
## ## 898 taxa:aad. Bacteria, aaf. Proteobacteria ... cho. wb1_P06, chs. Alcanivoracaceae
## 898 edges: NA->aad, aad->aaf, aad->aag, aad->aah ... bea->chi, anq->cho, apt->chs## 1 data sets:
## otu_counts:
## # A tibble: 33,641 x 1,707## taxon_id OTU_ID M1024P1833 M1551P81 M1551P57 M1551P85 M1551P28 M1551P29 M1551P38
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>## 1 ben 2 41. 22. 4. 726. 112492. 2. 413.
## 2 beo 3 67. 65. 12. 13514. 1. 4. 13314.
## 3 ben 4 3229. 13679. 1832. 951. 113. 2496. 567.
## # ... with 3.364e+04 more rows, and 1,698 more variables: M1551P90 <dbl>,
## # M1551P71 <dbl>, M1551P12 <dbl>, M1551P84 <dbl>, M1551P48 <dbl>, M1551P4 <dbl>,
## # M1551P52 <dbl>, M1551P3 <dbl>, M1551P15 <dbl>, M1551P31 <dbl>, …
## 0 functions:
注意使用subtaxa = TRUE; 这意味着通过过滤器的所有分类单元的子分类单元也应保留。如果不选择此选项，则只有一个分类单元，所有OTU都分配给它：
filter_taxa(obj, taxon_names == "Bacteria")
## <Taxmap>
## 1 taxa: aad. Bacteria
## 1 edges: NA->aad
## 1 data sets:
## otu_counts:
## # A tibble: 33,641 x 1,707
## taxon_id OTU_ID M1024P1833 M1551P81 M1551P57 M1551P85 M1551P28 M1551P29 M1551P38
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 aad 2 41. 22. 4. 726. 112492. 2. 413.
## 2 aad 3 67. 65. 12. 13514. 1. 4. 13314.
## 3 aad 4 3229. 13679. 1832. 951. 113. 2496. 567.
## # ... with 3.364e+04 more rows, and 1,698 more variables: M1551P90 <dbl>,
## # M1551P71 <dbl>, M1551P12 <dbl>, M1551P84 <dbl>, M1551P48 <dbl>, M1551P4 <dbl>,
## # M1551P52 <dbl>, M1551P3 <dbl>, M1551P15 <dbl>, M1551P31 <dbl>, …
## 0 functions:
接下来，让我们从丰度矩阵中删除在本节开始时从样本数据中删除的列。
obj$data$otu_counts <- obj$data$otu_counts[c("taxon_id", "OTU_ID", sample_data$SampleID)]
这使我们获得了292个样本中的33641个OTU的数据集，按898类分类：
print(obj)
## <Taxmap>
## 898 taxa: aad. Bacteria, aaf. Proteobacteria ... cho. wb1_P06, chs. Alcanivoracaceae
## 898 edges: NA->aad, aad->aaf, aad->aag, aad->aah ... bea->chi, anq->cho, apt->chs
## 1 data sets:
## otu_counts:
## # A tibble: 33,641 x 294
## taxon_id OTU_ID M1981P563 M1977P1709 M1980P502 M1981P606 M1981P599 M1981P611
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ben 2 7. 0. 2. 7. 1. 7.
## 2 beo 3 10251. 59. 15857. 309. 134. 90.
## 3 ben 4 1447. 2496. 2421. 3331. 143. 1847.
## # ... with 3.364e+04 more rows, and 286 more variables: M1981P587 <dbl>,
## # M1982P729 <dbl>, M1979P447 <dbl>, M1982P701 <dbl>, M1982P661 <dbl>,
## # M1978P293 <dbl>, M1980P462 <dbl>, M1555P221 <dbl>, M1980P471 <dbl>,
## # M1979P422 <dbl>, …
## 0 functions:
子分类数据分类
在上一节中，我们taxmap使用来对对象中的分类法进行子集化filter_taxa。默认情况下，删除分类单元时，分配给它们的OTU将重新分配给超级分类单元或也将其删除。我们还可以使用该函数直接子集OTU表并删除不再需要的任何分类单元filter_obs。函数名称中的“ obs”代表“observations”。与filter您以前用来对表进行子集化的不同，它处理对象中filter_obs的数据集taxmap，例如我们的OTU表。
在上一节中，我们删除了Wagner等人中未出现的样本（即OTU表中的列）。（2016）出版物。这可能意味着某些OTU仅存在于那些样本中，而在其余样本中现在没有读取。让我们删除那些以使数据集更易于管理。首先，我们需要找到没有读取的OTU（即行）。我们可以通过使用函数rowSums来获取每个OTU的所有样本中的读取次数来做到这一点。
has_no_reads <- rowSums(obj$data$otu_counts[, sample_data$SampleID]) == 0
由于该表otu_counts还具有分类单元ID和OTU ID列，因此我们需要在使用之前将这些列子集化为仅样本rowSums。上面的命令生成一个logical向量，该向量是可用于子集化事物的一系列TRUE/ FALSE值。我们可以sum用来计算TRUE值的数量，从而计算没有读取的OTU的数量。
sum(has_no_reads)
## [1] 11957
现在有33641个OTU中的11957个没有读取。现在我们知道了它们是什么，我们可以使用删除它们filter_obs。
filter_obs(obj, "otu_counts", ! has_no_reads) # note the ! negation operator
## ## 898 taxa:aad. Bacteria, aaf. Proteobacteria ... cho. wb1_P06, chs. Alcanivoracaceae
## 898 edges: NA->aad, aad->aaf, aad->aag, aad->aah ... bea->chi, anq->cho, apt->chs## 1 data sets:
## otu_counts:
## # A tibble: 21,684 x 294## taxon_id OTU_ID M1981P563 M1977P1709 M1980P502 M1981P606 M1981P599 M1981P611
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>## 1 ben 2 7. 0. 2. 7. 1. 7.
## 2 beo 3 10251. 59. 15857. 309. 134. 90.
## 3 ben 4 1447. 2496. 2421. 3331. 143. 1847.
## # ... with 2.168e+04 more rows, and 286 more variables: M1981P587 <dbl>,
## # M1982P729 <dbl>, M1979P447 <dbl>, M1982P701 <dbl>, M1982P661 <dbl>,
## # M1978P293 <dbl>, M1980P462 <dbl>, M1555P221 <dbl>, M1980P471 <dbl>,
## # M1979P422 <dbl>, …
## 0 functions:
请注意，这如何删除“ otu_counts”表中的行，但没有删除任何分类单元。要删除在任何OTU中不再观察到的分类单元，我们需要添加以下drop_taxa选项：
obj <- filter_obs(obj, "otu_counts", ! has_no_reads, drop_taxa = TRUE)
如果obj$data这些分类单元的数据中还有其他表格，那么除非drop_obs = FALSE添加，否则也将被删除。由于只有一个表，因此我们不必为此担心。
网友评论