[CHNOSZ-commits] r148 - in pkg/CHNOSZ: . R data inst inst/extdata/cpetc man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Feb 13 17:54:47 CET 2017
Author: jedick
Date: 2017-02-13 17:54:47 +0100 (Mon, 13 Feb 2017)
New Revision: 148
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/data/OBIGT.csv
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/inst/extdata/cpetc/rubisco.csv
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/man/protein.Rd
pkg/CHNOSZ/man/read.expr.Rd
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/equilibrium.Rnw
pkg/CHNOSZ/vignettes/equilibrium.lyx
pkg/CHNOSZ/vignettes/vig.bib
Log:
anintro.Rmd: add activity coefficients
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-02-13 06:42:44 UTC (rev 147)
+++ pkg/CHNOSZ/DESCRIPTION 2017-02-13 16:54:47 UTC (rev 148)
@@ -1,6 +1,6 @@
Date: 2017-02-13
Package: CHNOSZ
-Version: 1.0.8-37
+Version: 1.0.8-38
Title: Chemical Thermodynamics and Activity Diagrams
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2017-02-13 06:42:44 UTC (rev 147)
+++ pkg/CHNOSZ/R/diagram.R 2017-02-13 16:54:47 UTC (rev 148)
@@ -22,7 +22,7 @@
# colors
col=par("col"), col.names=par("col"), fill=NULL,
# labels
- names=NULL, main=NULL, legend.x="topright", format.names=TRUE,
+ names=NULL, main=NULL, legend.x=NA, format.names=TRUE,
# plotting controls
add=FALSE, plot.it=TRUE, tplot=TRUE, ...
) {
@@ -240,8 +240,6 @@
}
# draw the lines
for(i in 1:length(plotvals)) lines(xvalues, plotvals[[i]], col=col[i], lty=lty[i], lwd=lwd[i])
- # turn off legend if too many species
- if(ngroups > 10 & missing(legend.x)) legend.x <- NA
if(!add & !is.null(legend.x)) {
# 20120521: use legend.x=NA to label lines rather than make legend
if(is.na(legend.x)) {
Modified: pkg/CHNOSZ/data/OBIGT.csv
===================================================================
--- pkg/CHNOSZ/data/OBIGT.csv 2017-02-13 06:42:44 UTC (rev 147)
+++ pkg/CHNOSZ/data/OBIGT.csv 2017-02-13 16:54:47 UTC (rev 148)
@@ -1721,7 +1721,7 @@
MgADP-,MgADP-,MgC10H12N5O10P2-,aq,LH06b,NA,21.May.03,-567524,-730933,55.74,139.985659655832,189.17,2.4109,127.142,-19.16,-2.477,175.24,-3.8947,3,-1
MgHADP,MgHADP,MgC10H13N5O10P2,aq,LH06b,NA,10.Jun.03,-574465,-731643,76.54,161.472275334608,199.93,2.4625,137.497,-19.3997,-2.9051,244.57,-19.376,4.8,0
Mg2ADP+,Mg2ADP+,Mg2C10H12N5O10P2+,aq,LH06b,NA,6.Jun.03,-678301,-839029,41.31,159.440726577438,174.77,2.3411,113.105,-18.834,-1.8968,166.84,-1.3991,0.5,1
-MgATP-2,MgATP-2,MgC10H12N5O13P3-2,aq,LH06b,NA,24.Sep.06,-773608.986615679,-966278.680688336,56.5965583173996,150.143403441683,193.48,2.45028680688337,135.083891013384,-19.343690248566,-2.80544933078394,256.022944550669,-26.0994263862333,5.73613766730401,-1
+MgATP-2,MgATP-2,MgC10H12N5O13P3-2,aq,LH06b,NA,24.Sep.06,-773608.986615679,-966278.680688336,56.5965583173996,150.143403441683,193.48,2.45028680688337,135.083891013384,-19.343690248566,-2.80544933078394,256.022944550669,-26.0994263862333,5.73613766730401,-2
MgHATP-,MgHATP-,MgC10H13N5O13P3-,aq,LH06b,NA,11.Jun.03,-780995,-965191,85.01,210.803059273423,208.97,2.5011,145.27,-19.58,-3.2265,315.61,-25.582,5.8,-1
MgH2ATP,MgH2ATP,MgC10H14N5O13P3,aq,LH06b,NA,11.Jun.03,-786200,-970481,84.72,267.208413001912,208.63,2.4936,143.766,-19.545,-3.1643,362.6,-24.126,5.1,0
Mg2ATP,Mg2ATP,Mg2C10H12N5O13P3,aq,LH06b,NA,16.Jun.03,-884603,-1072754,48.32,204.134799235182,191.55,2.4505,135.089,-19.344,-2.8056,310.31,-23.113,6.5,0
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-02-13 06:42:44 UTC (rev 147)
+++ pkg/CHNOSZ/inst/NEWS 2017-02-13 16:54:47 UTC (rev 148)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.8-37 (2017-02-13)
+CHANGES IN CHNOSZ 1.0.8-38 (2017-02-13)
---------------------------------------
DOCUMENTATION:
@@ -58,6 +58,9 @@
- read.fasta(): change argument `i` to `iseq`; this is used to select
particular sequences to read from the file.
+- diagram(): default of `legend.x` changed to NA (instead of making a
+ legend, put labels near the lines).
+
CLEANUP AND BUG FIXES:
- subcrt() returns `loggam` using the common logarithm; add
@@ -76,6 +79,8 @@
- Correct charge (-2) of NAD(red)-2 in OBIGT.csv. Thanks to Peter
Canovas.
+- Also correct charge (-2) of MgATP-2.
+
- Remove msgout(), and replace previous calls to that function with
message() from base R. As a result, the messages don't appear in Sweave
vignettes, but can be turned on or off in knitr vignettes.
Modified: pkg/CHNOSZ/inst/extdata/cpetc/rubisco.csv
===================================================================
--- pkg/CHNOSZ/inst/extdata/cpetc/rubisco.csv 2017-02-13 06:42:44 UTC (rev 147)
+++ pkg/CHNOSZ/inst/extdata/cpetc/rubisco.csv 2017-02-13 16:54:47 UTC (rev 148)
@@ -10,7 +10,7 @@
L0RHZ1,35,35,B,"Desulfovibrio hydrothermalis",http://dx.doi.org/10.1099/ijs.0.02323-0
Q8THG2,35,40,A,"Methanosarcina acetivorans",http://aem.asm.org/content/47/5/971
F9ZLP0,45,45,B,"Acidithiobacillus caldus",http://dx.doi.org/10.1099/13500872-140-12-3451
-P37393,45,45,E,"Cyanidium caldarium",http://dx.doi.org/10.1007/BF00409031
+P37393,45,45,E,"Cyanidium caldarium",http://link.springer.com/article/10.1007/BF00409031
P72383,45,50,B,"Sulfobacillus acidophilus",http://dx.doi.org/10.1099/00221287-142-4-775
Q51856,52,52,B,"Pseudomonas hydrogenothermophila",http://dx.doi.org/10.1271/bbb1961.41.685
Q2JIP3,50,55,B,"Synechococcus sp. (strain JA-2-3B'a(2-13))",http://dx.doi.org/10.1128/AEM.72.1.544-550.2006
@@ -18,7 +18,7 @@
Q8DIS5,57,57,B,"Thermosynechococcus elongatus",http://dx.doi.org/10.1093/oxfordjournals.pcp.a075684
G8LZL2,60,60,B,"Clostridium clariflavum",http://dx.doi.org/10.1099/ijs.0.003483-0
F8IID7,65,65,B,"Bacillus acidocaldarius",http://dx.doi.org/10.1099/00221287-67-1-9
-A8F7V4,65,65,B,"Thermotoga lettingae",http://dx.doi.org/10.1099/ijs.0.02165-0
+A8F7V4,65,65,B,"Thermotoga lettingae",http://dx.doi.org/10.1099/00207713-52-4-1361
B9KXE5,70,70,B,"Thermomicrobium roseum",http://dx.doi.org/10.1371/journal.pone.0004207
O28635,76,76,A,"Archaeoglobus fulgidus",http://aem.asm.org/content/60/4/1227
Q58632,85,85,A,"Methanocaldococcus jannaschii",http://dx.doi.org/10.1007/BF00425213
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2017-02-13 06:42:44 UTC (rev 147)
+++ pkg/CHNOSZ/man/diagram.Rd 2017-02-13 16:54:47 UTC (rev 148)
@@ -15,7 +15,7 @@
cex=par("cex"), cex.names=1, cex.axis=par("cex"),
lty=NULL, lwd=par("lwd"), dotted=NULL,
col=par("col"), col.names=par("col"), fill=NULL,
- names=NULL, main=NULL, legend.x="topright", format.names=TRUE,
+ names=NULL, main=NULL, legend.x=NA, format.names=TRUE,
add=FALSE, plot.it=TRUE, tplot=TRUE, ...)
strip(affinity, ispecies = NULL, col = NULL, ns = NULL,
xticks = NULL, ymin = -0.2, xpad = 1, cex.names = 0.7)
@@ -78,8 +78,8 @@
The names of the list are used to label the lines or fields for the summed activities of the resulting groups.
For 1-D diagrams, the default setting for the y-axis is a logarithmic scale (unless \code{alpha} is TRUE) with limits corresponding to the range of logarithms of activities (or 0,1 if \code{alpha} is TRUE); these actions can be overridden by \code{ylog} and \code{ylim}.
-A \code{\link{legend}} is placed at the location identified by \code{legend.x}, or omitted if \code{legend.x} is \code{FALSE}.
-If \code{legend.x} is NA, instead of a legend, the lines are labeled with the names of the species near the maximum value.
+If \code{legend.x} is NA (the default), the lines are labeled with the names of the species near the maximum value.
+Otherwise, a \code{\link{legend}} is placed at the location identified by \code{legend.x}, or omitted if \code{legend.x} is NULL.
The line type and line width can be controlled with \code{lty} and \code{lwd}, respectively.
On 2-D diagrams, the fields represent the species with the highest equilibrium activity.
Modified: pkg/CHNOSZ/man/protein.Rd
===================================================================
--- pkg/CHNOSZ/man/protein.Rd 2017-02-13 06:42:44 UTC (rev 147)
+++ pkg/CHNOSZ/man/protein.Rd 2017-02-13 16:54:47 UTC (rev 148)
@@ -58,7 +58,7 @@
a <- affinity(O2=c(-100, -65))
# try normalize=FALSE to make Fig. 5a in the paper
e <- equilibrate(a, normalize=TRUE)
-d <- diagram(e, ylim=c(-5, -1), legend.x=NA, names=organisms)
+d <- diagram(e, ylim=c(-5, -1), names=organisms)
# add water stability line
abline(v=-83.1, lty=2)
title(main="Surface-layer proteins, after Dick, 2008")
@@ -87,7 +87,7 @@
col <- rep("red", length(prot))
col[ZC > 0] <- "blue"
e <- equilibrate(a, normalize=TRUE)
-d <- diagram(e, col=col, legend.x=NA)
+d <- diagram(e, col=col, format.names=FALSE)
title(main="Bovine proteins, GSH/GSSG redox buffer")
}
Modified: pkg/CHNOSZ/man/read.expr.Rd
===================================================================
--- pkg/CHNOSZ/man/read.expr.Rd 2017-02-13 06:42:44 UTC (rev 147)
+++ pkg/CHNOSZ/man/read.expr.Rd 2017-02-13 16:54:47 UTC (rev 148)
@@ -102,15 +102,10 @@
a <- affinity(O2=c(-82, -65))
e <- equilibrate(a, loga.balance=0, normalize=TRUE)
mycolor <- topo.colors(length(locations))
-diagram(e, names=locations, ylim=c(-5, -3), legend.x=NA, col=mycolor, lwd=2)
+diagram(e, names=locations, ylim=c(-5, -3), col=mycolor, lwd=2)
dp <- describe.property(c("T", "P"), c(25, 1))
db <- describe.basis(ibasis=(1:6)[-5])
legend("topright", legend=c(dp, db), bty="n")
-## Notable features include: proteins in the early Golgi, endoplasmic reticulum (ER)
-## and vacuole are chemically the most stable relative to those in other locations in the cell;
-## proteins in the vacuole are very oxidized; the stability of proteins decreases going from
-## early Golgi to Golgi to late Golgi; the least stable proteins (or most energetic)
-## are found in the microtubules and bud neck.
#############################
## examples using stress() ##
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2017-02-13 06:42:44 UTC (rev 147)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2017-02-13 16:54:47 UTC (rev 148)
@@ -79,8 +79,8 @@
For more information on R, see "An Introduction to R" and other documents in [The R Manuals](http://cran.r-project.org/manuals.html).
CHNOSZ has been developed since 2006 to support research projects in geochemistry and compositional biology.
-The package provides functions and an extensive thermodynamic database that can be used to calculate the stoichiometric and energetic properties of reactions among minerals and organic or inorganic aqueous species.
-These same functions enable calculations of chemical affinities and metastable equilibrium for systems of proteins.
+The package provides functions and a thermodynamic database that can be used to calculate the stoichiometric and energetic properties of reactions among minerals and organic or inorganic aqueous species.
+These same functions enable calculations of chemical affinities and metastable equilibrium distributions of proteins.
## Installing and loading CHNOSZ
@@ -667,14 +667,13 @@
```{r rainbow_diagram, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", echo=FALSE, message=FALSE, cache=TRUE, fig.cap="Affinities of organic synthesis in a hydrothermal system, after Shock and Canovas (2010).", pngquant=pngquant}
diagram(a, balance = 1, ylim = c(-100, 100), ylab = axis.label("A", prefix="k"),
- col = rainbow(8), lwd = 2, legend.x = NA, bg = "slategray3")
+ col = rainbow(8), lwd = 2, bg = "slategray3")
abline(h = 0, lty = 2, lwd = 2)
```
Finally, we use <span style="color:green">`diagram()`</span> to plot the results.
Although only temperature is shown on the *x* axis, pH and the activities of CO<sub>2</sub>, H<sub>2</sub>, NH<sub>4</sub><sup>+</sup>, and H<sub>2</sub>S are also varied according to the data in `rb`.
By default, <span style="color:green">`diagram()`</span> attempts to scale the affinities by dividing by the reaction coefficients of a shared basis species (in this case, CO<sub>2</sub>).
To override that behavior, we set `balance = 1` to plot the affinities of the formation reactions as written (per mole of the species being formed).
-Also, `legend.x = NA` is used to suppress making a legend (so the labels are placed next to the lines instead).
```{r rainbow_diagram, eval=FALSE}
```
@@ -712,7 +711,7 @@
unlist(affinity(T = 300, P = 100)$values)
```
-We use <span style="color:red">`mod.buffer()`</span> to choose the `cr2` phase of phyrrhotite, which is stable at this temperature.
+We use <span style="color:red">`mod.buffer()`</span> to choose the `cr2` phase of pyrrhotite, which is stable at this temperature.
Then, we set up H<sub>2</sub>S and `r o2` to be buffered by PPM, and inspect their buffered activities:
```{r PPM_setup, results="hide"}
mod.buffer("PPM", "pyrrhotite", "cr2")
@@ -778,7 +777,7 @@
diagram(a150, add = TRUE, col = "red")
e25 <- equilibrate(a25, loga.balance = -3)
e150 <- equilibrate(a150, loga.balance = -3)
-diagram(e25, ylim = c(-6, 0), bty = "n")
+diagram(e25, ylim = c(-6, 0), legend.x = "topright", bty = "n")
diagram(e150, add = TRUE, col = "red")
diagram(e25, alpha = TRUE, legend.x = "center", bty = "n")
diagram(e150, alpha = TRUE, add = TRUE, col = "red")
@@ -868,7 +867,7 @@
par(mfrow = c(1, 3))
groups <- list(inorganic = ii, alcohols = ia, ketones = ik,
`carboxylic acids` = c(ic, ica), alkenes = ie)
-diagram(e, alpha = TRUE, groups = groups, legend.x = NA)
+diagram(e, alpha = TRUE, groups = groups)
# plot only alcohols
names <- within(species(), name[-ia] <- "")$name
lty <- ifelse(names == "", 0, 1)
@@ -1038,7 +1037,7 @@
Then, in one line, we calculate the formula of the protein, followed by `r zc`.
Next, point symbols are assigned according to domain (Archaea, Bacteria, Eukaryota); numerals inside the symbols reflect the ordering by *T*<sub>opt</sub> in three temperature ranges (0--35 °C, 37.5--60 °C, and 65--100 °C).
-```{r rubisco_svg, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", fig.keep='none', fig.ext='svg', custom.plot=TRUE, embed.tag=TRUE, echo=FALSE, results="hide", message=FALSE, fig.cap='Average oxidation state of carbon in Rubisco compared with optimal growth temperature of organisms. **This is an interactive image.** Move the mouse over the points to show the names of the organisms, and click to open a reference in a new window. (Made using [RSVGTipsDevice](https://cran.r-project.org/package=RSVGTipsDevice); code differs from that shown to the left.)', cache=TRUE}
+```{r rubisco_svg, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, out.width="100%", fig.keep='none', fig.ext='svg', custom.plot=TRUE, embed.tag=TRUE, echo=FALSE, results="hide", message=FALSE, fig.cap='Average oxidation state of carbon in Rubisco compared with optimal growth temperature of organisms. **This is an interactive image.** Move the mouse over the points to show the names of the organisms, and click to open a reference in a new window. (Made using [RSVGTipsDevice](https://cran.r-project.org/package=RSVGTipsDevice) with code that can be seen in the source of this document.)', cache=TRUE}
if(require("RSVGTipsDevice")) {
datfile <- system.file("extdata/cpetc/rubisco.csv", package = "CHNOSZ")
fastafile <- system.file("extdata/fasta/rubisco.fasta", package = "CHNOSZ")
@@ -1180,9 +1179,9 @@
#mod.obigt("[Met]", G = -35245, H = -59310)
a <- affinity(O2 = c(-80, -73), iprotein = ip, loga.protein = logact)
e <- equilibrate(a)
-diagram(e, ylim = c(-5, -2), format.names = FALSE, legend.x = NA, col=1:5, lwd=2)
+diagram(e, ylim = c(-5, -2), format.names = FALSE, col = 1:5, lwd = 2)
e <- equilibrate(a, normalize = TRUE)
-diagram(e, ylim = c(-5, -2.5), format.names = FALSE, legend.x=NA, col=1:5, lwd=2)
+diagram(e, ylim = c(-5, -2.5), format.names = FALSE, col = 1:5, lwd = 2)
abline(h = logabundance, lty = 1:5, col = 1:5)
revisit(e, "DGinf", logabundance)
```
@@ -1363,8 +1362,138 @@
## Activity coefficients
-XXX Example using <span style="color:green">`nonideal()`</span> ... Alberty, 2003, Fig. 1.4-1.5 (p. 9-10, 235-236)
+<span style="color:green">`nonideal()`</span> uses the Debye-Hückel equation, described in Chapter 3 of @Alb03, to calculate activity coefficients of charged species.
+The calculations can be invoked by setting the `IS` argument in <span style="color:green">`subcrt()`</span> or <span style="color:green">`affinity()`</span>
+The activity coefficients are calculated as a function of ionic strength (*I*), temperature, and charge of the species, without any other species-specific parameters.
+Due to these assumptions as well as the limited applicability of the implemented equations to very high *I* and/or *T*, <span style="color:green">`nonideal()`</span> is presented as an experimental feature.
+Let's take a look at calculated activity coefficients at two temperatures and their effect on the standard Gibbs energies of formation (Δ*G*°<sub>*f*</sub>) of species with different charge:
+```{r subcrt_IS}
+subcrt(c("MgATP-2", "MgHATP-", "MgH2ATP"),
+ T = c(25, 100), IS = c(0, 0.25), property = "G")$out
+```
+
+The logarithms of the activity coefficients (`loggam`) are more negative for the higher-charged species, as well as at higher temperature, and have a stabilizing effect.
+That is, the apparent Gibbs energies at *I* > 0 are less than the standard Gibbs energies at *I* = 0.
+
+We can use these calculations to make some speciation plots, similar to Figures 1.2--1.5 in Alberty (2003).
+These figures show the distribution of differently charged species of adenosine triphosphate (ATP) as a function of pH, and the average number of H<sup>+</sup> and Mg<sup>+2</sup> bound to ATP in solution as a function of pH or pMg (-log*a*<sub>Mg<sup>+2</sup></sub>).
+
+Use <span style="color:green">`info()`</span> to see what ATP species are available.
+The sources of high-temperature thermodynamic data for these species are two papers by LaRowe and Helgeson [- at LH06a; - at LH06b].
+```{r info_ATP, results="hide"}
+info(" ATP")
+```
+
+The plots for this system in Alberty's book were made for *I* = 0.25 M and *T* = 25 °C.
+As a demonstration of CHNOSZ's capabilities, we set the temperature to 100 °C.
+```{r T_100}
+T <- 100
+```
+
+To begin, we set the basis species, add the variously protonated ATP species, calculate the affinities of the formation reactions, equilibrate the species, and make a degree of formation (α) or mole fraction diagram.
+This is similar to Figure 1.3 of Alberty (2003), but calculated for *I* = 0 and *T* = 100 °C:
+```{marginfigure}
+To make the code more readable, commands for plotting titles and legends are not shown.
+All of the commands are available in the source of this document.
+```
+
+```{r ATP, eval=FALSE, echo=2:6}
+par(mfrow = c(1, 4), mar = c(3.1, 3.6, 2.1, 1.6), mgp = c(1.8, 0.5, 0))
+basis("MgCHNOPS+")
+species(c("ATP-4", "HATP-3", "H2ATP-2", "H3ATP-", "H4ATP"))
+a <- affinity(pH = c(3, 9), T = T)
+e <- equilibrate(a)
+d <- diagram(e, alpha = TRUE, tplot = FALSE)
+title(main = describe.property("T", T))
+alphas <- do.call(rbind, d$plotvals)
+nH <- alphas * 0:4
+Hlab <- substitute(italic(N)[H^`+`])
+plot(a$vals[[1]], colSums(nH), type = "l", xlab = "pH", ylab=Hlab, lty=2, col=2)
+a <- affinity(pH = c(3, 9), IS = 0.25, T = T)
+e <- equilibrate(a)
+d <- diagram(e, alpha = TRUE, plot.it = FALSE)
+alphas <- do.call(rbind, d$plotvals)
+nH <- alphas * 0:4
+lines(a$vals[[1]], colSums(nH))
+legend("topright", legend = c("I = 0 M", "I = 0.25 M"), lty = 2:1, col = 2:1, cex = 0.8)
+ATP.H <- substitute("ATP and H"^`+`)
+title(main = ATP.H)
+species(c("MgATP-2", "MgHATP-", "MgH2ATP", "Mg2ATP"))
+Hplot <- function(pMg, IS = 0.25) {
+ basis("Mg+2", -pMg)
+ a <- affinity(pH = c(3, 9), IS = IS, T = T)
+ e <- equilibrate(a)
+ d <- diagram(e, alpha = TRUE, plot.it = FALSE)
+ alphas <- do.call(rbind, d$plotvals)
+ NH <- alphas * c(0:4, 0, 1, 2, 0)
+ lines(a$vals[[1]], colSums(NH), lty = 7 - pMg, col = 7 - pMg)
+}
+plot(c(3, 9), c(0, 2), type = "n", xlab = "pH", ylab = Hlab)
+lapply(2:6, Hplot)
+legend("topright", legend = paste("pMg = ", 2:6), lty = 5:1, col = 5:1, cex = 0.8)
+ATP.H.Mg <- substitute("ATP and H"^`+`~"and Mg"^`+2`)
+title(main = ATP.H.Mg)
+Mgplot <- function(pH, IS = 0.25) {
+ basis("pH", pH)
+ a <- affinity(`Mg+2` = c(-2, -7), IS = IS, T = T)
+ e <- equilibrate(a)
+ d <- diagram(e, alpha = TRUE, plot.it = FALSE)
+ alphas <- do.call(rbind, d$plotvals)
+ NMg <- alphas * species()$`Mg+`
+ lines(-a$vals[[1]], colSums(NMg), lty = 10 - pH, col = 10 - pH)
+}
+Mglab <- substitute(italic(N)[Mg^`+2`])
+plot(c(2, 7), c(0, 1.2), type = "n", xlab = "pMg", ylab = Mglab)
+lapply(3:9, Mgplot)
+legend("topright", legend = paste("pH = ", 3:9), lty = 7:1, col = 7:1, cex = 0.8)
+title(main = ATP.H.Mg)
+```
+
+Note that we have saved the numeric results of <span style="color:green">`diagram()`</span>, i.e. the degrees of formation of the species (α).
+Using that, we can calculate and plot the average number of protons bound per ATP molecule.
+To do that, we use R's `rbind()` with `do.call()` to turn `alpha` into a matrix, then multiply by the number of protons bound to each species, and sum the columns to get the total (i.e. average proton number, *N*<sub>H<sup>+</sup></sub>):
+```{r ATP, eval=FALSE, echo=8:11}
+```
+
+Adding the `IS` argument to <span style="color:green">`affinity()`</span>, we can now plot *N*<sub>H<sup>+</sup></sub> at the given ionic strength.
+Here we set `plot.it = FALSE` in `diagram()` because we use the computed α to make our own plot.
+This is similar to Figure 1.3 of Alberty (2003), but at higher temperature:
+```{r ATP, eval=FALSE, echo=12:17}
+```
+
+Next, we add the Mg<sup>+2</sup>-complexed ATP species:
+```{r ATP, eval=FALSE, echo=21}
+```
+
+Here is a function to calculate and plot *N*<sub>H<sup>+</sup></sub> for a given pMg:
+```{r ATP, eval=FALSE, echo=22:30}
+```
+
+With that function in hand, we plot the lines corresponding to pMg = 2 to 6.
+This is similar to Figure 1.4 of Alberty (2003):
+```{r ATP, eval=FALSE, echo=31:32}
+```
+
+The next function calculates and plots the average number of Mg<sup>+2</sup> bound to ATP (*N*<sub>Mg<sup>+2</sup></sub>) for a given pH.
+Here we multiply `alpha` by the coefficient on Mg<sup>+2</sup> in the formation reactions of each species, and negate log*a*<sub>Mg<sup>+2</sup></sub> (the variable used in <span style="color:green">`affinity()`</span>) to get pMg:
+```{r ATP, eval=FALSE, echo=36:44}
+```
+
+Using that function, we plot the lines corresponding to pH = 3 to 9.
+This is similar to Figure 1.5 of Alberty (2003):
+```{r ATP, eval=FALSE, echo=45:47}
+```
+
+<p>
+```{r ATP, fig.fullwidth=TRUE, fig.width=10, fig.height=2.5, dpi=ifelse(dpi==50, 50, 100), out.width="100%", echo=FALSE, message=FALSE, results="hide", fig.cap="Binding of H<sup>+</sup> and Mg<sup>+2</sup> to ATP at 100 °C and *I* = 0 M (first plot) or 0.25 M (third and fourth plots).", cache=TRUE, pngquant=pngquant}
+```
+</p>
+
+We have calculated the distribution of ATP species and average binding number of H<sup>+</sup> and Mg<sup>+2</sup> for given pH, pMg, ionic strength, and temperature.
+Other methods using equilibrium constants tabulated for specified ionic strength make an algebraic solution of the equilibrium distribution of species feasible, as shown by the Mathematica code in Alberty (2003).
+However, the required calculations can also be performed "on the fly" in CHNOSZ with reference to a conventional standard state.
+
## Optimization of chemical activities
What are the conditions that minimize the standard deviation of calculated activities of species?
@@ -1488,6 +1617,8 @@
# Functions outside of the main workflow
+XXX
+
transfer, wjd, eqdata, RH2obigt, EOSregress, anim, taxonomy
Gibbs energy minimization with amino acids: [Cob13]?
@@ -1609,7 +1740,7 @@
Here, <span style="color:green">`checkGHS()`</span> calculated the value of Δ*G*°<sub>*f*</sub> from those of Δ*H*°<sub>*f*</sub> and *S*° entered for the species and from the entropy of the elements [@CWM89] in the formula entered for the species.
If the difference between the entered and calculated values of Δ*G*°<sub>*f*</sub> is greater than 100 cal mol<sup>-1</sup>, <span style="color:green">`checkGHS()`</span> prints a message.
Here, the calculated value of Δ*G*°<sub>*f*</sub> was found to be 661 cal mol<sup>-1</sup> higher than the entered value.
-Besides double-checking the values entered for Δ*G*°<sub>*f*</sub>, Δ*H*°<sub>*f*</sub>, *S*°, and the chemical formula, the source should be consulted for further clarification.
+Besides checking for any typographical errors in the entries for Δ*G*°<sub>*f*</sub>, Δ*H*°<sub>*f*</sub>, *S*°, and the chemical formula, the source should be consulted for further clarification.
Some species in the main and secondary databases are known to have inconsistent entries.
For example, the entry for cyclohexane has *C*<sub>*P*</sub>° and *V*° that differ from those calculated from the HKF parameters.
@@ -1623,7 +1754,7 @@
The species in with detected inconsistencies in `OBIGT.csv` and `OBIGT-2.csv` are listed in the file `obigt_check.csv`.
```{r check_obigt}
file <- system.file("extdata/thermo/obigt_check.csv", package = "CHNOSZ")
-dat <- read.csv(file, as.is=TRUE)
+dat <- read.csv(file, as.is = TRUE)
nrow(dat)
```
@@ -1636,15 +1767,16 @@
As you get started writing your own code and functions that use CHNOSZ, it is not uncommon to encounter problems.
For example, mismatched data types can cause problems (the important difference between factor and character was mentioned above).
+There are many tricks to learn; becoming familiar with R's `lapply()` and `do.call()` will make your programming life easier.
-Many functions in CHNOSZ perform various degrees of "sanity checks" on the arguments and will report errors if an inconsistency is detected.
-For example, unequal lengths of variables in the transect mode of <span style="color:green">`affinity()`</span> (when the variables have more than 3 values) causes an error:
+Some functions in CHNOSZ perform "sanity checks" on the arguments and will report errors if an inconsistency is detected.
+For example, unequal lengths of variables in the transect mode of <span style="color:green">`affinity()`</span> (when the variables have more than 3 values) cause an error:
```{r affinity_error, error=TRUE, message=FALSE, results="hide"}
basis("CHNOS")
aa <- c("D", "T", "S", "E", "G", "A", "K", "H")
species(aminoacids("", aa))
-a <- affinity(O2=seq(-80, -50), T=seq(0, 100))
+a <- affinity(O2 = seq(-80, -50), T = seq(0, 100))
```
In normal operation, without errors, many functions print informative messages.
@@ -1654,15 +1786,15 @@
Checking these can help you decide if the system and calculations have the desired configuration.
Here, we fix the condition that caused the error above:
```{r message_example, results="hide"}
-a <- affinity(O2=seq(-80, -50, length.out=101), T=seq(0, 100))
+a <- affinity(O2 = seq(-80, -50, length.out = 101), T = seq(0, 100))
e <- equilibrate(a)
#diagram(e, alpha=TRUE, legend.x=NA)
```
-The messages give some useful information, such as the variables and their ranges, the "wetness" of the reactions (i.e. they contain water and/or aqueous species), the balanced basis species and balance coefficients, the logarithm of total activity of the balanced basis species, and the equilibration method.
+The messages give some useful information, such as the variables and their ranges, the "wetness" of the reactions (i.e. whether they contain water and/or aqueous species), the balanced basis species and balance coefficients, the logarithm of total activity of the balanced basis species, and the equilibration method.
The commented call to <span style="color:green">`diagram()`</span> would produce a diagram that, in this case, doesn't result in any messages.
-Here is another example, with the output hidden, but where the messages show the species in the reaction and the other states that are available for the same compounds in the database.
+Here is another example, with the results hidden, but where the messages show the species in the reaction and the other states that are available for the same compounds in the database.
The number of *T* and *P* values comes from the default arguments for <span style="color:green">`subcrt()`</span>:
```{r message_subcrt, results="hide"}
subcrt(c("C2H5OH", "O2", "CO2", "H2O"), c(-1, -3, 2, 3))
@@ -1670,7 +1802,7 @@
You may occasionally encounter programming bugs or limits of the algorithms.
The following setup breaks the reaction-matrix calculation for equilibration.
-Here, C<sub>5</sub>H<sub>9</sub>NO<sub>4</sub> (glutamic acid) is identified as the balance (the only basis species present in the formation reactions of all species).
+Here, C<sub>5</sub>H<sub>9</sub>NO<sub>4</sub> (glutamic acid) is identified as the balance (the only basis species present in all the formation reactions of species).
Both a warning and an error are generated due to missing values in a computation within <span style="color:green">`equil.reaction()`</span>:
<!-- 20170213
The warning and error don't get formatted correctly in html output.
@@ -1695,11 +1827,11 @@
# Citation and contact information
-If you use CHNOSZ for publications, it would be appreciated if you cite the following paper:
+If you use CHNOSZ for publications, please consider citing the following paper:
```{r citation_CHNOSZ, results="asis"}
-cref <- citation("CHNOSZ")
-print(cref, style = "html")
+pref <- citation("CHNOSZ")
+print(pref, style = "html")
```
If you have bug reports, questions that haven't been answered here or in the other documentation, or contributions of code or data, please contact the maintainer:
@@ -1720,16 +1852,21 @@
* 2010-09-30 Initial version.
* 2011-08-15 Add <span style="color:green">`browse.refs()`</span>; modifying database hint changed to <span style="color:blue">`help(thermo)`</span>.
```{marginfigure}
-<span style="color:green">`browse.refs()`</span> renamed to <span style="color:green">`thermo.refs()`</span> in 2017.
+<span style="color:green">`browse.refs()`</span> was renamed to <span style="color:green">`thermo.refs()`</span> in 2017.
```
* 2012-06-16 Add “More activity diagrams”.
* 2015-05-14 Add warning about internal consistency of thermodynamic data.
-* 2017-02-15 Complete rewrite; switch from Sweave to knitr (Tufte style).
+* 2017-02-15 Completely rewritten; switch from Sweave to knitr (Tufte style).
View the R Markdown source of this document [on R-Forge](https://r-forge.r-project.org/scm/viewvc.php/pkg/CHNOSZ/vignettes/anintro.Rmd?view=markup&root=chnosz) or in R:
```{r file_edit_anintro, eval=FALSE}
file.edit(system.file("doc/anintro.Rmd", package = "CHNOSZ"))
```
-
-***
+<p>
+```{r the_end}
+ ### ### ### ### ### ### ### ### ### ### ### ### ###
+### ### ### ### ### ### ### ### ### ### ### ### ### ###
+ ### ### ### ### ### ### ### ### ### ### ### ### ###
+```
+</p>
Modified: pkg/CHNOSZ/vignettes/equilibrium.Rnw
===================================================================
--- pkg/CHNOSZ/vignettes/equilibrium.Rnw 2017-02-13 06:42:44 UTC (rev 147)
+++ pkg/CHNOSZ/vignettes/equilibrium.Rnw 2017-02-13 16:54:47 UTC (rev 148)
@@ -435,7 +435,7 @@
prB <- function() {
a <- affinity(O2=c(-90, -70))
e <- equilibrate(a, balance="length", loga.balance=0)
- diagram(e, names=organisms, ylim=c(-5, -1), legend.x=NA)
+ diagram(e, names=organisms, ylim=c(-5, -1))
}
prC <- function() {
@@ -447,7 +447,7 @@
prD <- function() {
a <- affinity(O2=c(-90, -70))
e <- equilibrate(a, normalize=TRUE, loga.balance=0)
- diagram(e, names=organisms, ylim=c(-5, -1), legend.x=NA)
+ diagram(e, names=organisms, ylim=c(-5, -1))
}
prE <- function() {
@@ -459,7 +459,7 @@
prF <- function() {
a <- affinity(O2=c(-90, -70))
e <- equilibrate(a, as.residue=TRUE, loga.balance=0)
- diagram(e, names=organisms, ylim=c(-3, 1), legend.x=NA)
+ diagram(e, names=organisms, ylim=c(-3, 1))
}
@
@@ -1043,11 +1043,11 @@
a <- affinity(O2=c(-100, -65))
layout(matrix(1:2), heights=c(1, 2))
e <- equilibrate(a)
-diagram(e, ylim=c(-2.8, -1.6), legend.x=NA, names=organisms)
+diagram(e, ylim=c(-2.8, -1.6), names=organisms)
water.lines(xaxis="O2")
title(main="Equilibrium activities of proteins, normalize = FALSE")
e <- equilibrate(a, normalize=TRUE)
-diagram(e, ylim=c(-4, -2), legend.x=NA, names=organisms)
+diagram(e, ylim=c(-4, -2), names=organisms)
water.lines(xaxis="O2")
title(main="Equilibrium activities of proteins, normalize = TRUE")
@
Modified: pkg/CHNOSZ/vignettes/equilibrium.lyx
===================================================================
--- pkg/CHNOSZ/vignettes/equilibrium.lyx 2017-02-13 06:42:44 UTC (rev 147)
+++ pkg/CHNOSZ/vignettes/equilibrium.lyx 2017-02-13 16:54:47 UTC (rev 148)
@@ -1857,7 +1857,7 @@
\begin_layout Plain Layout
- diagram(e, names=organisms, ylim=c(-5, -1), legend.x=NA)
+ diagram(e, names=organisms, ylim=c(-5, -1))
\end_layout
\begin_layout Plain Layout
@@ -1915,7 +1915,7 @@
\begin_layout Plain Layout
- diagram(e, names=organisms, ylim=c(-5, -1), legend.x=NA)
+ diagram(e, names=organisms, ylim=c(-5, -1))
\end_layout
\begin_layout Plain Layout
@@ -1973,7 +1973,7 @@
\begin_layout Plain Layout
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 148
More information about the CHNOSZ-commits
mailing list