orp16S.Rmd

---
title: "Influence of redox potential on bacterial protein evolution (2023)"
output: 
  rmarkdown::html_vignette:
    mathjax: null
vignette: >
  %\VignetteIndexEntry{Influence of redox potential on bacterial protein evolution (2023)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
bibliography: [orp16S.bib, JMDplots.bib]
csl: elementa.csl
link-citations: true
---

<style>
/* https://gomakethings.com/how-to-break-an-image-out-of-its-parent-container-with-css/ */
@media (min-width: 700px) {
  .full-width {
    left: 50%;
    margin-left: -50vw;
    margin-right: -50vw;
    max-width: 100vw;
    position: relative;
    right: 50%;
    width: 100vw;
  }
}
@media (min-width: 900px) {
  .full-width {
    left: 50vw; /* fallback if needed */
    left: calc(50vw - 100px);
    width: 900px;
    position: relative;
  }
}
</style>

```{r setup, include=FALSE}
library(knitr)
library(kableExtra)
## use pngquant to reduce size of PNG images
knit_hooks$set(pngquant = hook_pngquant)
pngquant <- "--speed=1 --quality=0-25"
# in case pngquant isn't available
if (!nzchar(Sys.which("pngquant"))) pngquant <- NULL 

## colorize messages 20171031
## adapted from https://gist.github.com/yihui/2629886#file-knitr-color-msg-rnw
color_block = function(color) {
  function(x, options) sprintf('<pre style="color:%s">%s</pre>', color, x)
}
knit_hooks$set(warning = color_block('magenta'), error = color_block('red'), message = color_block('blue'))
```

```{r options, include=FALSE}
options(width = 80)
```

```{r HTML, include=FALSE}
Zc <- "<i>Z</i><sub>C</sub>"
```

This vignette runs the code to make the plots from the following paper:

> Dick JM, Meng D. 2023. Community- and genome-based evidence for a shaping influence of redox potential on bacterial protein evolution. *mSystems* **8**(3): e00014-23. doi: [10.1128/msystems.00014-23](https://doi.org/10.1128/msystems.00014-23)

**Post-publication update** on 2024-07-10: The reference database for constructing community reference proteomes has been changed from RefSeq 206 to GTDB 220.
Correspondingly, taxonomic classifications now use [RDP Classifier training files for 16S rRNA sequences from GTDB](https://doi.org/10.5281/zenodo.12703477) instead of the default RDP Classifier training set.
This update affects the numerical values in all the figures (except Fig. 4) but does not change the basic trends.

This vignette was compiled on `r Sys.Date()` with **[JMDplots](https://github.com/jedick/JMDplots)** `r packageDescription("JMDplots")$Version` and **[chem16S](https://github.com/jedick/chem16S)** `r packageDescription("chem16S")$Version`.

```{r library, message = FALSE, results = "hide"}
library(JMDplots)
```

## Thermodynamic model for the relationship between carbon oxidation state of reference proteomes and redox potential (Figure 1)
```{r orp16S_1, message = FALSE, results = "hide", out.width='100%', fig.align='center', fig.width = 7, fig.height = 6, pngquant = pngquant, dpi = 100}
orp16S_1()
```

## `r Zc` of reference proteomes compared with oxygen tolerance and with metaproteomes (Figure 2)
```{r orp16S_2, message = FALSE, out.width='100%', fig.width = 8, fig.height = 3.5, out.extra = 'class="full-width"', pngquant = pngquant, dpi = 100}
orp16S_2()
```

**Data sources:** List of strictly anaerobic and aerotolerant genera [@MR18], Metaproteomes: Manus Basin Inactive Chimney [@MPB+19], Manus Basin Active Chimneys [@PMM+18], Soda Lake Biomats [@KTS+17], Mock Communities [@KTS+17], Saanich Inlet [@HTZ+17].

## Methods overview and chemical depth profiles in Winogradsky columns (Figure 3)
```{r orp16S_3, message = FALSE, results = "hide", out.width='100%', fig.width = 11.25, fig.height = 4, out.extra = 'class="full-width"', pngquant = pngquant, dpi = 100}
orp16S_3()
```

**Data sources:** **(b)** ORP and pH [@DFYS19] to calculate Eh7. **(c)** 16S rRNA gene sequences [@RBW+14] to calculate `r Zc`.

## Sample locations on world map and Eh-pH diagram (Figure 4)
```{r orp16S_4, message = FALSE, results = "hide", out.width='100%', fig.width = 26, fig.height = 23, out.extra = 'class="full-width"', pngquant = pngquant, dpi = 100}
orp16S_4()
```

**Data sources:** **(a)** coastlineWorld dataset in R package oce [@KR22]; shapefiles for North American Great Lakes [@USGS10]. **(b)** Outline is from @BKM60.

## Associations between Eh7 and `r Zc` at local scales (Figure 5)
```{r orp16S_5, message = FALSE, results = "hide", out.width='100%', fig.width = 8, fig.height = 6, pngquant = pngquant, dpi = 100}
orp16S_5()
```

**Data sources:** **(a)** Daya Bay [@WHLH21a], Bay of Biscay [@LMBA21], Hunan Province [@MLL+19]. **(c)** **1** acidic and **2** circumneutral to alkaline New Zealand hot springs [@PCL+18], **3** Eastern Tibetan Plateau [@GWS+20], **4** Uzon Caldera [@PBU+20], **5** Southern Tibetan Plateau [@MWY+21].

## Tally of regression slopes and binomial *p*-values (Table 1a)
```{r orp16S_T1a}
colnames <- c("Ntot", "Npos", "Nneg", "p")
kable(orp16S_T1(), col.names = rep(colnames, 2)) %>%
  add_header_above(c(" " = 1, "Bacteria" = 4, "Archaea" = 4))
```

## Tally of regression slopes with same sign of minimum and maximum values in 95% confidence interval (Table 1b)
```{r orp16S_T1b}
colnames <- c("Ntot", "Npos", "Nneg", "p")
kable(orp16S_T1(samesign = TRUE), col.names = rep(colnames, 2)) %>%
  add_header_above(c(" " = 1, "Bacteria" = 4, "Archaea" = 4))
```

## Associations between Eh7 and `r Zc` at a global scale (Figure 6)
```{r orp16S_6, message = FALSE, results = "hide", out.width='100%', fig.width = 10, fig.height = 7, out.extra = 'class="full-width"', pngquant = pngquant, dpi = 100}
global.slopes <- orp16S_6()
```

## Comparison of Eh7, Eh, and O~2~ as predictors of carbon oxidation state (Figure S2)
```{r orp16S_S2, message = FALSE, results = "hide", out.width='100%', fig.width = 6, fig.height = 8, pngquant = pngquant, dpi = 100}
orp16S_S2()
```

## Global analysis including only datasets using 515F/806R primers from the Earth Microbiome Project (EMP) (Figure S3)
```{r orp16S_S3, message = FALSE, results = "hide", out.width='100%', fig.width = 10, fig.height = 7, out.extra = 'class="full-width"', pngquant = pngquant, dpi = 100}
orp16S_6(EMP_primers = TRUE)
```

## Genus-level taxonomic classifications (using 16S rRNA sequences from GTDB) mapped to the GTDB taxonomy
**NOTE**: With the switch to GTDB on both ends, this shows 100% mapping rate, higher than the values mentioned in the article for classifications from the RDP training set mapped to the NCBI taxonomy.
```{r mapperc}
# Read Dataset S3 (created by orp16S_D3())
dat <- read.csv(system.file("extdata/orp16S/Dataset_S3.csv", package = "JMDplots"))
perc_bac <- round(mean(subset(dat, lineage == "Bacteria")$mapperc))
perc_arc <- round(mean(subset(dat, lineage == "Archaea")$mapperc))
message(paste0(perc_bac, "% for bacteria, ", perc_arc, "% for archaea"))
```

## References