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

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

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

如果我想知道,在 phyloseq 对象(GlobalPattern)中计算每个家庭(或任何 Taxon)的相对丰度(百分比)将是正确的:

    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 & lt;-rowMeans(T [,c(3:ncol(T))])计算任何给定 Taxon 的平均值是正确的?

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

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

本站系公益性非盈利分享网址,本文来自用户投稿,不代表码文网立场,如若转载,请注明出处

(234)
BIM360是否具有用于FieldAPI的oAuth登录
上一篇
在一段时间内可能发生多少次 I/O中断
下一篇

相关推荐

发表评论

登录 后才能评论

评论列表(37条)