In jedick/JMDplots: Plots from Papers by Jeffrey M. Dick

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'))

options(width = 80)

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

Post-publication update on 2024-07-10: The reference database for constructing community reference proteomes has been changed from RefSeq 206 to GTDB 220. Likewise, taxonomic classifications now use RDP Classifier training files for 16S rRNA sequences from GTDB instead of the default RDP Classifier training set. This update affects the results shown in the plots (except for Fig. 4) but in most cases does not change the basic trends. The slope in Fig. 2b is now negative instead of positive, but the 95% confidence interval still contains zero. The number of samples with classifications for archaeal communities after filtering increased by more than 100, and the slope for Archaea in Fig. 6b is now significantly positive instead of non-significantly negative.

This vignette was compiled on r Sys.Date() with JMDplots r packageDescription("JMDplots")$Version and chem16S r packageDescription("chem16S")$Version.

library(JMDplots)

Thermodynamic model for the relationship between carbon oxidation state of reference proteomes and redox potential (Figure 1)

orp16S_1()

r Zc of reference proteomes compared with oxygen tolerance and with metaproteomes (Figure 2)

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)

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)

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)

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)

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)

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)

global.slopes <- orp16S_6()

Comparison of Eh7, Eh, and O~2~ as predictors of carbon oxidation state (Figure S2)

orp16S_S2()

Global analysis including only datasets using 515F/806R primers from the Earth Microbiome Project (EMP) (Figure S3)

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.

# 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

jedick/JMDplots documentation built on April 12, 2025, 1:35 p.m.