# MW-HR-SIP

#### Nick Youngblut

#### 2016-11-01

# Introduction

## HR-SIP method workflow:

- Remove sparsely occuring OTUs under the assumption that they are sequencing errors.
  - This is assuming that DNA from and OTU is likely going to be found in multiple gradient fractions due to differences in G+C content, diffusion, etc.
  - So, if an OTU is just detected in only one fraction of a gradient, we can’t tell whether it is a sequencing artifact or a real, but rare OTU.
- For each OTU, use DESeq2 to detect significant enrichment in ‘heavy’ gradient fractions.
  - ‘heavy’ gradients are selected ad hoc
  - Multiple hypotheses are adjusted for.
  - Significantly enriched OTUs are considered to have incorporated isotope (incorporators)

## MW-HR-SIP method workflow:

- The method is similar to HR-SIP, except instead of 1 ‘heavy’ BD window used, multiple windows are used, but the tradeoff is the higher number of multiple hypotheses that must be corrected for.
  - For each OTU, the BD window with highest log2 fold change (l2fc) value is used, so the BD window can vary among OTUs.
- This method has been shown to be generally more sensitive (more true positives) than HR-SIP and more specific (less false positives) than q-SIP. See the manuscript “Multiple window High Resolution Stable Isotope Probing (MW-HR-SIP)” for more information.

# Dataset

First, let’s load some packages including `HTSSIP`.

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

See HTSSIP introduction vignette for a description on why dataset parsing (all treatment-control comparisons) is needed.

Let’s see the already parsed dataset

    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 ]

# HR-SIP

Let’s set some parameters used later.

    # adjusted P-value cutoff
    padj_cutoff = 0.1
    # number of cores for parallel processing (increase depending on your computational hardware)
    ncores = 2

## One treatment-control comparison

First, we’ll just run HR-SIP on 1 treatment-control comparison. Let’s get the individual phyloseq object.

    physeq = physeq_S2D2_l[[1]]
    physeq

    ## 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 ]

Let’s check that the samples belong to either a 13C-treatment or 12C-control.

    physeq %>% sample_data %>% .$Substrate %>% table

    ## .
    ## 12C-Con 13C-Cel
    ##      23      23

OK, we should be ready to run HR-SIP!

Note that the `design` parameter for `HRSIP()` is the experimental design parameter for calculating log2 fold change (l2fc) values with DESeq. Here, it’s used to distinguish label-treatment and unlabel-control samples.

    df_l2fc = HRSIP(physeq,
                    design = ~Substrate,
                    padj_cutoff = padj_cutoff,
                    sparsity_threshold = c(0,0.15,0.3))   # just using 3 thresholds to reduce time

    ## Sparsity threshold: 0
    ## Density window: 1.7-1.75
    ## Sparsity threshold: 0.15
    ## Density window: 1.7-1.75
    ## Sparsity threshold: 0.3
    ## Density window: 1.7-1.75
    ## Sparsity threshold with the most rejected hypotheses: 0.3

    df_l2fc %>% head(n=3)

    ## # A tibble: 3 x 17
    ##        OTU log2FoldChange         p  padj    Rank1            Rank2
    ##      <chr>          <dbl>     <dbl> <dbl>   <fctr>           <fctr>
    ## 1  OTU.514     -0.2504081 0.8504578     1 Bacteria __Proteobacteria
    ## 2  OTU.816     -0.5615771 0.9137773     1 Bacteria __Proteobacteria
    ## 3 OTU.1099     -0.2104345 0.7761086     1 Bacteria  __Acidobacteria
    ## # ... with 11 more variables: Rank3 <fctr>, Rank4 <fctr>, Rank5 <fctr>,
    ## #   Rank6 <fctr>, Rank7 <fctr>, Rank8 <fctr>, density_min <dbl>,
    ## #   density_max <dbl>, sparsity_threshold <dbl>, sparsity_apply <chr>,
    ## #   l2fc_threshold <dbl>

How many “incorporators”" (rejected hypotheses)?

    df_l2fc %>%
      filter(padj < padj_cutoff) %>%
      group_by() %>%
      summarize(n_incorp_OTUs = OTU %>% unique %>% length)

    ## # A tibble: 1 x 1
    ##   n_incorp_OTUs
    ##           <int>
    ## 1             8

Let’s plot a breakdown of incorporators for each phylum.

    # summarizing
    df_l2fc_s = df_l2fc %>%
      filter(padj < padj_cutoff) %>%
      mutate(Rank2 = gsub('^__', '', Rank2)) %>%
      group_by(Rank2) %>%
      summarize(n_incorp_OTUs = OTU %>% unique %>% length) %>%
      ungroup()

    # plotting
    ggplot(df_l2fc_s, aes(Rank2, n_incorp_OTUs)) +
        geom_bar(stat='identity') +
        labs(x='Phylum', y='Number of incorporators') +
        theme_bw() +
        theme(
          axis.text.x = element_text(angle=45, hjust=1)
        )

## All treatment-control comparisons

Let’s now run HR-SIP on all treatment-control comparisons in the dataset:

    # Number of comparisons
    physeq_S2D2_l %>% length

    ## [1] 4

The function `plyr::ldply()` is useful (compared to `lapply()`) beccause it can be run in parallel and returns a data.frame object.

    # Running in parallel; you may need to change the backend for your environment.
    # Or you can just set .parallel=FALSE.
    doParallel::registerDoParallel(ncores)

    df_l2fc = plyr::ldply(physeq_S2D2_l,
                          HRSIP,
                          design = ~Substrate,
                          padj_cutoff = padj_cutoff,
                          sparsity_threshold = c(0,0.15,0.3),  # just using 3 thresholds to reduce run time
                          .parallel=TRUE)
    df_l2fc %>% head(n=3)

    ##                                                                       .id
    ## 1 (Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Cel' & Day == '3')
    ## 2 (Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Cel' & Day == '3')
    ## 3 (Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Cel' & Day == '3')
    ##        OTU log2FoldChange         p padj    Rank1            Rank2
    ## 1  OTU.514     -0.2504081 0.8504578    1 Bacteria __Proteobacteria
    ## 2  OTU.816     -0.5615771 0.9137773    1 Bacteria __Proteobacteria
    ## 3 OTU.1099     -0.2104345 0.7761086    1 Bacteria  __Acidobacteria
    ##                   Rank3                  Rank4            Rank5
    ## 1 __Deltaproteobacteria    __Desulfobacterales __Nitrospinaceae
    ## 2 __Deltaproteobacteria    __Desulfobacterales __Nitrospinaceae
    ## 3               __32-21 __uncultured_bacterium             <NA>
    ##          Rank6                  Rank7 Rank8 density_min density_max
    ## 1 __uncultured __uncultured_bacterium  <NA>         1.7        1.75
    ## 2 __uncultured __uncultured_bacterium  <NA>         1.7        1.75
    ## 3         <NA>                   <NA>  <NA>         1.7        1.75
    ##   sparsity_threshold sparsity_apply l2fc_threshold
    ## 1                0.3            all           0.25
    ## 2                0.3            all           0.25
    ## 3                0.3            all           0.25

Each specific phyloseq subset (treatment-control comparison) is delimited with the “.id” column.

    df_l2fc %>% .$.id %>% unique

    ## [1] "(Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Cel' & Day == '3')"
    ## [2] "(Substrate=='12C-Con' & Day=='14') | (Substrate=='13C-Cel' & Day == '14')"
    ## [3] "(Substrate=='12C-Con' & Day=='14') | (Substrate=='13C-Glu' & Day == '14')"
    ## [4] "(Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Glu' & Day == '3')"

For clarity, let’s edit these long strings to make them more readable when plotted.

    df_l2fc = df_l2fc %>%
      mutate(.id = gsub(' \\| ', '\n', .id))
    df_l2fc %>% .$.id %>% unique

    ## [1] "(Substrate=='12C-Con' & Day=='3')\n(Substrate=='13C-Cel' & Day == '3')"
    ## [2] "(Substrate=='12C-Con' & Day=='14')\n(Substrate=='13C-Cel' & Day == '14')"
    ## [3] "(Substrate=='12C-Con' & Day=='14')\n(Substrate=='13C-Glu' & Day == '14')"
    ## [4] "(Substrate=='12C-Con' & Day=='3')\n(Substrate=='13C-Glu' & Day == '3')"

How many incorporators (rejected hypotheses) & which sparsity cutoff was used for each comparison?

Note: you could set one sparsity cutoff for all comparisons by setting the `sparsity_cutoff` to a specific value.

    df_l2fc %>%
      filter(padj < padj_cutoff) %>%
      group_by(.id, sparsity_threshold) %>%
      summarize(n_incorp_OTUs = OTU %>% unique %>% length) %>%
      as.data.frame

    ##                                                                        .id
    ## 1 (Substrate=='12C-Con' & Day=='14')\n(Substrate=='13C-Cel' & Day == '14')
    ## 2 (Substrate=='12C-Con' & Day=='14')\n(Substrate=='13C-Glu' & Day == '14')
    ## 3   (Substrate=='12C-Con' & Day=='3')\n(Substrate=='13C-Cel' & Day == '3')
    ## 4   (Substrate=='12C-Con' & Day=='3')\n(Substrate=='13C-Glu' & Day == '3')
    ##   sparsity_threshold n_incorp_OTUs
    ## 1               0.15           161
    ## 2               0.30           234
    ## 3               0.30             8
    ## 4               0.30            50

How about a breakdown of incorporators for each phylum in each comparision.

    # summarizing
    df_l2fc_s = df_l2fc %>%
      filter(padj < padj_cutoff) %>%
      mutate(Rank2 = gsub('^__', '', Rank2)) %>%
      group_by(.id, Rank2) %>%
      summarize(n_incorp_OTUs = OTU %>% unique %>% length) %>%
      ungroup()

    # plotting
    ggplot(df_l2fc_s, aes(Rank2, n_incorp_OTUs)) +
        geom_bar(stat='identity') +
        labs(x='Phylum', y='Number of incorporators') +
        facet_wrap(~ .id, scales='free') +
        theme_bw() +
        theme(
          axis.text.x = element_text(angle=55, hjust=1)
        )

# MW-HR-SIP

MW-HR-SIP is run very similarly to HRSIP, but it uses multiple buoyant density (BD) windows. MW-HR-SIP is performed with the `HRSIP()` function, but multiple BD windows are specified.

Let’s use 3 buoyant density windows (g/ml):

> 1.70-1.73, 1.72-1.75, 1.74-1.77

    windows = data.frame(density_min=c(1.70, 1.72, 1.74),
                         density_max=c(1.73, 1.75, 1.77))
    windows

    ##   density_min density_max
    ## 1        1.70        1.73
    ## 2        1.72        1.75
    ## 3        1.74        1.77

Running HRSIP with all 3 BD windows. Let’s run this in parallel to speed things up.

You can turn off parallel processing by setting the `parallel` option to `FALSE`. Also, there’s 2 different levels that could be parallelized: either the `ldply()` or `HRSIP()`. Here, I’m running HRSIP in parallel, but it may make sense in other situations (eg., many comparisons but few density windows and/or sparsity cutoffs) to use ldply in parallel only.

    doParallel::registerDoParallel(ncores)

    df_l2fc = plyr::ldply(physeq_S2D2_l,
                          HRSIP,
                          density_windows = windows,
                          design = ~Substrate,
                          padj_cutoff = padj_cutoff,
                          sparsity_threshold = c(0,0.15,0.3), # just using 3 thresholds to reduce run time
                          .parallel = TRUE)
    df_l2fc %>% head(n=3)

    ##                                                                       .id
    ## 1 (Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Cel' & Day == '3')
    ## 2 (Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Cel' & Day == '3')
    ## 3 (Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Cel' & Day == '3')
    ##        OTU log2FoldChange         p      padj    Rank1            Rank2
    ## 1  OTU.514     -0.1952960 0.7648157 1.0000000 Bacteria __Proteobacteria
    ## 2 OTU.1136      1.5682247 0.0199097 0.8410301 Bacteria  __Acidobacteria
    ## 3 OTU.2119      0.4825207 0.3663796 1.0000000 Bacteria  __Acidobacteria
    ##                   Rank3                  Rank4            Rank5
    ## 1 __Deltaproteobacteria    __Desulfobacterales __Nitrospinaceae
    ## 2               __DA023 __uncultured_bacterium             <NA>
    ## 3               __DA023                   <NA>             <NA>
    ##          Rank6                  Rank7 Rank8 density_min density_max
    ## 1 __uncultured __uncultured_bacterium  <NA>         1.7        1.73
    ## 2         <NA>                   <NA>  <NA>         1.7        1.73
    ## 3         <NA>                   <NA>  <NA>         1.7        1.73
    ##   sparsity_threshold sparsity_apply l2fc_threshold
    ## 1                0.3            all           0.25
    ## 2                0.3            all           0.25
    ## 3                0.3            all           0.25

Let’s check that we have all treatment-control comparisons.

    df_l2fc %>% .$.id %>% unique

    ## [1] "(Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Cel' & Day == '3')"
    ## [2] "(Substrate=='12C-Con' & Day=='14') | (Substrate=='13C-Cel' & Day == '14')"
    ## [3] "(Substrate=='12C-Con' & Day=='14') | (Substrate=='13C-Glu' & Day == '14')"
    ## [4] "(Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Glu' & Day == '3')"

How many incorporators (rejected hypotheses) & which sparsity cutoff was used for each comparison?

Note: one sparsity cutoff could be set for all comparisons by setting the `sparsity_cutoff` to a specific value.

    df_l2fc %>%
      filter(padj < padj_cutoff) %>%
      group_by(.id, sparsity_threshold) %>%
      summarize(n_incorp_OTUs = OTU %>% unique %>% length) %>%
      as.data.frame

    ##                                                                         .id
    ## 1 (Substrate=='12C-Con' & Day=='14') | (Substrate=='13C-Cel' & Day == '14')
    ## 2 (Substrate=='12C-Con' & Day=='14') | (Substrate=='13C-Glu' & Day == '14')
    ## 3   (Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Cel' & Day == '3')
    ## 4   (Substrate=='12C-Con' & Day=='3') | (Substrate=='13C-Glu' & Day == '3')
    ##   sparsity_threshold n_incorp_OTUs
    ## 1               0.15           170
    ## 2               0.30           308
    ## 3               0.30             9
    ## 4               0.30            71

The density windows can vary for each OTU. Let’s plot which density windows were used for the OTUs in the dataset.

    # summarizing
    df_l2fc_s = df_l2fc %>%
      mutate(.id = gsub(' \\| ', '\n', .id)) %>%
      filter(padj < padj_cutoff) %>%
      mutate(density_range = paste(density_min, density_max, sep='-')) %>%
      group_by(.id, sparsity_threshold, density_range) %>%
      summarize(n_incorp_OTUs = OTU %>% unique %>% length)

    #plotting
    ggplot(df_l2fc_s, aes(.id, n_incorp_OTUs, fill=density_range)) +
        geom_bar(stat='identity', position='fill') +
        labs(x='Control-treatment comparision', y='Fraction of incorporators') +
        scale_y_continuous(expand=c(0,0)) +
        theme_bw() +
        theme(
          axis.text.x = element_text(angle=55, hjust=1)
        )

Different BD windows were used for different treatment-control comparisons because the amount of BD shift likely varied among taxa. For example, if a taxon incorporates 100% 13C isotope, then a very ‘heavy’ BD window may show a larger l2fc than a less ‘heavy’ BD window.

Let’s look at a breakdown of incorporators for each phylum in each comparision.

    # summarizing
    df_l2fc_s = df_l2fc %>%
      mutate(.id = gsub(' \\| ', '\n', .id)) %>%
      filter(padj < padj_cutoff) %>%
      mutate(Rank2 = gsub('^__', '', Rank2)) %>%
      group_by(.id, Rank2) %>%
      summarize(n_incorp_OTUs = OTU %>% unique %>% length) %>%
      ungroup()

    # plotting
    ggplot(df_l2fc_s, aes(Rank2, n_incorp_OTUs)) +
        geom_bar(stat='identity') +
        labs(x='Phylum', y='Number of incorporators') +
        facet_wrap(~ .id, scales='free') +
        theme_bw() +
        theme(
          axis.text.x = element_text(angle=55, hjust=1)
        )

Note that the MW-HR-SIP method identifies more incorporators than the HR-SIP method (which uses just one BD window).

MW-HR-SIP detects more taxa for 2 main reasons. First, taxa vary in G+C content, so using only 1 BD window likely encompasses BD shifts for taxa of certain G+C contents (eg., ~50% G+C), but may miss other taxa with higher or lower G+C content. Second, taxa can vary in how much isotope was incorporated, which will affect where each taxon’s DNA is in the density gradient.

# 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