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:

- Dataset exploration
  - Beta diversity ordinations of gradient fraction communities
  - Whole-community BD shift plots
- Dataset simulation
  - HTS-SIP dataset simulation
- Detecting incorporators
  - HR-SIP & MW-HR-SIP
  - q-SIP
- Quantifying the amount of isotope incorporated
  - quantifying incorporation

# HTS-SIP introduction

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:

- the microbial food web
  - e.g., Which bacteria contribute to each step of biomass breakdown in soil?
- taxonomic conservatism of ecological roles
  - e.g., Is taxonomic relatedness a good predictor of functional similarity (similar niche)?
  - e.g., Do OTU groupings accurately delimit functionally coherent individuals?
- *competition and mutualism *
  - e.g., Which microbes will (or will not) stabily coexist in a particular environment?
- how microbial physiology can dictate biogeographic patterns
  - e.g., Is the biogeography distribution of cellulose specialists more patchy than carbon source generalists?
- the prevalence of functional redundancy and community stability (resilience and resistance)
  - e.g., Do many microbes consume the same substrates? How can this affect stability of the microbial community and the ecosystem processes mediated by members of the community?

## Background

Note: I'm going to just refer to DNA-SIP, but most of this also applies to RNA-SIP.

### Stable Isotope Probing

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.

### HTS-SIP

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:

- DNA SIP
  - Neufeld, Josh D., Jyotsna Vohra, Marc G. Dumont, Tillmann Lueders, Mike Manefield, Michael W. Friedrich, and J. Colin Murrell. 2007. “DNA Stable-Isotope Probing.” Nature Protocols 2 (4): 860–66.
  - Buckley, Daniel H., Varisa Huangyutitham, Shi-Fang Hsu, and Tyrrell A. Nelson. 2007. “Stable Isotope Probing with 15N Achieved by Disentangling the Effects of Genome G+C Content and Isotope Enrichment on DNA Density.” Applied and Environmental Microbiology 73 (10): 3189–95.
  - Youngblut, Nicholas D., and Daniel H. Buckley. 2014. “Intra-Genomic Variation in G + C Content and Its Implications for DNA Stable Isotope Probing.” Environmental Microbiology Reports 6 (6): 767–75.
- HR-SIP
  - Pepe-Ranney, Charles, Chantal Koechli, Ruth Potrafka, Cheryl Andam, Erin Eggleston, Ferran Garcia-Pichel, and Daniel H. Buckley. 2015. “Non-Cyanobacterial Diazotrophs Mediate Dinitrogen Fixation in Biological Soil Crusts during Early Crust Formation.” The ISME Journal.
  - Pepe-Ranney, Charles, Ashley N. Campbell, Chantal N. Koechli, Sean Berthrong, and Daniel H. Buckley. 2016. “Unearthing the Ecology of Soil Microorganisms Using a High Resolution DNA-SIP Approach to Explore Cellulose and Xylose Metabolism in Soil.” Terrestrial Microbiology, 703.
- MW-HR-SIP
  - Youngblut, Nicholas D., and Daniel H. Buckley. “Evaluating the Accuracy of DNA Stable Isotope Probing”. in prep
- qSIP
  - Hungate, Bruce A., Rebecca L. Mau, Egbert Schwartz, J. Gregory Caporaso, Paul Dijkstra, Natasja van Gestel, Benjamin J. Koch, et al. 2015. “Quantitative Microbial Ecology Through Stable Isotope Probing.” Applied and Environmental Microbiology, August, AEM.02280–15.

# HTSSIP package data

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:         [ 10001 taxa and 139 samples ]
    ## sample_data() Sample Data:       [ 139 samples by 17 sample variables ]
    ## tax_table()   Taxonomy Table:    [ 10001 taxa by 8 taxonomic ranks ]
    ## phy_tree()    Phylogenetic Tree: [ 10001 tips and 10000 internal nodes ]

The “S2D2” part of the name stands for: 2 substrates assessed at 2 time points.

### Experiment overview

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:

1. 12C-Con, Day 3
2. 13C-Glu, Day 3
3. 13C-Cel, Day 3
4. 12C-Con, Day 14
5. 13C-Glu, Day 14
6. 13C-Cel, Day 14

Notes:

- “12C-Con” stands for the unlabeled control
- “13C-Glu” and “13C-Cel” stand for 13C-glucose and 13C-cellulose, respectively
- “Day” is the number of days from substrate addition to harvest

### Sample metadata

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>

## Parsing the dataset

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:

1. 12C-Con_Day3 vs 13C-Glu_Day3
2. 12C-Con_Day3 vs 13C-Cel_Day3
3. 12C-Con_Day14 vs 13C-Glu_Day14
4. 12C-Con_Day14 vs 13C-Cel_Day14

> 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:         [ 10001 taxa and 46 samples ]
    ## sample_data() Sample Data:       [ 46 samples by 17 sample variables ]
    ## tax_table()   Taxonomy Table:    [ 10001 taxa by 8 taxonomic ranks ]
    ## phy_tree()    Phylogenetic Tree: [ 10001 tips and 10000 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:         [ 10001 taxa and 46 samples ]
    ## sample_data() Sample Data:       [ 46 samples by 17 sample variables ]
    ## tax_table()   Taxonomy Table:    [ 10001 taxa by 8 taxonomic ranks ]
    ## phy_tree()    Phylogenetic Tree: [ 10001 tips and 10000 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:         [ 10001 taxa and 46 samples ]
    ## sample_data() Sample Data:       [ 46 samples by 17 sample variables ]
    ## tax_table()   Taxonomy Table:    [ 10001 taxa by 8 taxonomic ranks ]
    ## phy_tree()    Phylogenetic Tree: [ 10001 tips and 10000 internal nodes ]
    ##
    ## $`(Substrate=='12C-Con' & Day=='14') | (Substrate=='13C-Cel' & Day == '14')`
    ## phyloseq-class experiment-level object
    ## otu_table()   OTU Table:         [ 10001 taxa and 46 samples ]
    ## sample_data() Sample Data:       [ 46 samples by 17 sample variables ]
    ## tax_table()   Taxonomy Table:    [ 10001 taxa by 8 taxonomic ranks ]
    ## phy_tree()    Phylogenetic Tree: [ 10001 tips and 10000 internal nodes ]
    ##
    ## $`(Substrate=='12C-Con' & Day=='14') | (Substrate=='13C-Glu' & Day == '14')`
    ## phyloseq-class experiment-level object
    ## otu_table()   OTU Table:         [ 10001 taxa and 47 samples ]
    ## sample_data() Sample Data:       [ 47 samples by 17 sample variables ]
    ## tax_table()   Taxonomy Table:    [ 10001 taxa by 8 taxonomic ranks ]
    ## phy_tree()    Phylogenetic Tree: [ 10001 tips and 10000 internal nodes ]
    ##
    ## $`(Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Glu' & Day == '3')`
    ## phyloseq-class experiment-level object
    ## otu_table()   OTU Table:         [ 10001 taxa and 46 samples ]
    ## sample_data() Sample Data:       [ 46 samples by 17 sample variables ]
    ## tax_table()   Taxonomy Table:    [ 10001 taxa by 8 taxonomic ranks ]
    ## phy_tree()    Phylogenetic Tree: [ 10001 tips and 10000 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.

# References for methods in HTSSIP

- HTS-SIP data exploration

  - Multiple window High Resolution Stable Isotope Probing (MW-HR-SIP) as described in Youngblut and Buckley 2017

- HR-SIP

  - High Resolution Stable Isotope Probing (HR-SIP) as described in Pepe-Ranney et al., 2016

- MW-HR-SIP

  - Multiple window High Resolution Stable Isotope Probing (MW-HR-SIP) as described in Youngblut and Buckley 2017

- q-SIP

  - Quantitative stable isotope probing (q-SIP) as described in Hungate et al., 2015

# Session info

    sessionInfo()

    ## R version 3.4.2 (2017-09-28)
    ## Platform: x86_64-pc-linux-gnu (64-bit)
    ## Running under: Ubuntu 16.04.3 LTS
    ##
    ## Matrix products: default
    ## BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
    ## LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.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
    ##
    ## attached base packages:
    ## [1] stats     graphics  grDevices utils     datasets  methods   base
    ##
    ## other attached packages:
    ## [1] bindrcpp_0.2    phyloseq_1.20.0 HTSSIP_1.3.0    ggplot2_2.2.1
    ## [5] tidyr_0.7.2     dplyr_0.7.4
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] tidyselect_0.2.2    reshape2_1.4.2      purrr_0.2.4
    ##  [4] splines_3.4.2       lattice_0.20-35     rhdf5_2.20.0
    ##  [7] colorspace_1.3-2    stats4_3.4.2        mgcv_1.8-22
    ## [10] survival_2.41-3     rlang_0.1.2         glue_1.1.1
    ## [13] BiocGenerics_0.22.1 foreach_1.4.3       bindr_0.1
    ## [16] plyr_1.8.4          stringr_1.2.0       zlibbioc_1.22.0
    ## [19] Biostrings_2.44.2   munsell_0.4.3       gtable_0.2.0
    ## [22] codetools_0.2-15    evaluate_0.10.1     labeling_0.3
    ## [25] Biobase_2.36.2      knitr_1.17          permute_0.9-4
    ## [28] IRanges_2.10.5      biomformat_1.4.0    parallel_3.4.2
    ## [31] markdown_0.8        highr_0.6           Rcpp_0.12.13
    ## [34] scales_0.5.0        vegan_2.4-4         S4Vectors_0.14.7
    ## [37] jsonlite_1.5        XVector_0.16.0      mime_0.5
    ## [40] digest_0.6.12       stringi_1.1.5       grid_3.4.2
    ## [43] ade4_1.7-8          tools_3.4.2         magrittr_1.5
    ## [46] lazyeval_0.2.0      tibble_1.3.4        cluster_2.0.6
    ## [49] ape_4.1             pkgconfig_2.0.1     MASS_7.3-47
    ## [52] Matrix_1.2-11       data.table_1.10.4-2 assertthat_0.2.0
    ## [55] iterators_1.0.8     R6_2.2.2            multtest_2.32.0
    ## [58] igraph_1.1.2        nlme_3.1-131        compiler_3.4.2