Quantifying isotope incorporation

Introduction

There are two prevaling methods for using HTS-SIP data to estimate the amount of isotope that each OTU incorporated:

  • q-SIP
  • delta_BD (a complementary analysis to HR-SIP)
    • Note: delta_BD is formally written as: $\Delta\hat{BD}$

In this vignette, we are going to show how to run both analyses and also compare the results a bit.

Dataset

First, let’s load some packages including HTSSIP.

library(dplyr)
library(ggplot2)
library(HTSSIP)

OK. We’re going to be using 2 data files:

  • HTS-SIP data (phyloseq format)
  • qPCR data (total 16S rRNA gene copies per gradient fraction)

We’ll be using the dataset that we simulated in the HTSSIP_sim vignette.

The phyloseq object is similar to the dataset in the other vignettes.

# HTS-SIP data
physeq_rep3
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 6 taxa and 144 samples ]
## sample_data() Sample Data:       [ 144 samples by 5 sample variables ]

The associated qPCR data is a list of length = 2.

# qPCR data (list object)
physeq_rep3_qPCR %>% names
## [1] "raw"     "summary"

For the analyses in this vignette, we only need the ‘summary’ table.

# qPCR data (list object)
physeq_rep3_qPCR_sum = physeq_rep3_qPCR$summary
physeq_rep3_qPCR_sum %>% head(n=4)
##   IS_CONTROL                  Sample Buoyant_density qPCR_tech_rep_mean
## 1      FALSE 13C-Glu_rep1_1.671712_2        1.671712           46598172
## 2      FALSE 13C-Glu_rep1_1.671722_1        1.671722           42299449
## 3      FALSE 13C-Glu_rep1_1.680311_3        1.680311           53536519
## 4      FALSE 13C-Glu_rep1_1.683540_4        1.683540           51449609
##   qPCR_tech_rep_sd     Gradient Fraction Treatment Replicate
## 1          9055833 13C-Glu_rep1        2   13C-Glu         1
## 2         16369298 13C-Glu_rep1        1   13C-Glu         1
## 3         18265667 13C-Glu_rep1        3   13C-Glu         1
## 4          3796576 13C-Glu_rep1        4   13C-Glu         1

q-SIP

OK. Let’s quantify isotope incorporation witht the q-SIP method.

# transforming OTU counts
physeq_rep3_t = OTU_qPCR_trans(physeq_rep3, physeq_rep3_qPCR_sum)

# calculating atom fraction excess
atomX = qSIP_atom_excess(physeq_rep3_t,
                         control_expr='Treatment=="12C-Con"',
                         treatment_rep='Replicate')
atomX %>% names
## [1] "W" "A"

The resulting list object contains 2 data.frames. We are interested in the ‘A’ table, which contains estimated BD shifts (Z) and atom fraction excess (A).

atomX$A %>% head(n=4)
## # A tibble: 4 × 9
##   OTU    Wlab Wlight      Z    Gi Mlight Mheavymax  Mlab     A
##   <chr> <dbl>  <dbl>  <dbl> <dbl>  <dbl>     <dbl> <dbl> <dbl>
## 1 OTU.1  1.72   1.70 0.0222 0.636   308.      318.  312. 0.412
## 2 OTU.2  1.71   1.69 0.0117 0.586   308.      318.  310. 0.218
## 3 OTU.3  1.73   1.70 0.0347 0.608   308.      318.  314. 0.644
## 4 OTU.4  1.73   1.70 0.0261 0.705   308.      318.  313. 0.484

Next, let’s calculate bootstrap confidence intervales for the atom fraction excess estimations.

# calculating bootstrapped CI values
df_atomX_boot = qSIP_bootstrap(atomX, n_boot=100)
## Warning: `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
## `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(ret)`.
df_atomX_boot %>% head(n=4)
## # A tibble: 4 × 11
##   OTU    Wlab Wlight      Z    Gi Mlight Mheavymax  Mlab     A A_CI_low
##   <chr> <dbl>  <dbl>  <dbl> <dbl>  <dbl>     <dbl> <dbl> <dbl>    <dbl>
## 1 OTU.1  1.72   1.70 0.0222 0.636   308.      318.  312. 0.412  -0.0401
## 2 OTU.2  1.71   1.69 0.0117 0.586   308.      318.  310. 0.218   0.121 
## 3 OTU.3  1.73   1.70 0.0347 0.608   308.      318.  314. 0.644   0.431 
## 4 OTU.4  1.73   1.70 0.0261 0.705   308.      318.  313. 0.484   0.337 
## # ℹ 1 more variable: A_CI_high <dbl>

delta_BD

Now for delta_BD. The setup is easier because we are not using qPCR data, just relative abundances from 16S rRNA sequence data.

df_dBD = delta_BD(physeq_rep3, control_expr='Treatment=="12C-Con"')
## Warning: Supplying `...` without names was deprecated in tidyr 1.0.0.
## ℹ Please specify a name for each selection.
## ℹ Did you want `data = nest_cols`?
## ℹ The deprecated feature was likely used in the HTSSIP package.
##   Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: `cols` is now required when using `unnest()`.
## ℹ Please use `cols = c(data)`.
df_dBD %>% head(n=4)
## # A tibble: 4 × 4
##   OTU   CM_control CM_treatment delta_BD
##   <chr>      <dbl>        <dbl>    <dbl>
## 1 OTU.1       1.69         1.72   0.0260
## 2 OTU.2       1.69         1.70   0.0130
## 3 OTU.3       1.69         1.74   0.0562
## 4 OTU.4       1.71         1.73   0.0208

Comparing results

Let’s plot the data and compare all of the results. First, let’s join all of the data into one table for plotting. We’ll also format it for plotting.

# checking & joining data 
stopifnot(nrow(df_atomX_boot) == nrow(df_dBD))
df_j = dplyr::inner_join(df_atomX_boot, df_dBD, c('OTU'='OTU'))
stopifnot(nrow(df_atomX_boot) == nrow(df_j))

# formatting data for plotting
df_j = df_j %>%
  dplyr::mutate(OTU = reorder(OTU, -delta_BD))

OK. Time to plot!

# plotting BD shift (Z)
ggplot(df_j, aes(OTU)) +
  geom_point(aes(y=Z), color='blue') +
  geom_point(aes(y=delta_BD), color='red') +
  geom_hline(yintercept=0, linetype='dashed', alpha=0.5) +
  labs(x='OTU', y='BD shift (Z)') +
  theme_bw() +
  theme(
    axis.text.x = element_blank()
  )

In the figure, red points are delta_BD and blue points are q-SIP. It’s easy to see that delta_BD is a lot more variable than q-SIP. This is likely due to a high influence of compositional data artifacts on delta_BD versus q-SIP.

Let’s make a boxplot to show the difference in estimation variance between the two methods.

# plotting BD shift (Z): boxplots

## formatting the table
df_j_g = df_j %>%
  dplyr::select(OTU, Z, delta_BD) %>%
  tidyr::gather(Method, BD_shift, Z, delta_BD) %>%
  mutate(Method = ifelse(Method == 'Z', 'q-SIP', 'delta-BD'))

## plotting 
ggplot(df_j_g, aes(Method, BD_shift)) +
  geom_boxplot() +
  geom_hline(yintercept=0, linetype='dashed', alpha=0.5) +
  labs(x='Method', y='BD shift (Z)') +
  theme_bw()

The boxplot helps to summarize how much more variance delta_BD produces versus q-SIP.

Session info

sessionInfo()
## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] phyloseq_1.51.0 HTSSIP_1.4.1    ggplot2_3.5.1   tidyr_1.3.1    
## [5] dplyr_1.1.4     rmarkdown_2.29 
## 
## loaded via a namespace (and not attached):
##  [1] ade4_1.7-22                 tidyselect_1.2.1           
##  [3] farver_2.1.2                Biostrings_2.75.1          
##  [5] fastmap_1.2.0               lazyeval_0.2.2             
##  [7] digest_0.6.37               lifecycle_1.0.4            
##  [9] cluster_2.1.6               survival_3.7-0             
## [11] magrittr_2.0.3              compiler_4.4.2             
## [13] rlang_1.1.4                 sass_0.4.9                 
## [15] tools_4.4.2                 igraph_2.1.1               
## [17] utf8_1.2.4                  yaml_2.3.10                
## [19] data.table_1.16.2           knitr_1.49                 
## [21] S4Arrays_1.7.1              labeling_0.4.3             
## [23] DelayedArray_0.33.2         plyr_1.8.9                 
## [25] abind_1.4-8                 BiocParallel_1.41.0        
## [27] withr_3.0.2                 purrr_1.0.2                
## [29] BiocGenerics_0.53.3         sys_3.4.3                  
## [31] grid_4.4.2                  stats4_4.4.2               
## [33] fansi_1.0.6                 multtest_2.63.0            
## [35] biomformat_1.35.0           colorspace_2.1-1           
## [37] Rhdf5lib_1.29.0             scales_1.3.0               
## [39] iterators_1.0.14            MASS_7.3-61                
## [41] SummarizedExperiment_1.37.0 cli_3.6.3                  
## [43] vegan_2.6-8                 crayon_1.5.3               
## [45] generics_0.1.3              httr_1.4.7                 
## [47] reshape2_1.4.4              ape_5.8                    
## [49] cachem_1.1.0                rhdf5_2.51.0               
## [51] stringr_1.5.1               zlibbioc_1.52.0            
## [53] splines_4.4.2               parallel_4.4.2             
## [55] XVector_0.47.0              matrixStats_1.4.1          
## [57] vctrs_0.6.5                 Matrix_1.7-1               
## [59] jsonlite_1.8.9              IRanges_2.41.1             
## [61] S4Vectors_0.45.2            coenocliner_0.2-3          
## [63] maketools_1.3.1             locfit_1.5-9.10            
## [65] foreach_1.5.2               jquerylib_0.1.4            
## [67] glue_1.8.0                  codetools_0.2-20           
## [69] stringi_1.8.4               gtable_0.3.6               
## [71] GenomeInfoDb_1.43.0         GenomicRanges_1.59.0       
## [73] UCSC.utils_1.3.0            munsell_0.5.1              
## [75] tibble_3.2.1                pillar_1.9.0               
## [77] htmltools_0.5.8.1           rhdf5filters_1.19.0        
## [79] GenomeInfoDbData_1.2.13     R6_2.5.1                   
## [81] doParallel_1.0.17           evaluate_1.0.1             
## [83] lattice_0.22-6              Biobase_2.67.0             
## [85] bslib_0.8.0                 Rcpp_1.0.13-1              
## [87] SparseArray_1.7.2           nlme_3.1-166               
## [89] permute_0.9-7               mgcv_1.9-1                 
## [91] DESeq2_1.47.1               xfun_0.49                  
## [93] MatrixGenerics_1.19.0       buildtools_1.0.0           
## [95] pkgconfig_2.0.3