Analysis of Lloyd et al (2020) for manuscript case study

ref:lloyd2020evid differential abundance

Perform calculations for the Lloyd et al (2020) case study.

Michael R. McLaren true
2022-01-08

Setup

Load libraries,

library(here)
library(tidyverse)
# library(tidyxl)
library(unpivotr)

Our source data is Table 1 from Lloyd et al. (2020), which I’ve manually input into a CSV file. This table has the form of a ‘pivot table’ with nested headers; we can properly read it in and create a ‘tidy’ version table using the unpivotr package, following this example. First we’ll read it in with one row per spreadsheet cell,

all_cells <- here('notebook/_data/lloyd2020evid/table1.csv') %>%
  read_csv(col_names = FALSE) %>%
  as_cells %>%
  # filter(!is.na(hr), chr != 'Clade') %>%
  glimpse
Rows: 56
Columns: 4
$ row       <int> 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1,…
$ col       <int> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3,…
$ data_type <chr> "chr", "chr", "chr", "chr", "chr", "chr", "chr", "…
$ chr       <chr> NA, NA, "Clade", "Bathyarchaeota, or MCG", "Group …

Confirm that this data matches Table 1,

all_cells %>% rectify %>% knitr::kable()
row/col 1(A) 2(B) 3(C) 4(D) 5(E) 6(F) 7(G)
1 NA WOR marine sediment NA NA NA CLB marine sediment NA
2 NA Core 30 NA Core 32 NA Incubation 2 Incubation 3
3 Clade FRAxC qPCR FRAxC qPCR FRAxC FRAxC
4 Bathyarchaeota, or MCG 10.1 10.3 8.9 10 6.4 6.3
5 Group C3, MCG-related 16.2 NA 14.9 NA 2.1 1.8
6 Thermoprofundales (MBG-D) 4.6 NA 3.4 4.8 3.4 4.4
7 Thermoprofundales (MG-III) 3.7 NA 2.8 NA 1.4 1.5
8 Thermoprofundales (20a-9) 7.1 NA 3.7 NA 0.8 0.7

Note, the values of the table are doubling times.

Now, munge the data into tidy form, with one row per measurement. Additionally, compute the doubling rate as the reciprocal of the doubling time.

# abbreviations: exp for experiment, meas for measurement
x <- all_cells %>%
  behead("up-left", exp_type) %>%
  behead("up-left", exp_name) %>%
  behead("up", meas_type) %>%
  behead("left", clade) %>%
  select(clade, starts_with('exp'), meas_type, doubling_time = chr) %>%
  mutate(
    across(doubling_time, as.numeric),
    doubling_rate = 1 / doubling_time,
  ) %>%
  arrange(clade, exp_type, exp_name, meas_type) %>%
  glimpse
Rows: 30
Columns: 6
$ clade         <chr> "Bathyarchaeota, or MCG", "Bathyarchaeota, or …
$ exp_type      <chr> "CLB marine sediment", "CLB marine sediment", …
$ exp_name      <chr> "Incubation 2", "Incubation 3", "Core 30", "Co…
$ meas_type     <chr> "FRAxC", "FRAxC", "FRAxC", "qPCR", "FRAxC", "q…
$ doubling_time <dbl> 6.4, 6.3, 10.1, 10.3, 8.9, 10.0, 2.1, 1.8, 16.…
$ doubling_rate <dbl> 0.15625000, 0.15873016, 0.09900990, 0.09708738…
dir.create('_output')
saveRDS(x, '_output/table1-tidy.rds')

Analysis

Let’s compute the error in doubling rate, assuming that the qPCR rate is accurate.

y <- x %>%
  select(-doubling_time) %>%
  pivot_wider(names_from = meas_type, values_from = doubling_rate) %>%
  mutate(error = FRAxC - qPCR)
y %>% 
  filter(!is.na(qPCR)) %>%
  select(-exp_type) %>% 
  knitr::kable(digits = 4)
clade exp_name FRAxC qPCR error
Bathyarchaeota, or MCG Core 30 0.0990 0.0971 0.0019
Bathyarchaeota, or MCG Core 32 0.1124 0.1000 0.0124
Thermoprofundales (MBG-D) Core 32 0.2941 0.2083 0.0858

Core 30

y %>% 
  # filter(str_detect(exp_type, 'WOR')) %>%
  filter(exp_name == 'Core 30') %>%
  select(-exp_type, -exp_name) %>% 
  knitr::kable(digits = 4)
clade FRAxC qPCR error
Bathyarchaeota, or MCG 0.0990 0.0971 0.0019
Group C3, MCG-related 0.0617 NA NA
Thermoprofundales (20a-9) 0.1408 NA NA
Thermoprofundales (MBG-D) 0.2174 NA NA
Thermoprofundales (MG-III) 0.2703 NA NA

Core 32

y %>% 
  filter(exp_name == 'Core 32') %>%
  select(-exp_type, -exp_name) %>% 
  knitr::kable(digits = 4)
clade FRAxC qPCR error
Bathyarchaeota, or MCG 0.1124 0.1000 0.0124
Group C3, MCG-related 0.0671 NA NA
Thermoprofundales (20a-9) 0.2703 NA NA
Thermoprofundales (MBG-D) 0.2941 0.2083 0.0858
Thermoprofundales (MG-III) 0.3571 NA NA

Compute the calculations for the manuscript text,

bathy_core30_fraxc <- x %>%
  filter(exp_name == 'Core 30', str_detect(clade, 'Bathyarchaeota'), 
    meas_type == 'FRAxC') %>%
  pull(doubling_rate)
bathy_core30_qpcr <- x %>%
  filter(exp_name == 'Core 30', str_detect(clade, 'Bathyarchaeota'), 
    meas_type == 'qPCR') %>%
  pull(doubling_rate)
bathy_core32_fraxc <- x %>%
  filter(exp_name == 'Core 32', str_detect(clade, 'Bathyarchaeota'), 
    meas_type == 'FRAxC') %>%
  pull(doubling_rate)
bathy_core32_qpcr <- x %>%
  filter(exp_name == 'Core 32', str_detect(clade, 'Bathyarchaeota'), 
    meas_type == 'qPCR') %>%
  pull(doubling_rate)
mbgd_core32_fraxc <- x %>%
  filter(exp_name == 'Core 32', str_detect(clade, 'MBG-D'), 
    meas_type == 'FRAxC') %>%
  pull(doubling_rate)
mbgd_core32_qpcr <- x %>%
  filter(exp_name == 'Core 32', str_detect(clade, 'MBG-D'), 
    meas_type == 'qPCR') %>%
  pull(doubling_rate)

Draft manuscript text:

Our aim is to compare the growth rates for taxa with both FRAxC and qPCR values for the two sediment cores, as reported in Table 1. The first soil core included qPCR measurements of a single archaeal clade, Bathyarchaeota, for which growth rates by qPCR and FRAxC were nearly identical (doubling rates of 0.099/yr by FRAxC and 0.097/yr by qPCR). The second soil core included qPCR measurements of Bathyarchaeota and a second clade, Thermoprofundales/MBG-D. In this core, FRAxC and qPCR growth rates differed more substantially, with growth rates from FRAxC being larger by 0.012/yr for Bathyarchaeota (0.112/yr by FRAxC and 0.1/yr by qPCR) and by 0.086/yr for Thermoprofundales/MBG-D (0.294/yr by FRAxC and 0.208/yr by qPCR). A low number of experimental samples and noise in both the FRAxC and qPCR measurements place significant uncertainty in these measurements; however, the fact that FRAxC-derived growth rates are larger than qPCR-derived rates in all three cases is consistent with our hypothesis that mean efficiency decreases with depth in a manner that systematically biases FRAxC-derived rates to higher values. The differences in growth rate are small in absolute terms; however, the maximum observed difference of 0.086/yr suggests an error large enough to impact results for some taxa classified as positive growers by Lloyd et al. (2020), for which growth rates ranged between 0.04/yr and 0.5/yr. Overall, the comparison between FRAxC and qPCR measurements gives support to the study conclusions, but suggest that species at the lower end of this range of positive FRAxC-derived rates may in fact be merely persisting or even slowly declining in abundance.

Session info

Click for session info
sessioninfo::session_info()
─ Session info ───────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.1.2 (2021-11-01)
 os       Arch Linux
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/New_York
 date     2022-01-10
 pandoc   2.14.1 @ /usr/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────────
 package     * version date (UTC) lib source
 assertthat    0.2.1   2019-03-21 [1] CRAN (R 4.0.0)
 backports     1.4.1   2021-12-13 [1] CRAN (R 4.1.2)
 bit           4.0.4   2020-08-04 [1] CRAN (R 4.0.2)
 bit64         4.0.5   2020-08-30 [1] CRAN (R 4.0.2)
 bookdown      0.24    2021-09-02 [1] CRAN (R 4.1.1)
 broom         0.7.10  2021-10-31 [1] CRAN (R 4.1.2)
 bslib         0.3.1   2021-10-06 [1] CRAN (R 4.1.1)
 cachem        1.0.6   2021-08-19 [1] CRAN (R 4.1.1)
 cellranger    1.1.0   2016-07-27 [1] CRAN (R 4.0.0)
 cli           3.1.0   2021-10-27 [1] CRAN (R 4.1.1)
 colorspace    2.0-2   2021-08-11 [1] R-Forge (R 4.1.1)
 crayon        1.4.2   2021-10-29 [1] CRAN (R 4.1.1)
 DBI           1.1.1   2021-01-15 [1] CRAN (R 4.0.4)
 dbplyr        2.1.1   2021-04-06 [1] CRAN (R 4.0.5)
 digest        0.6.29  2021-12-01 [1] CRAN (R 4.1.2)
 distill       1.3     2021-10-13 [1] CRAN (R 4.1.1)
 downlit       0.4.0   2021-10-29 [1] CRAN (R 4.1.2)
 dplyr       * 1.0.7   2021-06-18 [1] CRAN (R 4.1.0)
 ellipsis      0.3.2   2021-04-29 [1] CRAN (R 4.1.0)
 evaluate      0.14    2019-05-28 [1] CRAN (R 4.0.0)
 fansi         0.5.0   2021-05-25 [1] CRAN (R 4.1.0)
 fastmap       1.1.0   2021-01-25 [1] CRAN (R 4.0.4)
 forcats     * 0.5.1   2021-01-27 [1] CRAN (R 4.0.4)
 fs            1.5.2   2021-12-08 [1] CRAN (R 4.1.2)
 generics      0.1.1   2021-10-25 [1] CRAN (R 4.1.1)
 ggplot2     * 3.3.5   2021-06-25 [1] CRAN (R 4.1.0)
 glue          1.5.1   2021-11-30 [1] CRAN (R 4.1.2)
 gtable        0.3.0   2019-03-25 [1] CRAN (R 4.0.0)
 haven         2.4.3   2021-08-04 [1] CRAN (R 4.1.1)
 here        * 1.0.1   2020-12-13 [1] CRAN (R 4.0.5)
 highr         0.9     2021-04-16 [1] CRAN (R 4.1.0)
 hms           1.1.1   2021-09-26 [1] CRAN (R 4.1.1)
 htmltools     0.5.2   2021-08-25 [1] CRAN (R 4.1.1)
 httr          1.4.2   2020-07-20 [1] CRAN (R 4.0.2)
 jquerylib     0.1.4   2021-04-26 [1] CRAN (R 4.1.0)
 jsonlite      1.7.2   2020-12-09 [1] CRAN (R 4.0.3)
 knitr         1.36    2021-09-29 [1] CRAN (R 4.1.1)
 lifecycle     1.0.1   2021-09-24 [1] CRAN (R 4.1.1)
 lubridate     1.8.0   2021-10-07 [1] CRAN (R 4.1.1)
 magrittr      2.0.1   2020-11-17 [1] CRAN (R 4.0.3)
 memoise       2.0.1   2021-11-26 [1] CRAN (R 4.1.2)
 modelr        0.1.8   2020-05-19 [1] CRAN (R 4.0.0)
 munsell       0.5.0   2018-06-12 [1] CRAN (R 4.0.0)
 nvimcom     * 0.9-102 2021-11-12 [1] local
 pillar        1.6.4   2021-10-18 [1] CRAN (R 4.1.1)
 pkgconfig     2.0.3   2019-09-22 [1] CRAN (R 4.0.0)
 purrr       * 0.3.4   2020-04-17 [1] CRAN (R 4.0.0)
 R6            2.5.1   2021-08-19 [1] CRAN (R 4.1.1)
 Rcpp          1.0.7   2021-07-07 [1] CRAN (R 4.1.0)
 readr       * 2.1.1   2021-11-30 [1] CRAN (R 4.1.2)
 readxl        1.3.1   2019-03-13 [1] CRAN (R 4.0.0)
 reprex        2.0.1   2021-08-05 [1] CRAN (R 4.1.1)
 rlang         0.4.12  2021-10-18 [1] CRAN (R 4.1.1)
 rmarkdown   * 2.11    2021-09-14 [1] CRAN (R 4.1.1)
 rprojroot     2.0.2   2020-11-15 [1] CRAN (R 4.0.3)
 rstudioapi    0.13    2020-11-12 [1] CRAN (R 4.0.3)
 rvest         1.0.2   2021-10-16 [1] CRAN (R 4.1.1)
 sass          0.4.0   2021-05-12 [1] CRAN (R 4.1.0)
 scales        1.1.1   2020-05-11 [1] CRAN (R 4.0.0)
 sessioninfo   1.2.2   2021-12-06 [1] CRAN (R 4.1.2)
 stringi       1.7.6   2021-11-29 [1] CRAN (R 4.1.2)
 stringr     * 1.4.0   2019-02-10 [1] CRAN (R 4.0.0)
 tibble      * 3.1.6   2021-11-07 [1] CRAN (R 4.1.2)
 tidyr       * 1.1.4   2021-09-27 [1] CRAN (R 4.1.1)
 tidyselect    1.1.1   2021-04-30 [1] CRAN (R 4.1.0)
 tidyverse   * 1.3.1   2021-04-15 [1] CRAN (R 4.1.0)
 tzdb          0.2.0   2021-10-27 [1] CRAN (R 4.1.2)
 unpivotr    * 0.6.2   2021-08-22 [1] CRAN (R 4.1.2)
 utf8          1.2.2   2021-07-24 [1] CRAN (R 4.1.0)
 vctrs         0.3.8   2021-04-29 [1] CRAN (R 4.1.0)
 vroom         1.5.7   2021-11-30 [1] CRAN (R 4.1.2)
 withr         2.4.3   2021-11-30 [1] CRAN (R 4.1.2)
 xfun          0.28    2021-11-04 [1] CRAN (R 4.1.2)
 xml2          1.3.3   2021-11-30 [1] CRAN (R 4.1.2)
 yaml          2.2.1   2020-02-01 [1] CRAN (R 4.0.0)

 [1] /home/michael/.local/lib/R/library
 [2] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────────
Lloyd, Karen G., Jordan T. Bird, Joy Buongiorno, Emily Deas, Richard Kevorkian, Talor Noordhoek, Jacob Rosalsky, and Taylor Roy. 2020. Evidence for a Growth Zone for Deep-Subsurface Microbial Clades in Near-Surface Anoxic Sediments.” Appl. Environ. Microbiol. 86 (19): 1–15. https://doi.org/10.1128/AEM.00877-20.

References