Introduction

HTSSIP has some functionality for simulating basic HTS-SIP datasets. With that said, I recommend using a more sophisticated simulation toolset such as SIPSim for applications other than simple testing of HTS-SIP data analysis functions.

HTSSIP relies heavily on the great R package coenocliner. See this tutorial for a short and simple introduction.

Simulating a HTS-SIP dataset

In this vignette, we're going to simulate gradient fraction communities for 6 gradients, with the basic experimental setup as follows:

First, let's load some packages including HTSSIP.

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

OK, let's set the parameters needed for community simulations. We are basically going to follow the coenocliner tutorial, but instead of a transect along an environmental gradient, we are simulating communities in each fraction of a density gradient.

# setting parameters for tests
set.seed(318)                              # reproduciblility
M = 6                                      # number of OTUs (species)
ming = 1.67                                # gradient minimum...
maxg = 1.78                                # ...and maximum
nfrac = 24                                 # number of gradient fractions
locs = seq(ming, maxg, length=nfrac)       # Fraction BD's
tol  = rep(0.005, M)                       # species tolerances
h    = ceiling(rlnorm(M, meanlog=11))      # max abundances

opt = rnorm(M, mean=1.71, sd=0.008)        # species optima (drawn from a normal dist.)
params = cbind(opt=opt, tol=tol, h=h)      # put in a matrix

With the current parameters, we can simulate the gradient fraction communities for 1 density gradient:

df_OTU = gradient_sim(locs, params)
df_OTU
##    OTU.1 OTU.2 OTU.3 OTU.4 OTU.5  OTU.6 Buoyant_density Fraction
## 1      0     0     0     0     0      0        1.670000        1
## 2      0     0     0     0     0      0        1.674783        2
## 3      0     0     5     0     0      0        1.679565        3
## 4      0     0   131     0     0      0        1.684348        4
## 5      0     3  2074     4    14     22        1.689130        5
## 6      0   134 13016   131   491    727        1.693913        6
## 7     15  1891 32630  1177  6248   9066        1.698696        7
## 8    408 11003 33250  4449 31193  48158        1.703478        8
## 9   3665 25564 13439  6521 62709 100021        1.708261        9
## 10 12256 23561  2218  3901 50332  83169        1.713043       10
## 11 16510  8423   151   948 16405  28235        1.717826       11
## 12  8946  1236     3   108  2052   3800        1.722609       12
## 13  1953    73     0     3   102    199        1.727391       13
## 14   163     1     0     0     4      5        1.732174       14
## 15     5     0     0     0     0      0        1.736957       15
## 16     0     0     0     0     0      0        1.741739       16
## 17     0     0     0     0     0      0        1.746522       17
## 18     0     0     0     0     0      0        1.751304       18
## 19     0     0     0     0     0      0        1.756087       19
## 20     0     0     0     0     0      0        1.760870       20
## 21     0     0     0     0     0      0        1.765652       21
## 22     0     0     0     0     0      0        1.770435       22
## 23     0     0     0     0     0      0        1.775217       23
## 24     0     0     0     0     0      0        1.780000       24

As you can see, the abundance distribution of each OTU is approximately Gaussian, with varying optima among OTUs.

Simulating the full SIP experiment

Let's make the communities of all gradient fractions for all gradients.

If all OTUs in the 13C-treatment incorporated labeled isotope, then their abundance distributions should be shifted to 'heavier' buoyant densities. Let's set the 13C-treatment gradients to have a higher mean species optima. For kicks, let's also increase the species optima variance (representing more variable isotope incorporation percentages).

opt1 = rnorm(M, mean=1.7, sd=0.005)      # species optima
params1 = cbind(opt=opt1, tol=tol, h=h)  # put in a matrix
opt2 = rnorm(M, mean=1.7, sd=0.005)     
params2 = cbind(opt=opt2, tol=tol, h=h) 
opt3 = rnorm(M, mean=1.7, sd=0.005)     
params3 = cbind(opt=opt3, tol=tol, h=h) 
opt4 = rnorm(M, mean=1.72, sd=0.008)    
params4 = cbind(opt=opt4, tol=tol, h=h) 
opt5 = rnorm(M, mean=1.72, sd=0.008)    
params5 = cbind(opt=opt5, tol=tol, h=h) 
opt6 = rnorm(M, mean=1.72, sd=0.008)    
params6 = cbind(opt=opt6, tol=tol, h=h) 

# we need a named list of parameters for the next step. The names in the list will be used as sample IDs
params_all = list(
  '12C-Con_rep1' = params1,
  '12C-Con_rep2' = params2,
  '12C-Con_rep3' = params3,
  '13C-Glu_rep1' = params4,
  '13C-Glu_rep2' = params5,
  '13C-Glu_rep3' = params6
)

Additional sample metadata

Additional sample metadata can be added to the resulting phyloseq object that we are going to simulate. To add metadata, just make a data.frame object with a Gradient column, which will be used for matching to the simulated samples. The Gradient column values should match the names in the parameter list.

meta = data.frame(
  'Gradient' = c('12C-Con_rep1', '12C-Con_rep2', '12C-Con_rep3',
                 '13C-Glu_rep1', '13C-Glu_rep2', '13C-Glu_rep3'),
  'Treatment' = c(rep('12C-Con', 3), rep('13C-Glu', 3)),
  'Replicate' = c(1:3, 1:3)
)

# do the names match?
all(meta$Gradient %in% names(params_all))
## [1] TRUE

OK. Let's make the phyloseq object with the parameters that we specified above.

## physeq object
physeq_rep3 = HTSSIP_sim(locs, params_all, meta=meta)
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 ]

How does the sample_data table look?

phyloseq::sample_data(physeq_rep3) %>% head
## Sample Data:        [6 samples by 5 sample variables]:
##                             Gradient Buoyant_density Fraction Treatment
## 12C-Con_rep1_1.668185_1 12C-Con_rep1        1.668185        1   12C-Con
## 12C-Con_rep1_1.680254_2 12C-Con_rep1        1.680254        2   12C-Con
## 12C-Con_rep1_1.679431_3 12C-Con_rep1        1.679431        3   12C-Con
## 12C-Con_rep1_1.682305_4 12C-Con_rep1        1.682305        4   12C-Con
## 12C-Con_rep1_1.689408_5 12C-Con_rep1        1.689408        5   12C-Con
## 12C-Con_rep1_1.691388_6 12C-Con_rep1        1.691388        6   12C-Con
##                         Replicate
## 12C-Con_rep1_1.668185_1         1
## 12C-Con_rep1_1.680254_2         1
## 12C-Con_rep1_1.679431_3         1
## 12C-Con_rep1_1.682305_4         1
## 12C-Con_rep1_1.689408_5         1
## 12C-Con_rep1_1.691388_6         1

Simulating qPCR data

The q-SIP analysis requires qPCR measurements of gene copies for each gradient fraction. We can simulate this data for the HTS-SIP dataset phyloseq object that we just created.

The qPCR value simulation function is fairly flexible. For input, functions are provided that define how qPCR values (averages & variability) relate to buoyant density. For example, you can set the 'peak' of qPCR values to be located at a BD of 1.7 for unlabeled gradient samples and a 'peak' at a BD of 1.75 for labeled samples. For this example, let's set the error among replicate qPCR reactions (technical replicates) to scale with the mean qPCR values for that gradient sample.

# unlabeled control qPCR values
control_mean_fun = function(x) dnorm(x, mean=1.70, sd=0.01) * 1e8
# error will scale with the mean 
control_sd_fun = function(x) control_mean_fun(x) / 3
# labeled treatment values
treat_mean_fun = function(x) dnorm(x, mean=1.75, sd=0.01) * 1e8
# error will scale with the mean 
treat_sd_fun = function(x) treat_mean_fun(x) / 3

OK. Let's simulate the qPCR data.

physeq_rep3_qPCR = qPCR_sim(physeq_rep3,
                control_expr='Treatment=="12C-Con"',
                control_mean_fun=control_mean_fun,
                control_sd_fun=control_sd_fun,
                treat_mean_fun=treat_mean_fun,
                treat_sd_fun=treat_sd_fun)

physeq_rep3_qPCR %>% names
## [1] "raw"     "summary"

The output is a list containing:

physeq_rep3_qPCR$raw %>% head
##   Buoyant_density IS_CONTROL qPCR_tech_rep1 qPCR_tech_rep2 qPCR_tech_rep3
## 1        1.668185       TRUE       38380656       37801138       22008529
## 2        1.680254       TRUE      560887106      659743071      354651067
## 3        1.679431       TRUE      618877197      770165585      538409720
## 4        1.682305       TRUE      554753523     1190798381     1037082632
## 5        1.689408       TRUE     3312072961     1398862657      707835246
## 6        1.691388       TRUE     1515455205     3901629542     2667702680
##                    Sample     Gradient Fraction Treatment Replicate
## 1 12C-Con_rep1_1.668185_1 12C-Con_rep1        1   12C-Con         1
## 2 12C-Con_rep1_1.680254_2 12C-Con_rep1        2   12C-Con         1
## 3 12C-Con_rep1_1.679431_3 12C-Con_rep1        3   12C-Con         1
## 4 12C-Con_rep1_1.682305_4 12C-Con_rep1        4   12C-Con         1
## 5 12C-Con_rep1_1.689408_5 12C-Con_rep1        5   12C-Con         1
## 6 12C-Con_rep1_1.691388_6 12C-Con_rep1        6   12C-Con         1
physeq_rep3_qPCR$summary %>% head
##   IS_CONTROL                  Sample Buoyant_density qPCR_tech_rep_mean
## 1      FALSE 13C-Glu_rep1_1.671712_2        1.671712       1.953729e-04
## 2      FALSE 13C-Glu_rep1_1.671722_1        1.671722       2.710671e-04
## 3      FALSE 13C-Glu_rep1_1.680311_3        1.680311       1.368177e-01
## 4      FALSE 13C-Glu_rep1_1.683540_4        1.683540       1.039858e+00
## 5      FALSE 13C-Glu_rep1_1.688831_5        1.688831       2.230201e+01
## 6      FALSE 13C-Glu_rep1_1.696594_7        1.696594       2.143897e+03
##   qPCR_tech_rep_sd     Gradient Fraction Treatment Replicate
## 1     5.435725e-05 13C-Glu_rep1        2   13C-Glu         1
## 2     3.818031e-05 13C-Glu_rep1        1   13C-Glu         1
## 3     6.395073e-02 13C-Glu_rep1        3   13C-Glu         1
## 4     6.074444e-01 13C-Glu_rep1        4   13C-Glu         1
## 5     1.446576e+01 13C-Glu_rep1        5   13C-Glu         1
## 6     6.352483e+02 13C-Glu_rep1        7   13C-Glu         1

Let's plot the data.

x_lab = bquote('Buoyant density (g '* ml^-1*')')
y_lab = '16S rRNA gene copies'

ggplot(physeq_rep3_qPCR$summary, aes(Buoyant_density, qPCR_tech_rep_mean,
                      ymin=qPCR_tech_rep_mean-qPCR_tech_rep_sd,
                      ymax=qPCR_tech_rep_mean+qPCR_tech_rep_sd,
                      color=IS_CONTROL, shape=Replicate)) +
  geom_pointrange() +
  scale_color_discrete('Unlabeled\ncontrol') +
  labs(x=x_lab, y=y_lab) +
  theme_bw()

plot of chunk unnamed-chunk-12

With this simulation, we made the separation between labeled and unlabeled DNA pretty easy to distinguish.

OK. Just to show the flexiblity of the qPCR value simulation function, let's try using some other functions as input.

# using the Cauchy distribution instead of normal distributions
control_mean_fun = function(x) dcauchy(x, location=1.70, scale=0.01) * 1e8
control_sd_fun = function(x) control_mean_fun(x) / 3
treat_mean_fun = function(x) dcauchy(x, location=1.75, scale=0.01) * 1e8
treat_sd_fun = function(x) treat_mean_fun(x) / 3
# simulating qPCR values
physeq_rep3_qPCR = qPCR_sim(physeq_rep3,
                control_expr='Treatment=="12C-Con"',
                control_mean_fun=control_mean_fun,
                control_sd_fun=control_sd_fun,
                treat_mean_fun=treat_mean_fun,
                treat_sd_fun=treat_sd_fun)

Now, how does the data look?

ggplot(physeq_rep3_qPCR$summary, aes(Buoyant_density, qPCR_tech_rep_mean,
                      ymin=qPCR_tech_rep_mean-qPCR_tech_rep_sd,
                      ymax=qPCR_tech_rep_mean+qPCR_tech_rep_sd,
                      color=IS_CONTROL, shape=Replicate)) +
  geom_pointrange() +
  scale_color_discrete('Unlabeled\ncontrol') +
  labs(x=x_lab, y=y_lab) +
  theme_bw()

plot of chunk unnamed-chunk-14

Session info

sessionInfo()
## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] phyloseq_1.22.3 HTSSIP_1.4.1    ggplot2_3.2.0   tidyr_0.8.3    
## [5] dplyr_0.8.0.1  
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_0.2.5    reshape2_1.4.3      purrr_0.3.2        
##  [4] splines_3.4.3       lattice_0.20-38     rhdf5_2.22.0       
##  [7] colorspace_1.4-1    stats4_3.4.3        mgcv_1.8-28        
## [10] survival_2.44-1.1   rlang_0.4.0         pillar_1.4.1       
## [13] glue_1.3.1          withr_2.1.2         BiocGenerics_0.24.0
## [16] foreach_1.4.4       plyr_1.8.4          stringr_1.4.0      
## [19] zlibbioc_1.24.0     Biostrings_2.46.0   munsell_0.5.0      
## [22] gtable_0.3.0        codetools_0.2-16    evaluate_0.14      
## [25] labeling_0.3        Biobase_2.38.0      knitr_1.18         
## [28] permute_0.9-5       IRanges_2.12.0      biomformat_1.6.0   
## [31] parallel_3.4.3      markdown_0.9        highr_0.8          
## [34] Rcpp_1.0.1          scales_1.0.0        vegan_2.5-1        
## [37] S4Vectors_0.16.0    jsonlite_1.6        XVector_0.18.0     
## [40] mime_0.6            digest_0.6.19       coenocliner_0.2-2  
## [43] stringi_1.4.3       grid_3.4.3          ade4_1.7-13        
## [46] tools_3.4.3         magrittr_1.5        lazyeval_0.2.2     
## [49] tibble_2.1.1        cluster_2.0.6       crayon_1.3.4       
## [52] ape_5.3             pkgconfig_2.0.2     MASS_7.3-51.4      
## [55] Matrix_1.2-17       data.table_1.10.4-3 assertthat_0.2.1   
## [58] iterators_1.0.10    R6_2.4.0            multtest_2.34.0    
## [61] igraph_1.2.4        nlme_3.1-131        compiler_3.4.3