[CHNOSZ-commits] r775 - in pkg/CHNOSZ: . R demo man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Mar 4 04:43:55 CET 2023
Author: jedick
Date: 2023-03-04 04:43:54 +0100 (Sat, 04 Mar 2023)
New Revision: 775
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/demo/NaCl.R
pkg/CHNOSZ/demo/ORP.R
pkg/CHNOSZ/demo/Shh.R
pkg/CHNOSZ/demo/TCA.R
pkg/CHNOSZ/demo/affinity.R
pkg/CHNOSZ/demo/aluminum.R
pkg/CHNOSZ/demo/arsenic.R
pkg/CHNOSZ/demo/buffer.R
pkg/CHNOSZ/demo/carboxylase.R
pkg/CHNOSZ/demo/dehydration.R
pkg/CHNOSZ/demo/density.R
pkg/CHNOSZ/demo/glycinate.R
pkg/CHNOSZ/demo/ionize.R
pkg/CHNOSZ/demo/lambda.R
pkg/CHNOSZ/demo/potassium.R
pkg/CHNOSZ/demo/protbuff.R
pkg/CHNOSZ/demo/protein.equil.R
pkg/CHNOSZ/demo/saturation.R
pkg/CHNOSZ/demo/sources.R
pkg/CHNOSZ/man/examples.Rd
Log:
Capitalize comments in demos
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2023-03-04 02:41:46 UTC (rev 774)
+++ pkg/CHNOSZ/DESCRIPTION 2023-03-04 03:43:54 UTC (rev 775)
@@ -1,6 +1,6 @@
Date: 2023-03-04
Package: CHNOSZ
-Version: 1.9.9-66
+Version: 1.9.9-67
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2023-03-04 02:41:46 UTC (rev 774)
+++ pkg/CHNOSZ/R/examples.R 2023-03-04 03:43:54 UTC (rev 775)
@@ -1,8 +1,8 @@
# CHNOSZ/examples.R
# Functions to run all examples and demos in the package
-examples <- function(save.png=FALSE) {
- # Run all the examples in CHNOSZ documentation
+examples <- function(save.png = FALSE) {
+ # Run examples in each of the CHNOSZ help pages
.ptime <- proc.time()
topics <- c("thermo", "examples",
"util.array", "util.data", "util.expression", "util.legend", "util.plot",
@@ -13,15 +13,17 @@
"solubility", "equilibrate",
"diagram", "mosaic", "mix",
"buffer", "nonideal", "NaCl",
- "add.protein", "ionize.aa")
+ "add.protein", "ionize.aa", "rank.affinity",
+ "DEW", "logB.to.OBIGT", "stack_mosaic"
+ )
plot.it <- FALSE
if(is.character(save.png))
- png(paste(save.png,"%d.png",sep=""),width=500,height=500,pointsize=12)
+ png(paste(save.png, "%d.png", sep = ""), width = 500, height = 500, pointsize = 12)
else if(save.png) plot.it <- TRUE
for(i in 1:length(topics)) {
- if(plot.it) png(paste(topics[i],"%d.png",sep=""),width=500,height=500,pointsize=12)
- myargs <- list(topic=topics[i],ask=FALSE)
- do.call(example,myargs)
+ if(plot.it) png(paste(topics[i], "%d.png", sep = ""), width = 500, height = 500, pointsize = 12)
+ myargs <- list(topic = topics[i], ask = FALSE)
+ do.call(example, myargs)
if(plot.it) dev.off()
}
if(is.character(save.png)) dev.off()
@@ -28,12 +30,13 @@
cat("Time elapsed: ", proc.time() - .ptime, "\n")
}
-demos <- function(which=c("sources", "protein.equil", "affinity", "NaCl", "density",
+demos <- function(which = c("sources", "protein.equil", "affinity", "NaCl", "density",
"ORP", "ionize", "buffer", "protbuff", "glycinate",
"mosaic", "copper", "arsenic", "solubility", "gold", "contour", "sphalerite", "minsol",
"Shh", "saturation", "DEW", "lambda", "potassium", "TCA", "aluminum",
- "AD", "comproportionation", "Pourbaix", "E_coli"), save.png=FALSE) {
- # Run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
+ "AD", "comproportionation", "Pourbaix", "E_coli", "yttrium", "rank.affinity"), save.png = FALSE) {
+ # Run one or more demos from CHNOSZ with ask = FALSE, and return the value of the last one
+ out <- NULL
for(i in 1:length(which)) {
# A message so the user knows where we are
message("------------")
@@ -40,15 +43,15 @@
if(which[i]=="dehydration" & !save.png) {
message("demos: skipping dehydration demo as save.png is FALSE")
next
- } else message(paste("demos: running '", which[i], "'", sep=""))
- if(save.png & !which[i]=="dehydration") {
+ } else message(paste("demos: running '", which[i], "'", sep = ""))
+ if(save.png & !which[i] == "dehydration") {
width <- 500
height <- 500
- if(which[i]=="comproportionation") width <- 600
- png(paste(which[i], "%d.png", sep=""), width = width, height = height, pointsize = 12)
+ if(which[i] == "comproportionation") width <- 600
+ png(paste(which[i], "%d.png", sep = ""), width = width, height = height, pointsize = 12)
}
- out <- demo(which[i], package="CHNOSZ", character.only=TRUE, echo=FALSE, ask=FALSE)
- if(save.png & !which[i]=="dehydration") dev.off()
+ out <- demo(which[i], package = "CHNOSZ", character.only = TRUE, echo = FALSE, ask = FALSE)
+ if(save.png & !which[i] == "dehydration") dev.off()
}
return(invisible(out))
}
Modified: pkg/CHNOSZ/demo/NaCl.R
===================================================================
--- pkg/CHNOSZ/demo/NaCl.R 2023-03-04 02:41:46 UTC (rev 774)
+++ pkg/CHNOSZ/demo/NaCl.R 2023-03-04 03:43:54 UTC (rev 775)
@@ -6,11 +6,11 @@
## J. Chem. Soc. Faraday Trans. 88, 803-826. https://doi.org/10.1039/FT9928800803 )
library(CHNOSZ)
-## uncomment these lines to make the plot with the g-function disabled
+## Uncomment these lines to make the plot with the g-function disabled
#mod.OBIGT("Cl-", z=0)
#mod.OBIGT("Na+", z=0)
-# start a new plot and show the experimental logK
+# Start a new plot and show the experimental logK
thermo.plot.new(xlim=c(0, 1000), ylim=c(-5.5, 1),
xlab=axis.label("T"), ylab=axis.label("logK"))
expt <- read.csv(system.file("extdata/cpetc/SOJSH.csv",
@@ -17,11 +17,11 @@
package="CHNOSZ"), as.is=TRUE)
points(expt$T,expt$logK, pch=expt$pch)
-# we'll be at 9 distinct pressure conditions, including Psat
+# We'll be at 9 distinct pressure conditions, including Psat
# Psat is repeated to show "not considered" region
# (T >= 355 degC; Fig. 6 of Shock et al., 1992)
P <- c(list("Psat", "Psat"), as.list(seq(500, 4000, by=500)))
-# for each of those what's the range of temperature
+# For each of those what's the range of temperature
T <- list()
T[[1]] <- seq(0, 354, 1)
T[[2]] <- seq(354, 370, 1)
@@ -30,7 +30,7 @@
T[[5]] <- seq(395, 920, 1)
T[[6]] <- T[[7]] <- T[[8]] <- T[[9]] <- T[[10]] <- seq(400, 1000, 1)
-# calculate and plot the logK
+# Calculate and plot the logK
species <- c("NaCl", "Na+", "Cl-")
coeffs <- c(-1, 1, 1)
logK <- numeric()
@@ -44,7 +44,7 @@
if(i > 2) logK <- c(logK, splinefun(s$out$T, s$out$logK)(Texpt))
}
-# add title, labels, and legends
+# Add title, labels, and legends
title(describe.reaction(s$reaction, states = 1))
text(150, -0.1, quote(italic(P)[sat]), cex=1.2)
text(462, -4, "500 bar")
@@ -58,6 +58,6 @@
l2 <- "Non-recommended region (Shock et al., 1992, Fig. 6)"
legend("topright", as.expression(c(l1, l2)), lty=c(1, 3), bty="n")
-# test for average divergence (excluding Psat)
+# Test for average divergence (excluding Psat)
expt <- expt[!expt$P %in% "Psat", ]
stopifnot(mean(abs(logK - expt$logK)) < 0.09)
Modified: pkg/CHNOSZ/demo/ORP.R
===================================================================
--- pkg/CHNOSZ/demo/ORP.R 2023-03-04 02:41:46 UTC (rev 774)
+++ pkg/CHNOSZ/demo/ORP.R 2023-03-04 03:43:54 UTC (rev 775)
@@ -1,175 +1,174 @@
# CHNOSZ/demo/ORP.R
-# first version 20100715 jmd
+# First version 20100715 jmd
library(CHNOSZ)
-# calculate the temperature dependence of
+# Calculate the temperature dependence of
# potentials vs. standard hydrogen electrode (SHE) of various electrodes (Ag/AgCl)
# and ORP standards (ZoBell, Light's, (tri)iodide)
-# use the extended Debye-Huckel equation with Alberty's parameters
+# Use the extended Debye-Huckel equation with Alberty's parameters
oldnon <- nonideal("Alberty")
# Bard et al.'s fit to the potential
# (Bard, Parson, Jordan, Standard Potentials In Aqueous Solution, 1985)
-AgAgCl.Bard <- function(T,high.T=TRUE) {
+AgAgCl.Bard <- function(T, high.T = TRUE) {
# we use the corrected high-T formula from wikipedia
if(high.T) return(0.23737 - 5.3783e-4 * T - 2.3728e-6 * T^2 - 2.2671e-9 * (T+273))
else return(0.23695 - 4.8564e-4 * T - 3.4205e-6 * T^2 - 5.869e-9 * (T+273))
}
-# function to calculate the potential of Ag/AgCl vs. SHE
+# Function to calculate the potential of Ag/AgCl vs. SHE
# Ag(s) + Cl- = AgCl(s) + e-
# logK = -pe - logaCl
# pe = -logK - logmCl - loggamCl
# ORP = RT/F * (logK - logmCl - loggamCl)
-AgAgCl <- function(T,mKCl=4) {
+AgAgCl <- function(T, mKCl = 4) {
# mKCl is the molality of KCl in the electrolyte
- # we take it as a first approximation to be equal to
+ # We take it as a first approximation to be equal to
# the molality of Cl- (and to the ionic strength)
logmCl <- log10(mKCl)
- # get the logK for the reaction
- logK <- subcrt(c("Ag","Cl-","AgCl","e-"),c(-1,-1,1,1),c("cr","aq","cr","aq"),T=T)$out$logK
- # get the activity coefficient for Cl-
- loggamCl <- log(10) * subcrt("Cl-",T=T,IS=mKCl)$out[[1]]$loggam
- # get the pe for the solution
+ # Get the logK for the reaction
+ logK <- subcrt(c("Ag", "Cl-", "AgCl", "e-"), c(-1, -1, 1, 1), c("cr", "aq", "cr", "aq"), T = T)$out$logK
+ # Get the activity coefficient for Cl-
+ loggamCl <- log(10) * subcrt("Cl-", T = T, IS = mKCl)$out[[1]]$loggam
+ # Get the pe for the solution
pe <- -logK - logmCl - loggamCl
- # convert that to Eh
- Eh <- convert(pe,"Eh",T=convert(T,"K"))
+ # Convert that to Eh
+ Eh <- convert(pe, "Eh", T = convert(T, "K"))
return(Eh)
}
ZoBell <- function(T) {
- # doesn't work very well because we ignore the
+ # Doesn't work very well because we ignore the
# ferricyanide and ferrocyanide complexes
# Fe+2 = Fe+3 + e-
# logK = logaFe3 - logaFe2 - pe
- # get the logK for the reaction
- logK <- subcrt(c("Fe+2","Fe+3","e-"),c(-1,1,1),T=T)$out$logK
- # we use the recipe from standard methods (table 2580:II)
+ # Get the logK for the reaction
+ logK <- subcrt(c("Fe+2", "Fe+3", "e-"), c(-1, 1, 1), T = T)$out$logK
+ # We use the recipe from standard methods (table 2580:II)
# 1.4080 g K4Fe(CN)6.3H2O -> 0.0033333 mol Fe+2
# 1.0975 g K3Fe(CN)6 -> 0.0033333 mol Fe+3
# 7.4555 g KCl -> 0.1 mol Cl-
logmFe2 <- logmFe3 <- log10(0.0033333)
- # get the loggam for the iron species
- loggamFe2 <- log(10) * subcrt("Fe+2",T=T,IS=1)$out[[1]]$loggam
- loggamFe3 <- log(10) * subcrt("Fe+3",T=T,IS=1)$out[[1]]$loggam
- # get the pe for the solution
+ # Get the loggam for the iron species
+ loggamFe2 <- log(10) * subcrt("Fe+2", T = T, IS = 1)$out[[1]]$loggam
+ loggamFe3 <- log(10) * subcrt("Fe+3", T = T, IS = 1)$out[[1]]$loggam
+ # Get the pe for the solution
pe <- -logK + logmFe3 + loggamFe3 - logmFe2 - loggamFe2
- # convert to Eh
- Eh <- convert(pe,"Eh",T=convert(T,"K"))
+ # Convert to Eh
+ Eh <- convert(pe, "Eh", T = convert(T, "K"))
return(Eh)
}
-ZoBell.table <- function(T=NULL,which=NULL) {
- # oxidation-reduction potential of ZoBell's solution
+ZoBell.table <- function(T = NULL, which = NULL) {
+ # Oxidation-reduction potential of ZoBell's solution
# from Standard Methods for Water and Wastewater or YSI
# (interpolated and/or extrapolated as necessary)
- # standard methods (1997) table 2580:I
+ # Standard methods (1997) table 2580:I
Eh.T.SMW <- 1:30
- Eh.SMW <- c(0.481,0.479,0.476,0.474,0.472,0.47,0.468,0.465,0.463,0.461,
- 0.459,0.457,0.454,0.452,0.45,0.448,0.446,0.443,0.441,0.439,0.437,
- 0.435,0.432,0.43,0.428,0.426,0.424,0.421,0.419,0.417)
- # from YSI (2005):
+ Eh.SMW <- c(0.481, 0.479, 0.476, 0.474, 0.472, 0.47, 0.468, 0.465, 0.463, 0.461,
+ 0.459, 0.457, 0.454, 0.452, 0.45, 0.448, 0.446, 0.443, 0.441, 0.439, 0.437,
+ 0.435, 0.432, 0.43, 0.428, 0.426, 0.424, 0.421, 0.419, 0.417)
+ # From YSI (2005):
# Measuring ORP on YSI 6-Series Sondes: Tips, Cautions and Limitations
# NOTE: these values are vs. Ag/AgCl (4 M KCl)
- Eh.T.YSI <- seq(-5,50,by=5)
- Eh.YSI <- c(267.0,260.5,254.0,247.5,241.0,234.5,228.0,221.5,215.0,208.5,202.0,195.5)/1000
- # spline function for each of the tables
- SMW <- splinefun(Eh.T.SMW,Eh.SMW)
- YSI <- splinefun(Eh.T.YSI,Eh.YSI)
- # just one of the tables
+ Eh.T.YSI <- seq(-5, 50, by = 5)
+ Eh.YSI <- c(267.0, 260.5, 254.0, 247.5, 241.0, 234.5, 228.0, 221.5, 215.0, 208.5, 202.0, 195.5)/1000
+ # Spline function for each of the tables
+ SMW <- splinefun(Eh.T.SMW, Eh.SMW)
+ YSI <- splinefun(Eh.T.YSI, Eh.YSI)
+ # Just one of the tables
Eh.fun <- get(which)
- Eh.T <- get(paste("Eh.T",which,sep="."))
+ Eh.T <- get(paste("Eh.T", which, sep = "."))
if(is.null(T)) T <- Eh.T
- return(data.frame(T=T,Eh=Eh.fun(T)))
+ return(data.frame(T = T, Eh = Eh.fun(T)))
}
Light <- function(T) {
- # this is going to look something like
+ # This is going to look something like
# Fe+2 = Fe+3 + e-
# logK = logaFe3 - logaFe2 - pe
- # get the logK for the reaction
- logK <- subcrt(c("Fe+2","Fe+3","e-"),c(-1,1,1),T=T)$out$logK
- # we use the recipe from standard methods (table 2580:II)
+ # Get the logK for the reaction
+ logK <- subcrt(c("Fe+2", "Fe+3", "e-"), c(-1, 1, 1), T = T)$out$logK
+ # We use the recipe from standard methods (table 2580:II)
# 39.21 g Fe(NH4)2(SO4)2(H2O)6 -> 0.1 mol Fe+2
# 48.22 g Fe(NH4)(SO4)2(H2O)12 -> 0.1 mol Fe+3
logmFe2 <- logmFe3 <- log10(0.1)
- # get the loggam for the iron species
- loggamFe2 <- log(10) * subcrt("Fe+2",T=T,IS=0.2)$out[[1]]$loggam
- loggamFe3 <- log(10) * subcrt("Fe+3",T=T,IS=0.2)$out[[1]]$loggam
- # get the pe for the solution
+ # Get the loggam for the iron species
+ loggamFe2 <- log(10) * subcrt("Fe+2", T = T, IS = 0.2)$out[[1]]$loggam
+ loggamFe3 <- log(10) * subcrt("Fe+3", T = T, IS = 0.2)$out[[1]]$loggam
+ # Get the pe for the solution
pe <- -logK + logmFe3 + loggamFe3 - logmFe2 - loggamFe2
- # convert to Eh
- Eh <- convert(pe,"Eh",T=convert(T,"K"))
+ # Convert to Eh
+ Eh <- convert(pe, "Eh", T = convert(T, "K"))
return(Eh)
}
Iodide.table <- function(T=NULL) {
- # oxidation-reduction potential of Thermo's iodide solution
+ # Oxidation-reduction potential of Thermo's iodide solution
# from thermo instruction sheet 255218-001 (articlesFile_18739)
- T.Iodide <- seq(0,50,5)
- Eh.Iodide <- c(438,435,431,428,424,420,415,411,406,401,396)/1000
- Iodide <- splinefun(T.Iodide,Eh.Iodide)
+ T.Iodide <- seq(0, 50, 5)
+ Eh.Iodide <- c(438, 435, 431, 428, 424, 420, 415, 411, 406, 401, 396)/1000
+ Iodide <- splinefun(T.Iodide, Eh.Iodide)
if(is.null(T)) T <- T.Iodide
- return(data.frame(T=T,Eh=Iodide(T)))
+ return(data.frame(T = T, Eh = Iodide(T)))
}
Iodide <- function(T) {
- # this is going to look something like
+ # This is going to look something like
# 3I- = I3- + 2e-
# logK = -2pe + logaI3 - 3logaI
- # get the logK for the reaction
- logK <- subcrt(c("I-","I3-","e-"),c(-3,1,2),T=T)$out$logK
- # could the activities be 0.1 M ... or something else?
+ # Get the logK for the reaction
+ logK <- subcrt(c("I-", "I3-", "e-"), c(-3, 1, 2), T = T)$out$logK
+ # Could the activities be 0.1 M ... or something else?
logmI <- log10(2)
logmI3 <- log10(0.01)
- # get the loggam for the iodine species
- loggamI <- log(10) * subcrt("I-",T=T,IS=0.2)$out[[1]]$loggam
- loggamI3 <- log(10) * subcrt("I3-",T=T,IS=0.2)$out[[1]]$loggam
- # get the pe for the solution
+ # Get the loggam for the iodine species
+ loggamI <- log(10) * subcrt("I-", T = T, IS = 0.2)$out[[1]]$loggam
+ loggamI3 <- log(10) * subcrt("I3-", T = T, IS = 0.2)$out[[1]]$loggam
+ # Get the pe for the solution
pe <- ( -logK + logmI3 + loggamI3 - 3 * (logmI - loggamI) ) / 2
- # convert to Eh
- Eh <- convert(pe,"Eh",T=convert(T,"K"))
+ # Convert to Eh
+ Eh <- convert(pe, "Eh", T = convert(T, "K"))
return(Eh)
}
figure <- function() {
- # make some figures
- # the temperatures we're interested in
- # in degrees C
- T <- seq(0,100,5)
- # temperature-Eh diagram for various electrodes
- thermo.plot.new(ylim=c(0,0.8),xlim=c(0,100),
- ylab=axis.label("Eh"),xlab=axis.label("T"))
- # the Ag/AgCl electrode (Bard et al. fit)
- points(T,AgAgCl.Bard(T),pch=0)
- # the Ag/AgCl electrode (equilibrium calculations)
- lines(T,AgAgCl(T))
+ # Make some figures
+ # Temperatures in degrees C
+ T <- seq(0, 100, 5)
+ # Temperature-Eh diagram for various electrodes
+ thermo.plot.new(ylim = c(0, 0.8), xlim = c(0, 100),
+ ylab = axis.label("Eh"), xlab = axis.label("T"))
+ # Ag/AgCl electrode (Bard et al. fit)
+ points(T, AgAgCl.Bard(T), pch = 0)
+ # Ag/AgCl electrode (equilibrium calculations)
+ lines(T, AgAgCl(T))
# ZoBell's solution (SMW table 2580)
- SMW <- ZoBell.table(which="SMW")
- points(SMW$T,SMW$Eh,pch=1)
+ SMW <- ZoBell.table(which = "SMW")
+ points(SMW$T, SMW$Eh, pch = 1)
# ZoBell's solution (YSI tech report table)
- YSI <- ZoBell.table(which="YSI")
- # make these values referenced to SHE instead of Ag/AgCl
+ YSI <- ZoBell.table(which = "YSI")
+ # Make these values referenced to SHE instead of Ag/AgCl
Eh.YSI <- YSI$Eh + AgAgCl(YSI$T)
- points(YSI$T,Eh.YSI,pch=2)
+ points(YSI$T, Eh.YSI, pch = 2)
# Light's solution (equilibrium values)
- lines(T,Light(T))
+ lines(T, Light(T))
# Light's solution (at 25 degrees only)
- points(25,0.475 + 0.200,pch=3)
+ points(25, 0.475 + 0.200, pch = 3)
# Thermo's I-/I3- solution
Thermo <- Iodide.table()
- points(Thermo$T,Thermo$Eh,pch=4)
- # calculated I-/I3- values
- lines(T,Iodide(T))
- # add some labels
- text(c(30,30,30,50),c(0.72,0.5,0.35,0.25),
- c("Light","ZoBell","(Tri)Iodide","Ag/AgCl"))
- title(main="Potentials vs standard hydrogen electrode (SHE)")
+ points(Thermo$T, Thermo$Eh, pch = 4)
+ # Calculated I-/I3- values
+ lines(T, Iodide(T))
+ # Add some labels
+ text(c(30, 30, 30, 50), c(0.72, 0.5, 0.35, 0.25),
+ c("Light", "ZoBell", "(Tri)Iodide", "Ag/AgCl"))
+ title(main = "Potentials vs standard hydrogen electrode (SHE)")
}
-# finally, make the plot
+# Finally, make the plot
figure()
-# reset the nonideality method to default
+# Reset the nonideality method to default
nonideal(oldnon)
Modified: pkg/CHNOSZ/demo/Shh.R
===================================================================
--- pkg/CHNOSZ/demo/Shh.R 2023-03-04 02:41:46 UTC (rev 774)
+++ pkg/CHNOSZ/demo/Shh.R 2023-03-04 03:43:54 UTC (rev 775)
@@ -2,7 +2,7 @@
# (Dick, 2015. Chemical integration of proteins in signaling and development. https://doi.org/10.1101/015826)
library(CHNOSZ)
-# to reproduce the calculations in the paper, use superseded data for [Gly] and [UPBB] 20190206
+# To reproduce the calculations in the paper, use superseded data for [Gly] and [UPBB] 20190206
mod.OBIGT("[Gly]", G = -6075, H = -5570, S = 17.31)
mod.OBIGT("[UPBB]", G = -21436, H = -45220, S = 1.62)
@@ -10,7 +10,7 @@
pname <- c("SHH", "OLIG2", "NKX22", "FOXA2", "IRX3", "PAX6", "NKX62", "DBX1",
"DBX2", "NKX61", "PAX7", "GLI1", "GLI2", "GLI3", "PTC1", "SMO", "GLI3R")[1:11]
-# colors modified from Dessaud et al., 2008 and Hui and Angers, 2011
+# Colors modified from Dessaud et al., 2008 and Hui and Angers, 2011
fill <- c(SHH = "#c8c7c8", OLIG2 = "#f9a330", NKX22 = "#ef2b2c", FOXA2 = "#6ab0b9",
NKX61 = "#b76775", DBX2 = "#35bcba", PAX7 = "#f6ef42",
PAX6 = "#4d2a59", IRX3 = "#63c54e", NKX62 = "#f24e33", DBX1 = "#d4e94e",
@@ -17,7 +17,7 @@
PTC1 = "#c7b0ee", SMO = "#8fd4ef", GLI1 = "#fcdc7e", GLI2 = "#c7e3b0", GLI3 = "#fcdc7e",
GLI3R = "#f1b1ae")
-# names for plotting
+# Names for plotting
names <- c(SHH = "Shh", OLIG2 = "Olig2", NKX22 = "Nkx2.2", FOXA2 = "Foxa2",
NKX61 = "Nkx6.1", DBX2 = "Dbx2", PAX7 = "Pax7",
PAX6 = "Pax6", IRX3 = "Irx3", NKX62 = "Nkx6.2", DBX1 = "Dbx1",
@@ -24,24 +24,24 @@
PTC1 = "Ptch1", SMO = "Smo", GLI1 = "Gli1A", GLI2 = "Gli2A", GLI3 = "Gli3A",
GLI3R = "Gli3R")
-# protein indices of Shh and the transcription factors
+# Protein indices of Shh and the transcription factors
ip <- match(pname, thermo()$protein$protein)
aa <- thermo()$protein[ip, ]
-# set up basis species
+# Set up basis species
basis("CHNOS")
basis("NH3", -7)
-# save as PDF file?
+# Save as PDF file?
pdf <- FALSE
-# draw interpretive legend?
+# Draw interpretive legend?
interp <- TRUE
-# set up ranges of logfO2 and logaH2O
+# Set up ranges of logfO2 and logaH2O
O2 <- seq(-70, -100, length.out = 500)
H2O <- seq(0.5, -4.5, length.out = 500)
A <- affinity(H2O = H2O, O2 = O2, iprotein = ip)
-# plot affinities per residue, compared to SHH
+# Plot affinities per residue, compared to SHH
pl <- protein.length(ip)
names(A$values) <- pname
for(i in 1:length(A$values)) A$values[[i]] <- A$values[[i]] / pl[i]
@@ -49,77 +49,77 @@
for(i in 1:length(A$values)) A$values[[i]] <- A$values[[i]] - A.SSH
ylab <- expression(bold(A)/2.303*italic(RT)*" vs Shh")
xlab <- expression(log*italic(a)[H[2]][O])
-# set up normal plot, or plot with interpretive drawings
+# Set up normal plot, or plot with interpretive drawings
opar <- par(no.readonly = TRUE)
par(mar = c(5.1, 4.1, 4.1, 2.1))
if(interp) {
- if(pdf) pdf("tfactor_interp.pdf", width=6, height=6)
+ if(pdf) pdf("tfactor_interp.pdf", width = 6, height = 6)
plot.new()
- plot.window(rev(range(H2O)), c(-0.7, 4.7), xaxs="i", yaxs="i")
+ plot.window(rev(range(H2O)), c(-0.7, 4.7), xaxs = "i", yaxs = "i")
par(xpd=FALSE)
clip(range(H2O)[1], range(H2O)[2], -0.3, 4.2)
} else {
- if(pdf) pdf("tfactor_affinity.pdf", width=6, height=6)
- plot(range(H2O), c(-0.5, 4.5), type="n", xaxs="i", yaxs="i", xlab=xlab, ylab=ylab, xlim=rev(range(H2O)), mgp=c(2.5, 1, 0))
+ if(pdf) pdf("tfactor_affinity.pdf", width = 6, height = 6)
+ plot(range(H2O), c(-0.5, 4.5), type = "n", xaxs = "i", yaxs = "i", xlab = xlab, ylab = ylab, xlim = rev(range(H2O)), mgp = c(2.5, 1, 0))
}
for(i in 1:length(pl)) {
lty <- 3
lwd <- 1
- # highlight SHH with solid line
- if(pname[i]=="SHH") {
+ # Highlight SHH with solid line
+ if(pname[i] == "SHH") {
lty <- 1
lwd <- 3
}
- lines(H2O, A$values[[i]], lty=lty, lwd=lwd)
+ lines(H2O, A$values[[i]], lty = lty, lwd = lwd)
}
-# highlight lines for reaction sequence from OLIG2 to SHH
+# Highlight lines for reaction sequence from OLIG2 to SHH
A <- A$values
names(A) <- pname
# Olig2
iOLIG2 <- A$OLIG2 > A$SHH
-lines(H2O[iOLIG2], A$OLIG2[iOLIG2], col=fill["OLIG2"], lwd=3)
+lines(H2O[iOLIG2], A$OLIG2[iOLIG2], col = fill["OLIG2"], lwd = 3)
# Foxa2 with offset to distinguish it from Nkx6.1/Dbx2
iFOXA2 <- A$FOXA2 > A$OLIG2 & A$FOXA2 > A$NKX22
-lines(H2O[iFOXA2], A$FOXA2[iFOXA2], col=fill["FOXA2"], lwd=3)
+lines(H2O[iFOXA2], A$FOXA2[iFOXA2], col = fill["FOXA2"], lwd = 3)
# Nkx2.2
iNKX22 <- A$NKX22 > A$FOXA2 & A$NKX22 > A$SHH
-lines(H2O[iNKX22], A$NKX22[iNKX22], col=fill["NKX22"], lwd=3)
+lines(H2O[iNKX22], A$NKX22[iNKX22], col = fill["NKX22"], lwd = 3)
# Pax6
iPAX6 <- A$PAX6 > A$SHH
-lines(H2O[iPAX6], A$PAX6[iPAX6], col=fill["PAX6"], lwd=3)
+lines(H2O[iPAX6], A$PAX6[iPAX6], col = fill["PAX6"], lwd = 3)
# Nkx6.1 with overstepping
iNKX61 <- A$NKX61 > A$DBX2
imax <- max(which(iNKX61))
iNKX61[1: (imax-20)] <- FALSE
-lines(H2O[iNKX61], A$NKX61[iNKX61], col=fill["NKX61"], lwd=3)
+lines(H2O[iNKX61], A$NKX61[iNKX61], col = fill["NKX61"], lwd = 3)
# Dbx2
iDBX2 <- A$DBX2 > A$NKX61 & A$DBX2 < A$IRX3
-lines(H2O[iDBX2], A$DBX2[iDBX2], col=fill["DBX2"], lwd=3)
+lines(H2O[iDBX2], A$DBX2[iDBX2], col = fill["DBX2"], lwd = 3)
# Irx3
iIRX3 <- A$IRX3 > A$NKX62 & A$IRX3 > A$OLIG2
-lines(H2O[iIRX3], A$IRX3[iIRX3], col=fill["IRX3"], lwd=3)
+lines(H2O[iIRX3], A$IRX3[iIRX3], col = fill["IRX3"], lwd = 3)
# Nkx6.2
iNKX62 <- A$NKX62 > A$IRX3 & A$NKX62 > A$DBX1
-lines(H2O[iNKX62], A$NKX62[iNKX62], col=fill["NKX62"], lwd=3)
+lines(H2O[iNKX62], A$NKX62[iNKX62], col = fill["NKX62"], lwd = 3)
# Dbx1
iDBX1 <- A$DBX1 > A$NKX62 & A$DBX1 > A$SHH
-lines(H2O[iDBX1], A$DBX1[iDBX1], col=fill["DBX1"], lwd=3)
+lines(H2O[iDBX1], A$DBX1[iDBX1], col = fill["DBX1"], lwd = 3)
# Shh
#iSHH <- A$SHH > A$DBX1
-#lines(H2O[iSHH], A$SHH[iSHH], col=fill["SHH"], lwd=3)
-# the lines need names
+#lines(H2O[iSHH], A$SHH[iSHH], col = fill["SHH"], lwd = 3)
+# The lines need names
if(interp) {
- # remove plot clip region
- par(xpd=NA)
- text(-2.12, -0.48, "Nkx2.2", srt=90)
- text(-0.87, -0.65, "Olig2", srt=90)
- text(0, -0.72, "Pax6", srt=90)
- text(0.06, 4.3, "Olig2", srt=90)
- text(-0.13, 4.3, "Irx3", srt=90)
- text(-0.77, 4.3, "Nkx6.1", srt=90)
- text(-.97, 4.3, "Dbx2", srt=90)
- text(-3.45, 4.3, "Nkx6.2", srt=90)
- text(-3.65, 4.3, "Dbx1", srt=90)
+ # Remove plot clip region
+ par(xpd = NA)
+ text(-2.12, -0.48, "Nkx2.2", srt = 90)
+ text(-0.87, -0.65, "Olig2", srt = 90)
+ text(0, -0.72, "Pax6", srt = 90)
+ text(0.06, 4.3, "Olig2", srt = 90)
+ text(-0.13, 4.3, "Irx3", srt = 90)
+ text(-0.77, 4.3, "Nkx6.1", srt = 90)
+ text(-.97, 4.3, "Dbx2", srt = 90)
+ text(-3.45, 4.3, "Nkx6.2", srt = 90)
+ text(-3.65, 4.3, "Dbx1", srt = 90)
} else {
text(-1.5, 0.5, "Nkx2.2")
text(-0.1, 0.3, "Pax6")
@@ -132,53 +132,53 @@
}
text(0.3, -0.15, "Shh")
text(-4.25, 0.15, "Shh")
-text(-0.47, 0.5, "Pax7", srt=-35)
-text(-0.22, 2, "Foxa2", srt=-61)
-# are we making an interpretive plot?
+text(-0.47, 0.5, "Pax7", srt = -35)
+text(-0.22, 2, "Foxa2", srt = -61)
+# Are we making an interpretive plot?
if(interp) {
- # the left,bottom x,y-position and horizontal width of the bottom gradient wedge
+ # The left,bottom x,y-position and horizontal width of the bottom gradient wedge
xbot <- H2O[1]
ybot <- -1.4
wbot <- 3
- # the height of the bottom gradient wedge as a function of the x position
+ # The height of the bottom gradient wedge as a function of the x position
hbot <- function(x) 0.3 + 0.08*(xbot - x)
lines(c(xbot, xbot-wbot), ybot+hbot(c(xbot, xbot-wbot)))
- # draw the base and sides of the bottom gradient wedge
+ # Draw the base and sides of the bottom gradient wedge
lines(c(xbot, xbot-wbot), c(ybot, ybot))
lines(rep(xbot, 2), c(ybot, ybot+hbot(xbot)))
lines(rep(xbot, 2)-wbot, c(ybot, ybot+hbot(xbot-wbot)))
- # draw drop lines from reactions between TFs and Shh
+ # Draw drop lines from reactions between TFs and Shh
xNKX22 <- H2O[max(which(A$NKX22 > A$SHH))]
- lines(rep(xNKX22, 2), c(ybot, 0), lty=2)
+ lines(rep(xNKX22, 2), c(ybot, 0), lty = 2)
xOLIG2 <- H2O[max(which(A$OLIG2 > A$SHH))]
- lines(rep(xOLIG2, 2), c(ybot, 0), lty=2)
+ lines(rep(xOLIG2, 2), c(ybot, 0), lty = 2)
xPAX6 <- H2O[max(which(A$PAX6 > A$SHH))]
- lines(rep(xPAX6, 2), c(ybot, 0), lty=2)
- # the left,bottom x,y-position and horizontal width of the top gradient wedge
+ lines(rep(xPAX6, 2), c(ybot, 0), lty = 2)
+ # The left,bottom x,y-position and horizontal width of the top gradient wedge
xtop <- H2O[2]
ytop <- 4.7
wtop <- 5
- # the height of the top gradient wedge as a function of the x position
+ # The height of the top gradient wedge as a function of the x position
htop <- function(x) 0.4 - 0.08*(xtop - x)
lines(c(xtop, xtop-wtop), ytop+htop(c(xtop, xtop-wtop)))
- # draw the base and sides of the top gradient wedge
+ # Draw the base and sides of the top gradient wedge
lines(c(xtop, xtop-wtop), c(ytop, ytop))
lines(rep(xtop, 2), c(ytop, ytop+htop(xtop)))
lines(rep(xtop, 2)-wtop, c(ytop, ytop+htop(xtop-wtop)))
- # draw drop lines to reactions between TFs
+ # Draw drop lines to reactions between TFs
iIRX3 <- min(which(A$IRX3 > A$OLIG2))
- lines(rep(H2O[iIRX3], 2), c(A$IRX3[iIRX3], ytop+htop(H2O[iIRX3])), lty=2)
+ lines(rep(H2O[iIRX3], 2), c(A$IRX3[iIRX3], ytop+htop(H2O[iIRX3])), lty = 2)
iDBX2 <- min(which(A$DBX2 > A$NKX61))
- lines(rep(H2O[iDBX2], 2), c(A$DBX2[iDBX2], ytop+htop(H2O[iDBX2])), lty=2)
+ lines(rep(H2O[iDBX2], 2), c(A$DBX2[iDBX2], ytop+htop(H2O[iDBX2])), lty = 2)
iDBX1 <- min(which(A$DBX1 > A$NKX62))
- lines(rep(H2O[iDBX1], 2), c(A$DBX1[iDBX1], ytop+htop(H2O[iDBX1])), lty=2)
- # indicate plot variables
+ lines(rep(H2O[iDBX1], 2), c(A$DBX1[iDBX1], ytop+htop(H2O[iDBX1])), lty = 2)
+ # Indicate plot variables
arrows(-2.7, 2, -2.7, 2.5, 0.1)
- text(-2.8, 2.3, "affinity\nvs. Shh", adj=0)
+ text(-2.8, 2.3, "affinity\nvs. Shh", adj = 0)
arrows(-2.7, 2, -2.27, 2, 0.1)
- text(-2.3, 1.8, expression(list(log*italic(f)[O[2]], )), adj=0)
- text(-2.3, 1.6, expression(list(log*italic(a)[H[2]*O])), adj=0)
- # label neural progenitors
+ text(-2.3, 1.8, expression(list(log*italic(f)[O[2]], )), adj = 0)
+ text(-2.3, 1.6, expression(list(log*italic(a)[H[2]*O])), adj = 0)
+ # Label neural progenitors
text(-2.37, -1.55, "FP")
text(-1.6, -1.55, "p3")
text(-0.53, -1.55, "pMN")
@@ -187,20 +187,20 @@
text(-0.47, 5.25, "p2")
text(-2.2, 5.25, "p1")
text(-4, 5.25, "p0")
- # label Shh gradient and stages
- text(1.2, -1.2, "Shh\ngradient", adj=0)
- text(1.2, 4.9, "Shh\ngradient", adj=0)
- text(-2.6, -0.6, expression(italic("Stage 1: loading")), adj=0)
- text(-1.5, 4.3, expression(italic("Stage 2: unloading")), adj=0)
+ # Label Shh gradient and stages
+ text(1.2, -1.2, "Shh\ngradient", adj = 0)
+ text(1.2, 4.9, "Shh\ngradient", adj = 0)
+ text(-2.6, -0.6, expression(italic("Stage 1: loading")), adj = 0)
+ text(-1.5, 4.3, expression(italic("Stage 2: unloading")), adj = 0)
} else {
- # add second axis: logfO2
+ # Add second axis: logfO2
pu <- par("usr")
pu[1:2] <- rev(range(O2))
- par(usr=pu)
- axis(3, at=seq(-75, -105, by=-5))
- mtext(expression(log*italic(f)[O[2]]), line=2)
+ par(usr = pu)
+ axis(3, at = seq(-75, -105, by = -5))
+ mtext(expression(log*italic(f)[O[2]]), line = 2)
}
-# all done!
+# All done!
par(opar)
if(pdf) dev.off()
reset()
Modified: pkg/CHNOSZ/demo/TCA.R
===================================================================
--- pkg/CHNOSZ/demo/TCA.R 2023-03-04 02:41:46 UTC (rev 774)
+++ pkg/CHNOSZ/demo/TCA.R 2023-03-04 03:43:54 UTC (rev 775)
@@ -7,7 +7,7 @@
# These plots use calories 20220325
E.units("cal")
-# species in reactions
+# Species in reactions
NADox <- "NAD(ox)-"; NADred <- "NAD(red)-2"
ADP <- "ADP-3"; ATP <- "ATP-4"
species <- list(
@@ -21,7 +21,7 @@
c("malate-2", NADox, "oxaloacetate-2", NADred, "H+"),
c("pyruvate", NADox, ADP, "HPO4-2", "H2O", "CO2", NADred, "H+", ATP, "H2")
)
-# reaction coefficients
+# Reaction coefficients
coeffs <- list(
c(-1, -1, -1, -1, 1, 1, 1, 1),
c(-1, 1, 1),
@@ -33,7 +33,7 @@
c(-1, -1, 1, 1, 1),
c(-1, -4, -1, -1, -2, 3, 4, 2, 1, 1)
)
-# species names
+# Species names
oxal <- quote(Oxaloacetate^-2)
pyr <- quote(Pyruvate^-"")
h2o <- quote(H[2]*O)
@@ -54,10 +54,10 @@
mal <- quote(Malate^-2)
# the reaction double arrow
eq <- "\u21cc"
-sublist <- list(oxal=oxal, pyr=pyr, h2o=h2o, nox=nox, cit=cit, nred=nred,
- co2=co2, hplus=hplus, aco=aco, iso=iso, ket=ket, adp=adp,
- hpo4=hpo4, suc=suc, atp=atp, fum=fum, h2=h2, mal=mal, eq=eq)
-# reaction titles
+sublist <- list(oxal = oxal, pyr = pyr, h2o = h2o, nox = nox, cit = cit, nred = nred,
+ co2 = co2, hplus = hplus, aco = aco, iso = iso, ket = ket, adp = adp,
+ hpo4 = hpo4, suc = suc, atp = atp, fum = fum, h2 = h2, mal = mal, eq = eq)
+# Reaction titles
rtitle <- list(
c(substitute(" "*oxal + pyr + h2o + nox ~eq~ "", sublist), substitute(cit + nred + co2 + hplus, sublist)),
substitute(cit ~eq~ aco + h2o, sublist),
@@ -70,34 +70,34 @@
c(substitute(pyr + 4*nox + adp + hpo4 + 2*h2o ~eq~ " ", sublist),
substitute(3*co2 + 4*nred + 2*hplus + atp + h2 * " ", sublist))
)
-# set up plot
+# Set up plot
opar <- par(no.readonly = TRUE)
-par(mfrow=c(3, 3))
+par(mfrow = c(3, 3))
ylims <- list(
c(-10, 45), c(1, 6), c(-2.5, 7.5),
c(-35, 5), c(-9, 5), c(5, 28),
c(-1.5, 6), c(14, 18), c(20, 80)
)
-# loop over reactions
+# Loop over reactions
for(i in seq_along(species)) {
- thermo.plot.new(xlim=c(0, 500), ylim=ylims[[i]], xlab=axis.label("T"),
- ylab=axis.label("DrG0", prefix="k"), mar=c(3.0, 3.5, 3.5, 2.0))
- # loop over isobars
+ thermo.plot.new(xlim = c(0, 500), ylim = ylims[[i]], xlab = axis.label("T"),
+ ylab = axis.label("DrG0", prefix = "k"), mar = c(3.0, 3.5, 3.5, 2.0))
+ # Loop over isobars
for(P in seq(500, 5000, 500)) {
T <- seq(0, 500, 10)
if(P==500) T <- seq(0, 350, 10)
if(P==5000) T <- seq(100, 500, 10)
- # calculate and plot standard Gibbs energy
- sout <- subcrt(species[[i]], coeffs[[i]], T=T, P=P)$out
+ # Calculate and plot standard Gibbs energy
+ sout <- subcrt(species[[i]], coeffs[[i]], T = T, P = P)$out
lines(T, sout$G/1000)
}
- if(is.list(rtitle[[i]])) mtitle(as.expression(rtitle[[i]]), spacing = 1.6, cex=0.8)
- else mtitle(as.expression(rtitle[[i]]), line=0.4, cex=0.8)
+ if(is.list(rtitle[[i]])) mtitle(as.expression(rtitle[[i]]), spacing = 1.6, cex = 0.8)
+ else mtitle(as.expression(rtitle[[i]]), line = 0.4, cex = 0.8)
}
-# make an overall title
-par(xpd=NA)
-text(-70, 284, "Citric Acid Cycle, after Canovas and Shock, 2016", font=2, cex=1.5)
-par(xpd=FALSE)
+# Make an overall title
+par(xpd = NA)
+text(-70, 284, "Citric Acid Cycle, after Canovas and Shock, 2016", font = 2, cex = 1.5)
+par(xpd = FALSE)
par(opar)
# Reset the units
Modified: pkg/CHNOSZ/demo/affinity.R
===================================================================
--- pkg/CHNOSZ/demo/affinity.R 2023-03-04 02:41:46 UTC (rev 774)
+++ pkg/CHNOSZ/demo/affinity.R 2023-03-04 03:43:54 UTC (rev 775)
@@ -1,5 +1,5 @@
-## affinities of metabolic reactions
-## after Amend and Shock, 2001, Fig. 7
+## Affinities of metabolic reactions
+## After Amend and Shock, 2001, Fig. 7
## Amend, J. P. and Shock, E. L. (2001) Energetics of overall metabolic reactions of thermophilic and hyperthermophilic Archaea and Bacteria.
## FEMS Microbiol. Rev. 25, 175--243. https://doi.org/10.1016/S0168-6445(00)00062-0
library(CHNOSZ)
Modified: pkg/CHNOSZ/demo/aluminum.R
===================================================================
--- pkg/CHNOSZ/demo/aluminum.R 2023-03-04 02:41:46 UTC (rev 774)
+++ pkg/CHNOSZ/demo/aluminum.R 2023-03-04 03:43:54 UTC (rev 775)
@@ -3,38 +3,38 @@
# 20190415 add diagrams from Tutolo et al., 2014
library(CHNOSZ)
-## set up plotting area
+## Set up plotting area
opar <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 775
More information about the CHNOSZ-commits
mailing list