R使用phyloseq对象计算最丰富的分类群

Question

我想知道我计算任何分类单元相对丰度平均值的方法是否正确!!!

如果我想知道计算 phyloseq 对象 (GlobalPattern) 中每个科(或任何分类单元)的相对丰度(百分比)是否正确，如下所示:

    data("GlobalPatterns")

T <- GlobalPatterns %>% 
    tax_glom(., "Family") %>% 
    transform_sample_counts(function(x)100* x / sum(x)) %>% psmelt() %>% 
    arrange(OTU) %>% rename(OTUsID = OTU) %>% 
    select(OTUsID, Family, Sample, Abundance) %>%
    spread(Sample, Abundance)

T$Mean <- rowMeans(T[, c(3:ncol(T))])

FAM <- T[, c("Family", "Mean" ) ]

#order data frame  
FAM <- FAM[order(dplyr::desc(FAM$Mean)),]
rownames(FAM) <- NULL

head(FAM)
          Family          Mean
1    Bacteroidaceae     7.490944
2   Ruminococcaceae     6.038956
3   Lachnospiraceae     5.758200
4 Flavobacteriaceae     5.016402
5  Desulfobulbaceae     3.341026
6            ACK-M1     3.242808

在这种情况下，拟杆菌科是 GlobalPattern 所有样本(26 个样本和 19216 个 OTU)中最丰富的科，它在 26 个样本中的平均含量为 7.49% !!!! 使 T$Mean <- rowMeans(T[, c(3:ncol(T))]) 计算任何给定分类单元的平均值是否正确？

Answer 1

如果将所有样本汇集在一起​​，拟杆菌科的丰度最高。 然而，它仅在 2 个样本中具有最高丰度。 然而，没有其他分类单元在平均样本中具有更高的丰度。

让我们对所有步骤使用 dplyr 动词，以获得更具描述性和一致性的代码:

library(tidyverse)
library(phyloseq)
#> Creating a generic function for 'nrow' from package 'base' in package 'biomformat'
#> Creating a generic function for 'ncol' from package 'base' in package 'biomformat'
#> Creating a generic function for 'rownames' from package 'base' in package 'biomformat'
#> Creating a generic function for 'colnames' from package 'base' in package 'biomformat'
data(GlobalPatterns)

data <-
  GlobalPatterns %>%
  tax_glom("Family") %>%
  transform_sample_counts(function(x)100* x / sum(x)) %>%
  psmelt() %>%
  as_tibble()

# highest abundance: all samples pooled together
data %>%
  group_by(Family) %>%
  summarise(Abundance = mean(Abundance)) %>%
  arrange(-Abundance)
#> # A tibble: 334 × 2
#>    Family             Abundance
#>    <chr>                  <dbl>
#>  1 Bacteroidaceae          7.49
#>  2 Ruminococcaceae         6.04
#>  3 Lachnospiraceae         5.76
#>  4 Flavobacteriaceae       5.02
#>  5 Desulfobulbaceae        3.34
#>  6 ACK-M1                  3.24
#>  7 Streptococcaceae        2.77
#>  8 Nostocaceae             2.62
#>  9 Enterobacteriaceae      2.55
#> 10 Spartobacteriaceae      2.45
#> # … with 324 more rows

# sanity check: is total abundance of each sample 100%?
data %>%
  group_by(Sample) %>%
  summarise(Abundance = sum(Abundance)) %>%
  pull(Abundance) %>%
  `==`(100) %>%
  all()
#> [1] TRUE

# get most abundant family for each sample individually
data %>%
  group_by(Sample) %>%
  arrange(-Abundance) %>%
  slice(1) %>%
  select(Family) %>%
  ungroup() %>%
  count(Family, name = "n_samples") %>%
  arrange(-n_samples)
#> Adding missing grouping variables: `Sample`
#> # A tibble: 18 × 2
#>    Family              n_samples
#>    <chr>                   <int>
#>  1 Desulfobulbaceae            3
#>  2 Bacteroidaceae              2
#>  3 Crenotrichaceae             2
#>  4 Flavobacteriaceae           2
#>  5 Lachnospiraceae             2
#>  6 Ruminococcaceae             2
#>  7 Streptococcaceae            2
#>  8 ACK-M1                      1
#>  9 Enterobacteriaceae          1
#> 10 Moraxellaceae               1
#> 11 Neisseriaceae               1
#> 12 Nostocaceae                 1
#> 13 Solibacteraceae             1
#> 14 Spartobacteriaceae          1
#> 15 Sphingomonadaceae           1
#> 16 Synechococcaceae            1
#> 17 Veillonellaceae             1
#> 18 Verrucomicrobiaceae         1

^{由 reprex package 创建于 2022-06-10 (v2.0.0)}