HTSSIP is a package for analyzing high throughput sequence (HTS) data from DNA- or RNA-based stable isotope probing (SIP) experiments.
I recommend starting with this introductory vignette, then checking out the other vignettes:
The goal of HTS-SIP methods is to accurately determine which taxa (e.g., bacterial 16S rRNA OTUs) have incorporated isotopically labed substrate into biomass. By correctly identifying these 'incorporators', a number of fundamental questions about microbial ecology can be investigated, such as:
Note: I'm going to just refer to DNA-SIP, but most of this also applies to RNA-SIP.
In theory, if a taxon incorporates isotopically labeled substrate, then its DNA with be 'heavier' than if that same taxon had instead incorporated unlabeled substrate. By 'heavier', I mean that the DNA with have a higher buoyant density (BD) in a CsCl gradient. So, if you were to load isotopically labeled DNA into one CsCl gradient and unlabeled DNA (from the same taxon) in another CsCl gradient, then the labeled DNA with appear shifted to a heavier BD than the unlabeled DNA. The DNA should be more or less normally distributed in each gradient, so picture a Gaussian distribution that shifts from a 'light' BD to a 'heavy' BD between the control and treatment gradients. This is the basics of stable isotope probing, which is a powerful method because it can be used to detect which taxa in an entire community consumed particular substrates (e.g., cellulose).
If 'heavy' DNA is shifted from 'light' DNA, why are both a control and treatment gradient needed instead of just carefully separating out the 'heavy' DNA band from the 'light' DNA band in just one gradient? While this technique used be done in older SIP studies, the main problem is the BD of DNA is also a factor of G+C content (& some other factors). DNA with a higher G+C content has a 'heavier' BD. So, labeled, low G+C DNA can easily overlap with unlabeled high G+C DNA.
Let's consider a simple DNA-SIP experiment. A basic SIP experiment could consist of incubating aliquots of soil with either 13C-glucose (treatment) or 12C-glucose (the corresponding control). The DNA is then extracted from each soil sample, and heavy DNA is separated (at least partially) from light DNA using CsCl gradient ultra centrifugation. Now, the challenge is to accurately determine the location of each taxon's DNA in each of the 2 gradients in order to identify which taxa 'shifted' in the labeled treatment gradient due to isotope incorporation. Prior to the rise of high throughput sequencing (HTS) methods, this task was very difficult due to the limited throughput and taxonomic resolution of molecular fingerprinting and Sanger sequencing.
With high throughput sequencing, the taxonomic composition can be assessed at many points along a CsCl gradient by fractionating the gradient and sequencing the DNA in each fraction (note: I'm just going to focus on 16S rRNA sequencing, but fungal-ITS or metagenomic sequencing are also effective HTS-SIP methods). A large number of gradient fractions can be sequenced, which can help pinpoint the exact BD density range that each OTU occupies. For a typical HR-SIP experiment, ~20-30 fractions from each gradient are often sequenced.
Now remember, we are trying to use HTS-SIP to identify 'BD shifts', which indicate isotope incorporation. This is done by comparing labeled-treatment gradients to their corresponding unlabeled control gradients. Optimally, we'd be able to measure the absolute abundance of each taxon's OTU in each gradient fraction in order to get an idea where the Gaussian distribution of that OTU's DNA is in the gradient.
The major challenge of analyzing HTS data is that the data is compositional, in that the total number of sequences in a sample does not provide any information on the total numbers of each taxon in the sample. This is why relative abundances of OTUs are shown in most microbiome papers. Therefore, we have to try to determine the absolute abundance Gaussian distributions of each OTU in each gradient based on compositional 16S rRNA sequence data. Various HTS-SIP data analysis methods have be developed to address this challenge (e.g., HR-SIP and q-SIP). This R package was developed to easily implement on your own data, so that you can compare the results (and possible develop new methods). For more information on SIP and HTS-SIP, check out these references:
The phyloseq
R package provides some great functions for importing and analyzing microbiome data in R. So instead of reinventing the wheel, HTSSIP utilizes Phyloseq objects as input. See the GitHub page or the
package website for more info about the phyloseq.
The basic HTS-SIP dataset in the HTSSIP
package consists of three CsCl 3 gradients, each with ~24 samples (gradient fractions). There are 2 treatment gradients, which contained DNA extracted from microcosms recieving both cellulose and glucose, but only one substrate was 13C-labeled for each treatment. For example, the 13C-cellulose (13C-Cel) treatment recieved 13C-cellulose and 12C-glucose. This experimental design had the benefit of requiring only one 12C-control that recieved 12C-cellulose and 12C-glucose.
First, let's load some packages including HTSSIP
.
library(dplyr)
library(HTSSIP)
library(phyloseq)
Now, let's check out the HTSSIP phyloseq object that we will be using in the other vignettes.
physeq_S2D2
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1072 taxa and 139 samples ]
## sample_data() Sample Data: [ 139 samples by 17 sample variables ]
## tax_table() Taxonomy Table: [ 1072 taxa by 8 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1072 tips and 1071 internal nodes ]
The “S2D2” part of the name stands for: 2 substrates assessed at 2 time points.
The dataset is from a larger SIP experiment in which labeled substrates were incubated with 15 g of soil in microcosms. The microcosms were destructively sampled at certain numbers of days from the initial substrate addition. The DNA was extracted and loaded into CsCl gradients. Importantly, each microcosm recieved the same amount of all substrates used in the experiment, but only one of the substrates was 13C labeled for each labeled treatment, while the unlabeled controls recieved all 12C unlabeled versions. This design allowed for only one 12C control microcosm per time point (instead of 1 control per labeled substrate). Approximately 24 gradient fractions from each gradient were sequenced with MiSeq 16S rRNA sequencing.
The dataset subset that we'll be using consists of samples 6 CsCl gradients:
Important note: The dataset used in this example has been filtered such that all OTUs have over 300 reads. This is not a standard operation but has been done here to reduce the size of the dataset and the processing time.
Notes:
Here's a look at (some) of the sample metadata. Note that Bouyant_density
is
one of the columns. This is necessary for many of the analyses in HTSSIP}
physeq_S2D2 %>% sample_data %>% .[1:4,1:10]
## Sample Data: [4 samples by 10 sample variables]:
## X.Sample primer_number fwd_barcode
## 12C-Con.D14.R1_F17 12C-Con.D14.R1_F17 38 CGTGAGTG
## 13C-Cel.D3.R1_F21 13C-Cel.D3.R1_F21 7 GGATATCT
## 13C-Cel.D3.R1_F20 13C-Cel.D3.R1_F20 6 CGTGAGTG
## 13C-Cel.D14.R1_F12 13C-Cel.D14.R1_F12 57 ATCGTACG
## rev_barcode Substrate Day Microcosm_replicate Fraction
## 12C-Con.D14.R1_F17 TGAGTACG 12C-Con 14 1 17
## 13C-Cel.D3.R1_F21 CGAGAGTT 13C-Cel 3 1 21
## 13C-Cel.D3.R1_F20 CGAGAGTT 13C-Cel 3 1 20
## 13C-Cel.D14.R1_F12 CGAGCGAC 13C-Cel 14 1 12
## Buoyant_density Sample_type
## 12C-Con.D14.R1_F17 1.715475 unknown
## 13C-Cel.D3.R1_F21 1.705640 unknown
## 13C-Cel.D3.R1_F20 1.706186 unknown
## 13C-Cel.D14.R1_F12 1.735145 unknown
To check the formatting of the phyloseq object:
# Our phyloseq object should be formatted correctly (no errors)
physeq_S2D2 = physeq_format(physeq_S2D2)
physeq_S2D2 %>% sample_data %>% .[1:4,1:10]
## Sample Data: [4 samples by 10 sample variables]:
## X.Sample primer_number fwd_barcode
## 12C-Con.D14.R1_F17 12C-Con.D14.R1_F17 38 CGTGAGTG
## 13C-Cel.D3.R1_F21 13C-Cel.D3.R1_F21 7 GGATATCT
## 13C-Cel.D3.R1_F20 13C-Cel.D3.R1_F20 6 CGTGAGTG
## 13C-Cel.D14.R1_F12 13C-Cel.D14.R1_F12 57 ATCGTACG
## rev_barcode Substrate Day Microcosm_replicate Fraction
## 12C-Con.D14.R1_F17 TGAGTACG 12C-Con 14 1 17
## 13C-Cel.D3.R1_F21 CGAGAGTT 13C-Cel 3 1 21
## 13C-Cel.D3.R1_F20 CGAGAGTT 13C-Cel 3 1 20
## 13C-Cel.D14.R1_F12 CGAGCGAC 13C-Cel 14 1 12
## Buoyant_density Sample_type
## 12C-Con.D14.R1_F17 1.715475 unknown
## 13C-Cel.D3.R1_F21 1.705640 unknown
## 13C-Cel.D3.R1_F20 1.706186 unknown
## 13C-Cel.D14.R1_F12 1.735145 unknown
# This phyloseq object should NOT be formatted correctly
data(GlobalPatterns)
GlobalPatterns %>% sample_data %>% .[1:4,1:ncol(.)]
## Sample Data: [4 samples by 7 sample variables]:
## X.SampleID Primer Final_Barcode Barcode_truncated_plus_T
## CL3 CL3 ILBC_01 AACGCA TGCGTT
## CC1 CC1 ILBC_02 AACTCG CGAGTT
## SV1 SV1 ILBC_03 AACTGT ACAGTT
## M31Fcsw M31Fcsw ILBC_04 AAGAGA TCTCTT
## Barcode_full_length SampleType
## CL3 CTAGCGTGCGT Soil
## CC1 CATCGACGAGT Soil
## SV1 GTACGCACAGT Soil
## M31Fcsw TCGACATCTCT Feces
## Description
## CL3 Calhoun South Carolina Pine soil, pH 4.9
## CC1 Cedar Creek Minnesota, grassland, pH 6.1
## SV1 Sevilleta new Mexico, desert scrub, pH 8.3
## M31Fcsw M3, Day 1, fecal swab, whole body study
tryCatch(
physeq_format(GlobalPatterns),
error = function(e) e
)
## <simpleError in physeq_format(GlobalPatterns): No "Buoyant_density" column found in metadata>
For most of the analyses in HTSSIP
, the gradient fraction samples of each labeled-treatment gradient is compared to its corresponding control. For example, the dataset that we've been looking at has 4 comparisons:
Note that the same 12C-control is used for both treatments because both recieved glucose & cellulose, but only 1 substrate was isotopically labeled in each.
So to subset the data with the phyloseq
package, we can use the prune_samples()
function to subset the dataset:
m = sample_data(physeq_S2D2)
physeq_13C.Glu_D3 = prune_samples((m$Substrate=='12C-Con' & m$Day==3) | (m$Substrate=='13C-Glu' & m$Day==3), physeq_S2D2)
physeq_13C.Glu_D3
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1072 taxa and 46 samples ]
## sample_data() Sample Data: [ 46 samples by 17 sample variables ]
## tax_table() Taxonomy Table: [ 1072 taxa by 8 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1072 tips and 1071 internal nodes ]
physeq_13C.Cel_D14 = prune_samples((m$Substrate=='12C-Con' & m$Day==3) | (m$Substrate=='13C-Cel' & m$Day==14), physeq_S2D2)
physeq_13C.Cel_D14
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1072 taxa and 46 samples ]
## sample_data() Sample Data: [ 46 samples by 17 sample variables ]
## tax_table() Taxonomy Table: [ 1072 taxa by 8 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1072 tips and 1071 internal nodes ]
Writing a bunch of redundant code to make each subset is not a good method to subset the data if you have a bunch of subsets to make, so HTSSIP
includes some functions to help.
The expression used in the prune_samples()
to subset the data can be generalized as:
ex = "(Substrate=='12C-Con' & Day=='${Day}') | (Substrate=='${Substrate}' & Day == '${Day}')"
Now we need the parameters for each individual prune_samples()
call. Then, we can iteratively insert these parameters into the expression in order create all of the subsets. This can be accomplished with:
params = get_treatment_params(physeq_S2D2, c('Substrate', 'Day'), "Substrate != '12C-Con'")
params
## Substrate Day
## 1 13C-Cel 3
## 2 13C-Cel 14
## 3 13C-Glu 14
## 4 13C-Glu 3
Let's use the expression and parameters to subset the phyloseq object into a list of phyloseq object subsets:
physeq_S2D2_l = phyloseq_subset(physeq_S2D2, params, ex)
physeq_S2D2_l
## $`(Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Cel' & Day == '3')`
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1072 taxa and 46 samples ]
## sample_data() Sample Data: [ 46 samples by 17 sample variables ]
## tax_table() Taxonomy Table: [ 1072 taxa by 8 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1072 tips and 1071 internal nodes ]
##
## $`(Substrate=='12C-Con' & Day=='14') | (Substrate=='13C-Cel' & Day == '14')`
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1072 taxa and 46 samples ]
## sample_data() Sample Data: [ 46 samples by 17 sample variables ]
## tax_table() Taxonomy Table: [ 1072 taxa by 8 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1072 tips and 1071 internal nodes ]
##
## $`(Substrate=='12C-Con' & Day=='14') | (Substrate=='13C-Glu' & Day == '14')`
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1072 taxa and 47 samples ]
## sample_data() Sample Data: [ 47 samples by 17 sample variables ]
## tax_table() Taxonomy Table: [ 1072 taxa by 8 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1072 tips and 1071 internal nodes ]
##
## $`(Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Glu' & Day == '3')`
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1072 taxa and 46 samples ]
## sample_data() Sample Data: [ 46 samples by 17 sample variables ]
## tax_table() Taxonomy Table: [ 1072 taxa by 8 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1072 tips and 1071 internal nodes ]
With a list of phyloseq objects (1 list item per treatment-control comparison), we can use lapply()
or a similar function to run analyses on each subset of the dataset.
HTS-SIP data exploration
HR-SIP
MW-HR-SIP
q-SIP
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 stringi_1.4.3
## [43] grid_3.4.3 ade4_1.7-13 tools_3.4.3
## [46] magrittr_1.5 lazyeval_0.2.2 tibble_2.1.1
## [49] cluster_2.0.6 crayon_1.3.4 ape_5.3
## [52] pkgconfig_2.0.2 MASS_7.3-51.4 Matrix_1.2-17
## [55] data.table_1.10.4-3 assertthat_0.2.1 iterators_1.0.10
## [58] R6_2.4.0 multtest_2.34.0 igraph_1.2.4
## [61] nlme_3.1-131 compiler_3.4.3