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

## some frequently used HTML expressions
# use lowercase here because these tend to be variable names in the examples
zc <- "<i>Z</i><sub>C</sub>"
o2 <- "O<sub>2</sub>"
h2o <- "H<sub>2</sub>O"

options(width = 90)

This vignette runs the code to make selected plots from the following papers:

Dick JM. 2016. Proteomic indicators of oxidation and hydration state in colorectal cancer. PeerJ 4: e2238. doi: 10.7717/peerj.2238

Dick JM. 2017. Chemical composition and the potential for proteomic transformation in cancer, hypoxia, and hyperosmotic stress. PeerJ 5: e3421. doi: 10.7717/peerj.3421

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

library(JMDplots)

Microbial proteins in colorectal cancer (2016 Figure 4)

Stability fields represent the ranges of oxygen fugacity and water activity where a protein with the mean amino acid composition from the labeled microbial species has a higher per-residue affinity (lower Gibbs energy) of formation than the others. Blue and red shading designate microbes relatively enriched in samples from healthy donors and cancer patients, respectively. Plot (E) is a composite figure in which the intensity of shading corresponds to the number of overlapping healthy- or cancer-enriched microbes in the preceding diagrams

# To reproduce the calculations in the paper, use superseded data for [Gly] and [UPBB] 20190206
OldAAfile <- "extdata/OBIGT/OldAA.csv"
add.OBIGT(system.file(OldAAfile, package = "JMDplots"), c("[Gly]", "[UPBB]"))
# Set up graphics device
layout(cbind(matrix(sapply(list(c(1, 2), c(3, 4)), function(x) rep(rep(x, each=3), 3)), nrow=6, byrow=TRUE),
             matrix(rep(c(0, 5, 5, 5, 5, 0), each=4), nrow=6, byrow=TRUE)))
opar <- par(mar=c(3.3, 3.3, 1.5, 1.5), mgp=c(2.1, 0.7, 0), xaxs="i", yaxs="i", las=1, cex=0.9)
# Resolution for plots
res <- 500
# Basis can be "QEC" or "CHNOS"
basis <- "QEC"
# Read bioproject ids, species names
mfile <- system.file("extdata/Dic16/microbes.csv", package="JMDplots")
bugs <- read.csv(mfile, as.is=TRUE)
# Where to keep the locations of healthy zones
healthbugs <- list()
if(basis=="CHNOS") {
  O2 <- c(-85, -65, res)
  H2O <- c(-15, 5, res)
} else {
  O2 <- c(-75, -55, res)
  H2O <- c(-10, 10, res)
}
# The datasets that are considered
datasets <- c("WCQ+12", "ZTV+14", "CTB+14", "FLJ+15")
# colors
blue <- "#4A6FE3"
lightblue <- "#9DA8E2"
neutral <- "#E2E2E2"
lightred <- "#E495A5"
red <- "#D33F6A"
# Labels
logfO2lab <- expression(log~italic("f")*O[2]*"(g)")
logaH2Olab <- expression(log~italic("a")*H[2]*O*"(liq)")
for(i in 1:4) {
  ibug <- bugs$study==datasets[i] & (is.na(bugs$logodds) | abs(bugs$logodds) > 0.15)
  # Get amino acid compositions
  aafile <- system.file("extdata/Dic16/microbial.aa.csv", package="JMDplots")
  microbial.aa <- read.csv(aafile, stringsAsFactors=FALSE)
  aa <- microbial.aa[match(bugs$bioproject[ibug], microbial.aa$organism), ]
  ip <- add.protein(aa)
  # Set up system, calculate relative stabilities
  col <- ifelse(bugs$upcan[ibug], lightred, lightblue)
  basis(basis)
  a <- affinity(O2=O2, H2O=H2O, iprotein=ip, T=37)
  names <- bugs$abbrv[ibug]
  d <- diagram(a, names=names, fill=col, as.residue=TRUE, tplot=FALSE, xlab=logfO2lab, ylab=logaH2Olab, format.names=FALSE)
  if(i==1) title(main="Fecal 16S rRNA", cex.main=1)
  if(i==2) title(main="Fecal metagenome (ZTV+14)", cex.main=1)
  if(i==3) title(main="Co-abundance groups", cex.main=1)
  if(i==4) title(main="Fecal metagenome (FLJ+15)", cex.main=1)
  box()
  label.figure(LETTERS[i], yfrac=0.96, paren=FALSE, font=2, cex=1)
  # Store locations of healthy bug zones
  p <- d$predominant
  p[p %in% which(bugs$upcan[ibug])] <- 0
  p[p != 0] <- 1
  healthbugs[[i]] <- p
}
# Now show the healthy zones
xs <- seq(O2[1], O2[2], length.out=O2[3])
ys <- seq(H2O[1], H2O[2], length.out=H2O[3])
hhh <- healthbugs[[1]] + healthbugs[[2]] + healthbugs[[3]] + healthbugs[[4]]
image(xs, ys, hhh, col=c(red, lightred, neutral, lightblue, blue), useRaster=TRUE, xlab=logfO2lab, ylab=logaH2Olab)
title(main="Cumulative stability count")
box()
label.figure("E", yfrac=0.96, paren=FALSE, font=2, cex=1)
# Reset thermodynamic database
reset()

Data sources: A. @WCQ+12. B. @ZTV+14. C. @CTB+14. D. @FLJ+15.

Potential diagrams: Pancreatic cancer (2017 Figure S3 and 3E)

The potential diagrams show the weighted rank difference of chemical affinities between up- and down-expressed proteins in each dataset. Groups of datasets are considered that have similar chemical features, i.e. changes in r zc and nr h2o.

Here we make plots for datasets for pancreatic cancer having a mean difference of nr h2o that is > 0.01 and a small r zc, as judged by the p-value and common language effect size (CLES). Red and blue correspond to greater potential for formation of the up- and down-expressed proteins, respectively; the line of equipotential is shown in white:

gpresult <- groupplots("pancreatic_H2O_up", res = 25)

Now let's make a merged diagram. The red-white-blue shading is computed from the mean of the previous diagrams. The black lines show the median and quartiles for the y-positions of the equipotential lines in the previous diagrams. The second plot shows effective values of Eh (redox potential) as a function of the same variables (oxygen fugacity and water activity) (see 2016 Figure 6I).

par(mfrow = c(1, 2))
mergedplot(gpresult, res = 25)
Ehplot(xlim = c(-70, -62), ylim = c(-6, 2), dy = 0.1)

Data sources: @LHE+04, @MLC+11, @KHO+13, @KPC+13, @PKB+13, @KKC+16.

Basis species comparison (2017 Figure S1)

These plots show projections of elemental composition of proteins made using two sets of basis species. Using the CHNOS basis species (CO₂, NH₃, H₂S, r h2o, r o2), the plots show that nr h2o and nr o2, i.e. the number of r h2o and r o2 in the formation per residue of the proteins from basis species, are both moderately correlated with r zc (average oxidation state of carbon). Using the QEC basis species (glutamine, glutamic acid, cysteine, r h2o, r o2), we find that nr o2 is strongly correlated with r zc, but nr h2o shows very little correlation. Accordingly, the QEC basis more clearly exposes two chemical variables -- oxidation state and hydration state -- in proteomic data.

Here we define some labels used in the plot.

QEClab <- CHNOSZ::syslab(c("glutamine", "glutamic acid", "cysteine", "H2O", "O2"))
CHNOSlab <- CHNOSZ::syslab(c("CO2", "NH3", "H2S", "H2O", "O2"))

Next, get the amino acid compositions of all proteins in the UniProt human proteome and calculate the protein formulas and r zc. Note that r zc is a sum of elemental ratios and is independent of the choice of basis species.

aa <- get("human.base", canprot)
protein.formula <- CHNOSZ::protein.formula(aa)
ZC <- CHNOSZ::ZC(protein.formula)

Now set up the figure and plot the per-residue elemental compositions of the proteins projected into different sets of basis species.

par(mfrow = c(2, 2))
par(mar = c(4, 4, 2.5, 1))
par(cex = 1.1)
par(mgp = c(2.5, 1, 0))
for(basis in c("QEC", "CHNOS")) {
  CHNOSZ::basis(basis)
  protein.basis <- CHNOSZ::protein.basis(aa)
  protein.length <- CHNOSZ::protein.length(aa)
  residue.basis <- protein.basis / protein.length
  smoothScatter(ZC, residue.basis[, "O2"], xlab = cplab$Zc, ylab = cplab$nO2)
  smoothScatter(ZC, residue.basis[, "H2O"], xlab = cplab$Zc, ylab = cplab$nH2O)
  if(basis=="QEC") mtext(QEClab, outer = TRUE, cex = 1.2, line = -1.5)
  if(basis=="CHNOS") mtext(CHNOSlab, outer = TRUE, cex = 1.2, line = -15)
}