# q-SIP

#### Nick Youngblut

#### 2016-11-01

# Introduction

## q-SIP method workflow

Quantitative stable isotope probing (q-SIP) attempts to deal with the compositional problems of HTS-SIP data by transforming relative abundances with qPCR counts of total gene copies.

More details can be found at:

> 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` can both simulate q-SIP datasets and run the q-SIP analysis (with parallel computation of bootstrap confidence intervals).

# 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 4 sample variables ]

The associated qPCR data is, in this case a list of length = 2.

    # qPCR data (list object)
    physeq_rep3_qPCR %>% names

    ## [1] "raw"     "summary"

Really, only the `summary` table is needed for the next step. You could just use the `summary` data.frame object and that would be fine. However, you should have just 1 qPCR value per sample. So, if you have ‘technical’ qPCR replicates (eg., triplicate qPCR reactions per DNA sample), then you need to summarize those replicates to get one value (eg., the mean of all replicate values).

    physeq_rep3_qPCR$raw %>% head(n=3)

    ##   Buoyant_density IS_CONTROL qPCR_tech_rep1 qPCR_tech_rep2 qPCR_tech_rep3
    ## 1        1.668185       TRUE      374460160      326925692      191209173
    ## 2        1.680254       TRUE      418974851      846274648     1057133579
    ## 3        1.679431       TRUE      456070748      431437845      307250348
    ##                  Sample     Gradient Treatment Replicate
    ## 1 12C-Con_rep1_1.668185 12C-Con_rep1   12C-Con         1
    ## 2 12C-Con_rep1_1.680254 12C-Con_rep1   12C-Con         1
    ## 3 12C-Con_rep1_1.679431 12C-Con_rep1   12C-Con         1

    physeq_rep3_qPCR$summary %>% head(n=3)

    ##   IS_CONTROL                Sample Buoyant_density qPCR_tech_rep_mean
    ## 1      FALSE 13C-Glu_rep1_1.671712        1.671712           46598172
    ## 2      FALSE 13C-Glu_rep1_1.671722        1.671722           42299449
    ## 3      FALSE 13C-Glu_rep1_1.680311        1.680311           53536519
    ##   qPCR_tech_rep_sd     Gradient Treatment Replicate
    ## 1          9055833 13C-Glu_rep1   13C-Glu         1
    ## 2         16369298 13C-Glu_rep1   13C-Glu         1
    ## 3         18265667 13C-Glu_rep1   13C-Glu         1

# q-SIP

## Abundance count transformation

First, we’ll transform the OTU counts. The following function will do the following:

- Total sum scaling (TSS) on OTU counts (relative abundances)
  - It’s OK if OTU counts are already TSS-transformed
- Mapping qPCR data to the OTU table
- The `sample_idx` parameter determines how the two datasets are matched to each other.
- There should be
- Multiplying the counts by the qPCR values.

    phyloseq::otu_table(physeq_rep3) %>% .[1:5, 1:5]

    ## OTU Table:          [5 taxa and 5 samples]
    ##                      taxa are rows
    ##       12C-Con_rep1_1.668185 12C-Con_rep1_1.680254 12C-Con_rep1_1.679431
    ## OTU.1                    57                   925                  5809
    ## OTU.2                     0                     1                    73
    ## OTU.3                    19                   492                  5078
    ## OTU.4                     0                     0                     0
    ## OTU.5                     0                     0                     3
    ##       12C-Con_rep1_1.682305 12C-Con_rep1_1.689408
    ## OTU.1                 14645                 15089
    ## OTU.2                  1077                  7499
    ## OTU.3                 21603                 36426
    ## OTU.4                     0                    11
    ## OTU.5                    43                  1017

    physeq_rep3_t = OTU_qPCR_trans(physeq_rep3, physeq_rep3_qPCR)
    phyloseq::otu_table(physeq_rep3_t) %>% .[1:5, 1:5]

    ## OTU Table:          [5 taxa and 5 samples]
    ##                      taxa are rows
    ##       12C-Con_rep1_1.668185 12C-Con_rep1_1.680254 12C-Con_rep1_1.679431
    ## OTU.1             223148756             504984567             207913325
    ## OTU.2                     0                545929               2612786
    ## OTU.3              74382919             268597197             181749675
    ## OTU.4                     0                     0                     0
    ## OTU.5                     0                     0                107375
    ##       12C-Con_rep1_1.682305 12C-Con_rep1_1.689408
    ## OTU.1             252808193             327552736
    ## OTU.2              18591630             162788652
    ## OTU.3             372920136             790737356
    ## OTU.4                     0                238789
    ## OTU.5                742284              22077085

## BD shift (Z) and atom fraction excess (A)

With the transformed OTU abundance counts, let’s calculate the BD shift (Z) and atom fraction excess (A) for each OTU. We’ll be comparing controls (designated by the `control_expr` option) and the labeled treatment (all samples not identified by the `control_expr` expression). In this case the `control_expr` is set to `Treatment=="12C-Con"`, which will select all samples where the `Treatment` column in the phyloseq `sample_data` table equals `12C-Con`.

The `treatment_rep` option designates which column in the `sample_data` table defines the replicate gradients (eg., `12C-Con_rep1`, `12C-Con_rep2`, etc).

    # The 'Treatment' column designates treatment vs control
    # The 'Replicate' column designates treatment replicates
    physeq_rep3 %>% sample_data %>% dplyr::select(Treatment, Replicate) %>% distinct

    ## Sample Data:        [6 samples by 2 sample variables]:
    ##   Treatment Replicate
    ## 1   12C-Con         1
    ## 2   12C-Con         2
    ## 3   12C-Con         3
    ## 4   13C-Glu         1
    ## 5   13C-Glu         2
    ## 6   13C-Glu         3

    atomX = qSIP_atom_excess(physeq_rep3_t,
                             control_expr='Treatment=="12C-Con"',
                             treatment_rep='Replicate')
    atomX %>% names

    ## [1] "W" "A"

The output is a list of data.frames. The `W` table lists the weighted mean BD for each OTU in each gradient. This can be considered a sort-of ‘raw data’ table. The real meat of the output is the `A` table, which lists (among other things) the BD shift (Z) and atom fraction excess (A) as defined by Hungate et al., (2015).

> Note: a more complicated expression can be used to designate control samples. For example: ‘Treatment == “12C-Con” & Day == 1’

## Bootstrap confidence intervals

OK. The last step is to calculate atom % excess confidence intervals (CI). This is necessary for getting some sort of indication on how accurate the atom fraction excess estimations are. Also, Hungate et al., (2015) uses the CIs to call incorporators, in which an OTU was considered an incorporator if the CI was greater than, and did not span, zero.

We’ll use 20 bootstrap replicates for this tutorial, but I recommend 100 for a real analysis. Bootstrap replicates can be run in parallel (see the function documentation).

    df_atomX_boot = qSIP_bootstrap(atomX, n_boot=20)
    df_atomX_boot %>% head

    ## # A tibble: 6 x 11
    ##     OTU     Wlab   Wlight          Z        Gi   Mlight Mheavymax     Mlab
    ##   <chr>    <dbl>    <dbl>      <dbl>     <dbl>    <dbl>     <dbl>    <dbl>
    ## 1 OTU.1 1.721357 1.699178 0.02217880 0.6361398 308.0065  317.6638 312.0268
    ## 2 OTU.2 1.706713 1.694977 0.01173549 0.5858298 307.9816  317.6640 310.1139
    ## 3 OTU.3 1.731520 1.696813 0.03470620 0.6078171 307.9925  317.6639 314.2921
    ## 4 OTU.4 1.730946 1.704892 0.02605395 0.7045621 308.0405  317.6636 312.7479
    ## 5 OTU.5 1.729513 1.705316 0.02419674 0.7096354 308.0430  317.6636 312.4138
    ## 6 OTU.6 1.727713 1.699375 0.02833873 0.6384898 308.0077  317.6638 313.1440
    ## # ... with 3 more variables: A <dbl>, A_CI_low <dbl>, A_CI_high <dbl>

Let’s use the same threshold as Hungate et al., (2015) for defining incorporators.

    CI_threshold = 0
    df_atomX_boot = df_atomX_boot %>%
      mutate(Incorporator = A_CI_low > CI_threshold,
                    OTU = reorder(OTU, -A))

How many incorporators?

    n_incorp = df_atomX_boot %>%
      filter(Incorporator == TRUE) %>%
      nrow
    cat('Number of incorporators:', n_incorp, '\n')

    ## Number of incorporators: 6

It’s time to plot the results!

    ggplot(df_atomX_boot, aes(OTU, A, ymin=A_CI_low, ymax=A_CI_high, color=Incorporator)) +
      geom_pointrange(size=0.25) +
      geom_linerange() +
      geom_hline(yintercept=0, linetype='dashed', alpha=0.5) +
      labs(x='OTU', y='Atom fraction excess') +
      theme_bw() +
      theme(
        axis.text.x = element_blank()
      )

# Session info

    sessionInfo()

    ## R version 3.2.4 Revised (2016-03-16 r70336)
    ## Platform: x86_64-pc-linux-gnu (64-bit)
    ## Running under: Ubuntu precise (12.04.5 LTS)
    ##
    ## 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] phyloseq_1.14.0 HTSSIP_1.0.3    ggplot2_2.1.0   tidyr_0.5.1
    ## [5] dplyr_0.5.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] Biobase_2.30.0             splines_3.2.4
    ##  [3] foreach_1.4.3              Formula_1.2-1
    ##  [5] assertthat_0.1             stats4_3.2.4
    ##  [7] latticeExtra_0.6-28        yaml_2.1.13
    ##  [9] RSQLite_1.0.0              lattice_0.20-33
    ## [11] chron_2.3-47               digest_0.6.10
    ## [13] GenomicRanges_1.22.4       RColorBrewer_1.1-2
    ## [15] XVector_0.10.0             colorspace_1.2-6
    ## [17] htmltools_0.3.5            Matrix_1.2-6
    ## [19] plyr_1.8.4                 XML_3.98-1.4
    ## [21] DESeq2_1.10.1              genefilter_1.52.1
    ## [23] zlibbioc_1.16.0            xtable_1.8-2
    ## [25] scales_0.4.0               BiocParallel_1.4.3
    ## [27] annotate_1.48.0            tibble_1.1
    ## [29] mgcv_1.8-12                IRanges_2.4.8
    ## [31] SummarizedExperiment_1.0.2 nnet_7.3-12
    ## [33] BiocGenerics_0.16.1        lazyeval_0.2.0
    ## [35] survival_2.39-5            RJSONIO_1.3-0
    ## [37] magrittr_1.5               evaluate_0.9
    ## [39] doParallel_1.0.10          nlme_3.1-128
    ## [41] MASS_7.3-45                RcppArmadillo_0.7.100.3.1
    ## [43] foreign_0.8-66             vegan_2.4-0
    ## [45] tools_3.2.4                data.table_1.9.6
    ## [47] formatR_1.4                stringr_1.0.0
    ## [49] S4Vectors_0.8.11           locfit_1.5-9.1
    ## [51] munsell_0.4.3              cluster_2.0.4
    ## [53] AnnotationDbi_1.32.3       lambda.r_1.1.7
    ## [55] Biostrings_2.38.4          ade4_1.7-4
    ## [57] compiler_3.2.4             GenomeInfoDb_1.6.3
    ## [59] futile.logger_1.4.1        grid_3.2.4
    ## [61] coenocliner_0.2-2          iterators_1.0.8
    ## [63] biom_0.3.12                igraph_1.0.1
    ## [65] labeling_0.3               rmarkdown_1.0
    ## [67] gtable_0.2.0               codetools_0.2-14
    ## [69] multtest_2.26.0            DBI_0.4-1
    ## [71] reshape2_1.4.1             R6_2.1.2
    ## [73] gridExtra_2.2.1            knitr_1.13
    ## [75] Hmisc_3.17-4               futile.options_1.0.0
    ## [77] permute_0.9-0              ape_3.5
    ## [79] stringi_1.1.1              parallel_3.2.4
    ## [81] Rcpp_0.12.5                geneplotter_1.48.0
    ## [83] rpart_4.1-10               acepack_1.3-3.3