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

library(knitr)
## 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 functions to make the plots from the following manuscript:

Dick JM. 2025. Earth's oxygenation recorded in genomes. Manuscript in preparation.

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

library(JMDplots)

Genome-wide differences of oxidation state between two lineages of methanogens (Figure 1)

genoGOE_1()

Data sources: List of Class I and II methanogen species [@LL18], archaeal marker genes used in GTDB [@RCM+21], and protein abundance in Methanococcus maripaludis [@XHZ+06].

Carbon oxidation state of proteins in eukaryotic gene age groups (Figure 2)

genoGOE_2()

Data sources: Consensus gene ages [@LMM16], protein sequences [@Uni23], and divergence times [@KSC+22].

Evolutionary oxidation of ancestral Rubiscos and thermodynamic prediction of redox boundaries around the GOE (Figure 3)

genoGOE_3()

Data sources: Reconstructed ancestral sequences of Rubisco [@KHAB17] and lists of methanogen [@LL18] and Nitrososphaeria [@RFH+19] genomes.

Evolutionary oxidation and relative stabilities for genomes with S-cycling genes (Figure 4)

genoGOE_4()

Data sources: Genomes with different S-cycling genes [@MCK+23].

Carbon oxidation state of ancestral and extant nitrogenases (Figure 5)

genoGOE_5()

Data source: @GMKK20.

Composition reaction for ancestral Rubisco I/II/III (Reaction 1)

Show code

# Setup basis species and swap in e- for O2
# NOTE: This sets pH to 7; changing it will change the charge (Z) and affinity below
basis("QEC+")
swap.basis("O2", "e-")

# Calculate Eh from pe
pe <- -basis()$logact[rownames(basis()) == "e-"]
Eh <- round(convert(pe, "Eh"), 3)

# Read amino acid compositions and add the first protein to CHNOSZ
fasta_file <- system.file("extdata/fasta/KHAB17.fasta", package = "canprot")
aa <- read_fasta(fasta_file)[1, ]
ip <- add.protein(aa)

# Run affinity() to calculate the affinity of the composition reaction for the first protein
affinity_out <- affinity(iprotein = ip)
affinity_protein <- affinity_out$values[[1]]

# Run protein.basis() to calculate the stoichiometry of the composition reaction
reaction_protein <- protein.basis(ip)

# NOTE: Affinity and reaction stoichiometry are already corrected for pH
# Next we construct the elemental formula of the uncharged protein ...
elements_protein <- protein.formula(ip)
# ... and make the pH correction
Z_protein <- ionize.aa(aa)[[1]]
elements_protein <- cbind(elements_protein, Z = 2)
elements_protein[, "H"] <- elements_protein[, "H"] + Z_protein

# Divide values by protein length and do some rounding
plength <- protein.length(ip)
affinity_residue <- round(affinity_protein / plength, 3)
reaction_residue <- round(reaction_protein / plength, 3)
elements_residue <- round(elements_protein / plength, 3)
# Format the elemental formula
formula_residue <- as.chemical.formula(round(elements_residue, 3))

# Print the results
print(paste("Eh:", Eh, "volt"))
print(paste("Elemental formula of protein per residue:", formula_residue))
print("Stoichiometry of composition reaction:")
print(reaction_residue)
print(paste("Affinity/2.303RT of composition reaction:", affinity_residue))

print(paste("Eh:", Eh, "volt"))
print(paste("Elemental formula of protein per residue:", formula_residue))
print("Stoichiometry of composition reaction:")
print(reaction_residue)
print(paste("Affinity/2.303RT of composition reaction:", affinity_residue))

References

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