[CHNOSZ-commits] r539 - in pkg/CHNOSZ: . R inst man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sat Jul 4 04:26:50 CEST 2020
Author: jedick
Date: 2020-07-04 04:26:49 +0200 (Sat, 04 Jul 2020)
New Revision: 539
Added:
pkg/CHNOSZ/.Rbuildignore
Removed:
pkg/CHNOSZ/vignettes/hotspring.Rnw
pkg/CHNOSZ/vignettes/hotspring.lyx
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/diagram.Rd
Log:
Remove hotspring.Rnw vignette (plots will be moved to JMDplots)
Added: pkg/CHNOSZ/.Rbuildignore
===================================================================
--- pkg/CHNOSZ/.Rbuildignore (rev 0)
+++ pkg/CHNOSZ/.Rbuildignore 2020-07-04 02:26:49 UTC (rev 539)
@@ -0,0 +1 @@
+vignettes/.*.html
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2020-07-04 02:00:42 UTC (rev 538)
+++ pkg/CHNOSZ/DESCRIPTION 2020-07-04 02:26:49 UTC (rev 539)
@@ -1,6 +1,6 @@
-Date: 2020-07-02
+Date: 2020-07-04
Package: CHNOSZ
-Version: 1.3.6-13
+Version: 1.3.6-14
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/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2020-07-04 02:00:42 UTC (rev 538)
+++ pkg/CHNOSZ/NAMESPACE 2020-07-04 02:26:49 UTC (rev 539)
@@ -38,8 +38,6 @@
"label.figure", "syslab",
# equilibrium vignette
"usrfig",
-# hotspring vignette
- "strip",
# ecipex package
"count.elements",
# (no other functions are used in the tests)
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2020-07-04 02:00:42 UTC (rev 538)
+++ pkg/CHNOSZ/R/diagram.R 2020-07-04 02:26:49 UTC (rev 539)
@@ -665,78 +665,6 @@
return(invisible(out))
}
-
-strip <- function(affinity, ispecies=NULL, col=NULL, ns=NULL,
- xticks=NULL, ymin=-0.2, xpad=1, cex.names=0.7) {
- # make strip chart(s) showing the degrees of formation
- # of species as color bars of varying widths
- # extracted from bison/plot.R 20091102 jmd
- # figure out defaults
- a <- affinity
- xlim <- range(a$vals[[1]])
- xlab <- axis.label(a$vars[1])
- if(is.null(ispecies)) ispecies <- list(1:nrow(a$species))
- if(!is.list(ispecies)) ispecies <- list(ispecies)
- if(!is.null(ns) & !is.list(ns)) ns <- list(ns)
- if(is.null(col)) col <- rainbow(length(ispecies[[1]]))
- # how many strip charts on this plot
- # determined by the length of the ispecies list
- nstrip <- length(ispecies)
- # start up the plot
- plot.xlim <- c(xlim[1]-xpad,xlim[2]+xpad)
- ymax <- nstrip+0.3
- thermo.plot.new(xlim=plot.xlim,ylim=c(ymin,ymax),xlab=xlab,ylab="",
- side=c(1,3),mar=par('mar'),plot.box=FALSE)
- if(!is.null(xticks)) {
- # mark the positions of the sites on the x-axis
- for(i in 1:5) lines(rep(xticks[i],2),c(ymin,ymin+0.1),lwd=6,col=col[i])
- for(i in 1:5) lines(rep(xticks[i],2),c(ymax,ymax-0.1),lwd=6,col=col[i])
- }
- for(j in 1:nstrip) {
- # get the degrees of formation
- e <- equilibrate(a, normalize=TRUE, ispecies=ispecies[[j]])
- # depict the relative stabilities of the proteins as color bars
- # make vertical color bars sizes proportional to activities
- xs <- seq(xlim[1],xlim[2],length.out=length(e$loga.equil[[1]]))
- # total height of the stack
- ly <- 0.5
- # where to start plotting on the y-axis
- y0 <- j-ly
- for(i in 1:length(e$loga.equil[[1]])) {
- # create a vector of activities
- loga <- numeric()
- for(k in 1:length(ispecies[[j]])) loga <- c(loga, e$loga.equil[[k]][i])
- # turn the log activities into alphas
- act <- 10^loga
- alpha <- act/sum(act)
- alpha.order <- order(alpha, decreasing=FALSE)
- # plot the bars, least abundant first
- dy <- y0 # keep track of the height of the stack
- for(k in 1:length(ispecies[[j]])) {
- y1 <- dy
- y2 <- y1 + alpha[alpha.order[k]] * ly
- dy <- y2
- # note: lwd should be lowered here for very high-resolution plots
- lines(c(xs[i],xs[i]),c(y1,y2),col=col[alpha.order[k]],lwd=3,lend="butt")
- }
- }
- # label the color bar
- text(xlim[1],j-ly*1.4,names(ispecies[j]),adj=0,cex=cex.names)
- # add inset plot showing the relative numbers of species
- if(!is.null(ns)) {
- ys1 <- y0 - 0.85*ly
- ys2 <- y0 - 0.15*ly
- xmax <- xlim[2]
- xmin <- xlim[1] + 17/22*diff(xlim)
- xss <- seq(xmin,xmax,length.out=length(ispecies[[j]]))
- yss <- numeric()
- for(i in 1:length(ispecies[[j]])) yss <- c(yss,ys1+(ys2-ys1)*(ns[[j]][i]/max(ns[[j]])))
- points(xss,yss,pch=20,cex=0.5)
- lines(xss,yss)
- }
- }
-}
-
find.tp <- function(x) {
# find triple points in an matrix of integers 20120525 jmd
# these are the locations closest to the greatest number of different values
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2020-07-04 02:00:42 UTC (rev 538)
+++ pkg/CHNOSZ/inst/NEWS 2020-07-04 02:26:49 UTC (rev 539)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.3.6-13 (2020-07-04)
+CHANGES IN CHNOSZ 1.3.6-14 (2020-07-04)
--------------------------------------
- Use updated data for the 2019 Deep Earth Water (DEW) model.
@@ -32,6 +32,8 @@
- Revise obigt.Rmd to reduce the size of the HTML file and make
deep linking to individual sections work.
+- Remove vignete hotspring.Rnw.
+
- TODO: thermo.refs(): add current date, write 'thermo()$obigt'.
- TODO: use degree sign in messages from subcrt().
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2020-07-04 02:00:42 UTC (rev 538)
+++ pkg/CHNOSZ/man/diagram.Rd 2020-07-04 02:26:49 UTC (rev 539)
@@ -1,7 +1,6 @@
\encoding{UTF-8}
\name{diagram}
\alias{diagram}
-\alias{strip}
\alias{find.tp}
\title{Chemical Activity Diagrams}
\description{
@@ -34,8 +33,6 @@
main = NULL, legend.x = NA,
# plotting controls
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)
find.tp(x)
}
@@ -67,7 +64,7 @@
\item{spline.method}{character, method used in \code{\link{splinefun}}}
\item{contour.method}{character, labelling method used in \code{\link{contour}} (use NULL for no labels).}
\item{levels}{numeric, levels at which to draw contour lines}
- \item{col}{character, color of activity lines (1D diagram) or predominance field boundaries (2D diagram), or colors of bars in a strip diagram (\code{strip})}
+ \item{col}{character, color of activity lines (1D diagram) or predominance field boundaries (2D diagram)}
\item{col.names}{character, colors for labels of species}
\item{fill}{character, colors used to fill predominance fields}
\item{fill.NA}{character, color for grid points with NA values}
@@ -87,11 +84,6 @@
\item{plot.it}{logical, make a plot?}
\item{tplot}{logical, set up plot with \code{\link{thermo.plot.new}}?}
\item{affinity}{list, object returned by \code{\link{affinity}}}
- \item{ispecies}{numeric, which species to consider (default of \code{NULL} is to consider all species)}
- \item{ns}{numeric, numbers of species, used to make inset plots for strip diagrams}
- \item{xticks}{numeric, location of supplemental tick marks on x-axis}
- \item{ymin}{numeric, lower limit of y-axis}
- \item{xpad}{numeric, amount to extend x-axis on each side}
\item{x}{matrix, value of the \code{predominant} list element from \code{diagram}}
\item{...}{additional arguments passed to \code{\link{plot}} or \code{\link{barplot}}}
}
@@ -171,17 +163,6 @@
}
\section{Other Functions}{
-A different incarnation of 1-D speciation diagrams is provided by \code{strip}.
-This function generates any number of strip diagrams in a single plot.
-The diagrams are made up of colors bars whose heights represent the relative abundances of species; the color bars are arranged in order of abundance and the total height of the stack of colors bars is constant.
-If \code{ispecies} is a list, the number of strip diagrams is equal to the number of elements of the list, and the elements of this list are numeric vectors that identify the species to consider for each diagram.
-The strips are labeled with the \code{\link{names}} of \code{ispecies}.
-If \code{col} is NULL, the colors of the bars are generated using \code{\link{rainbow}}.
-Supplemental ticks can be added to the x-axis at the locations specified in \code{xtick}; they are larger than the standard ticks and have colors corresponding to those of the color bars.
-\code{ymin} can be decreased in order to add more space at the bottom of the plot, and \code{xpad} can be changed in order to increase or decrease the size of the x-axis relative to the width of the strips.
-An inset dot-and-line plot is created below each strip if \code{ns} is given.
-This argument has the same format as \code{ispecies}, and can be used e.g. to display the relative numbers of species for comparison with the stability calculations.
-
\code{find.tp} finds the locations in a matrix of integers that are surrounded by the greatest number of different values.
The function counts the unique values in a 3x3 grid around each point and returns a matrix of indices (similar to \code{\link{which}(..., arr.ind = TRUE)}) for the maximum count (ties result in more than one pair of indices).
It can be used with the output from \code{diagram} for calculations in 2 dimensions to approximately locate the triple points on the diagram.
@@ -198,7 +179,6 @@
\seealso{
Other examples are present in the help for \code{\link{protein}} and \code{\link{buffer}}, and even more can be found in \code{\link{demos}}.
-See the vignette \emph{Hot-spring proteins in CHNOSZ} for an example of the \code{strip} charts.
}
\examples{
Deleted: pkg/CHNOSZ/vignettes/hotspring.Rnw
===================================================================
--- pkg/CHNOSZ/vignettes/hotspring.Rnw 2020-07-04 02:00:42 UTC (rev 538)
+++ pkg/CHNOSZ/vignettes/hotspring.Rnw 2020-07-04 02:26:49 UTC (rev 539)
@@ -1,727 +0,0 @@
-%% LyX 2.3.1 created this file. For more info, see http://www.lyx.org/.
-%% Do not edit unless you really know what you are doing.
-\documentclass[english]{article}
-\usepackage{mathpazo}
-\renewcommand{\sfdefault}{lmss}
-\renewcommand{\ttdefault}{lmtt}
-\usepackage[T1]{fontenc}
-\usepackage[latin9]{inputenc}
-\usepackage{geometry}
-\geometry{verbose,tmargin=2.5cm,bmargin=2.5cm,lmargin=2.5cm,rmargin=2.5cm}
-\usepackage{color}
-\usepackage{babel}
-\usepackage{amsbsy}
-\usepackage{amssymb}
-\usepackage[numbers]{natbib}
-\usepackage[unicode=true,pdfusetitle,
- bookmarks=true,bookmarksnumbered=false,bookmarksopen=false,
- breaklinks=false,pdfborder={0 0 0},pdfborderstyle={},backref=false,colorlinks=true]
- {hyperref}
-\hypersetup{
- citecolor=magenta, urlcolor=blue}
-
-\makeatletter
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
-%\VignetteIndexEntry{Hot-spring proteins in CHNOSZ}
-%\VignetteEngine{knitr::knitr}
-
-% so DOIs in bibliography show up as hyperlinks
-\newcommand*{\doi}[1]{\href{https://doi.org/#1}{doi: #1}}
-
-\makeatother
-
-\begin{document}
-<<setup, include=FALSE, cache=FALSE>>=
-library(knitr)
-## set global chunk options
-opts_chunk$set(fig.path='figure/hotspring-', cache.path='cache/hotspring-', fig.align='center', fig.show='hold', par=TRUE)
-## set code/output width
-options(width=85)
-## set number of digits
-options(digits=5)
-## tune details of base graphics (http://yihui.name/knitr/hooks)
-knit_hooks$set(par=function(before, options, envir){
-if (before && options$fig.show!='none') par(mar=c(4,4,1,1),cex.lab=.95,cex.axis=.9,mgp=c(2,.7,0),tcl=-.3)
-})
-@
-\title{Hot-spring proteins in CHNOSZ}
-\author{Jeffrey M. Dick and Everett L. Shock}
-\maketitle
-
-\section{Introduction}
-
-This document is intended to demonstrate the calculations described
-in two recent papers \citep{DS11,DS13} dealing with the distribution
-and abundances of proteins in ``Bison Pool'', a hot spring in Yellowstone
-National Park. The calculations use metastable equilibrium to interrelate
-the compositions of proteins (from metagenomic data) with environmental
-conditions, particularly pH, temperature, and redox chemistry. This
-document is focused on the details of the calculations and has little
-in the way of introduction of concepts or interpretation and discussion;
-please see the papers for more details.
-
-There is no hidden code in this document; all code used to make the
-figures is shown in the blocks, but much of the text output (particularly
-in the figure-generating code) has been suppressed. The R code here
-is based on the Supporting Information for the papers, with modifications.
-In order to keep the code here as short and efficient as possible,
-some aspects of the published figures including manual labeling and
-higher resolution are not all reproduced; these changes are noted
-where possible. Also, some shortcuts (e.g. using \texttt{xtabs()}
-to create the BLAST frequency table, and \texttt{revisit()} to calculate
-the $\Delta G_{tr}$) have been introduced here in order to streamline
-the code.
-
-Load CHNOSZ and the thermodynamic database. In order to reproduce
-the calculations from the 2011 paper, we load old values of standard
-Gibbs energy and enthalpy of the methionine sidechain group from \citep{DLH06};
-these are inaccurate values that were updated by \citep{LD12} and
-are available starting in CHNOSZ\_0.9-9 (the current values are used
-further below). This step also loads values for the glycine group
-and the protein backone {[}UPBB{]} that were superseded in 2019.
-
-<<libraryCHNOSZ>>=
-library(CHNOSZ)
-reset()
-add.obigt("OldAA")
-@
-
-Some values that are shared among different calculations: measured
-temperature ($^{\circ}$C) and pH. These are representative values
-only; the actual values are not constant but vary due to water flow,
-weather, animals, etc.
-
-<<TpH>>=
-bison.T <- c(93.3, 79.4, 67.5, 65.3, 57.1)
-bison.pH <- c(7.350, 7.678, 7.933, 7.995, 8.257)
-@
-
-Now let's plot the measured temperature and pH in the hot spring.
-Distances (in meters) of the sampling sites are measured from the
-source of the hot spring. The \texttt{Tfun()} and \texttt{pHfun()}
-are also used further below.
-
-<<TpHplot, fig.width=6, fig.height=3>>=
-distance <- c(0, 6, 11, 14, 22)
-par(mfrow=c(1, 2), mar=c(4, 4, 3, 2))
-xpoints <- seq(0, 22, length.out=128)
-# T plot
-plot(distance, bison.T, xlab="distance, m", ylab=axis.label("T"))
-Tfun <- splinefun(distance, bison.T, method="mono")
-lines(xpoints, Tfun(xpoints))
-
-# pH plot
-plot(distance, bison.pH, xlab="distance, m", ylab="pH")
-pHfun <- splinefun(distance, bison.pH, method="mono")
-lines(xpoints, pHfun(xpoints))
-@
-
-\href{https://doi.org/10.1371/journal.pone.0022782.g005}{DOI link}
-to published figure (panels A and B). The length of \texttt{xpoints}
-here, 128, is the resolution used for the figures in the 2011 paper
-(defined in the \texttt{mkargs()} function of the Supporting Information).
-
-\section{General setup}
-
-Load the proteins. These are ``model proteins'', i.e. average amino
-acid compositions of sequences classified according to their functional
-annotation or by their phylum association, at each of the five sampling
-sites located in the hot spring.
-
-<<proteins>>=
-# read the amino acid compositions
-aa.annot <- read.csv(system.file("extdata/protein/DS11.csv", package="CHNOSZ"), as.is=TRUE)
-aa.phyla <- read.csv(system.file("extdata/protein/DS13.csv", package="CHNOSZ"), as.is=TRUE)
-@
-
-Here are the site names for the sampling locations (also referred
-to as sites 1\textendash 5).
-
-<<sitename>>=
-sites <- c("N", "S", "R", "Q", "P")
-sitenames <- paste("bison", sites, sep="")
-@
-
-Here are the classifications according to functional annotation:
-
-<<classes>>=
-classes <- unique(aa.annot$protein)
-classes
-@
-
-Here are names of phyla and colors and line types used for plotting
-(colors based on Figure 2 of Wu and Eisen, 2008 \citep{WE08}, with
-modifications):
-
-<<names>>=
-# the names of the phyla in alphabetical order (except Deinococcus-Thermus at end)
-phyla.abc <- sort(unique(aa.phyla$organism))[c(1:7,9:11,8)]
-# an abbreviation for Dein.-Thermus
-phyla.abbrv <- phyla.abc
-phyla.abbrv[[11]] <- "Dein.-Thermus"
-phyla.cols <- c("#f48ba5", "#f2692f", "#cfdd2a",
- "#962272", "#87c540", "#66c3a2", "#12a64a", "#f58656",
- "#ee3237", "#25b7d5", "#3953a4")
-phyla.lty <- c(1:6, 1:5)
-phyla.abbrv
-@
-
-\section{Average oxidation state of carbon}
-
-The average oxidation state of carbon ($\overline{Z}_{\mathrm{C}}$)
-is calculated using
-\begin{equation}
-\overline{Z}_{\mathrm{C}}=\frac{Z-n_{\mathrm{H}}+2\left(n_{\mathrm{O}}+n_{\mathrm{S}}\right)+3n_{\mathrm{N}}}{n_{\mathrm{C}}}\,,\label{eq:ZC}
-\end{equation}
-where $n_{\mathrm{C}}$, $n_{\mathrm{H}}$, $n_{\mathrm{N}}$, $n_{\mathrm{O}}$
-and $n_{\mathrm{S}}$ are the numbers of the indicated subscripted
-elements in the chemical formula of a protein or other chemical species,
-and $Z$ is the net charge of the chemical species.
-
-<<ZCplot, fig.width=5, fig.height=5, out.width='.49\\textwidth'>>=
-# 2011 plot
-ylab <- expression(bar(italic(Z))[C])
-plot(0, 0, xlim=c(-0.5, 5), ylim=c(-0.27, -0.11), xlab="location", xaxt="n", ylab=ylab)
-axis(1, at=1:5)
-col <- c("green", rep("black", 20))
-lwd <- c(3, rep(1, 20))
-clab <- c("hydrolase", "overall", "protease", "oxidoreductase", "transport", "membrane", "permease")
-pf.annot <- protein.formula(aa.annot)
-ZC.annot <- ZC(pf.annot)
-for(i in 1:length(classes)) {
- lines(1:5, ZC.annot[(1:5)+5*(i-1)], col=col[i], lwd=lwd[i])
- if(classes[i] %in% clab) text(0.8, ZC.annot[1+5*(i-1)], classes[i], adj=1)
-}
-title(main="annotations")
-
-# 2013 plot
-pf.phyla <- protein.formula(aa.phyla)
-ZC.phyla <- ZC(pf.phyla)
-# set up plot
-plot(0, 0, xlim=c(1, 5), ylim=c(-0.27, -0.11), xlab="location", ylab=ylab)
-for(i in 1:length(phyla.abc)) {
- # which of the model proteins correspond to this phylum
- iphy <- which(aa.phyla$organism==phyla.abc[i])
- # the locations (of 1, 2, 3, 4, 5) where this phylum is found
- ilocs <- match(aa.phyla$protein[iphy], sitenames)
- # the plotting symbol: determined by alphabetical position of the phylum
- points(ilocs, ZC.phyla[iphy], pch=i-1, cex=1.2)
- # a line to connect same phyla occurring at adjacent sites
- inlocs <- rep(NA, 5)
- inlocs[ilocs] <- ilocs
- lines(inlocs, ZC.phyla[iphy][match(1:5, ilocs)])
-}
-legend("bottomright", pch=0:10, legend=phyla.abbrv, bg="white", cex=0.9)
-title(main="major phyla")
-@
-
-DOI links for these plots as they appeared in the publications: \href{https://doi.org/10.1371/journal.pone.0022782.g004}{plot 1},
-\href{https://doi.org/10.1371/journal.pone.0072395.g001}{plot 2}.
-Compared to the papers, a different \emph{y}-axis scale is used here,
-in order to have same scale on both plots; permease at lower $\overline{Z}_{\mathrm{C}}$
-does not appear in plot 1 here.
-
-\section{Relative stabilities of ``overall'' model proteins}
-
-\subsection{Formation reactions of proteins from basis species}
-
-Function to setup the basis species. The basis species consist of
-$\mathrm{HCO_{3}^{-}}$, $\mathrm{H_{2}O}$, $\mathrm{NH_{3}}$, $\mathrm{HS^{-}}$,
-$\mathrm{H_{2}}$ and $\mathrm{H^{+}}$ (all aqueous species except
-for $\mathrm{H_{2}O}$ liquid).
-
-<<setup.basis>>=
-setup.basis <- function() {
- basis(c("HCO3-", "H2O", "NH3", "HS-", "H2", "H+"))
- basis(c("HCO3-", "NH3", "HS-", "H+"), c(-3, -4, -7, -7.933))
-}
-@
-
-Set up the basis species and the species. Here, we add the proteins
-using the previously read amino acid compositions (\texttt{aa.annot}),
-and we save their index number (\texttt{ip.annot}) for use later.
-For now, we load the species corresponding to the ``overall'' model
-proteins (the first 5 in aa.annot).
-
-<<overall>>=
-setup.basis()
-ip.annot <- add.protein(aa.annot)
-species("overall", sitenames)
-@
-
-The first 6 columns there represent the stoichiometry of the reactions
-to form the proteins from the basis species. The reactions can be
-divided by the lengths (number of amino acid residues) of the proteins
-to write per-residue formation reactions.
-
-<<residue>>=
-pl <- protein.length(ip.annot[1:5])
-mysp <- species()
-mysp[, 1:6]/pl
-@
-
-Note e.g. the higher coefficient on $\mathrm{H_{2}}$ for sites 1
-and 2; increasing the activity of $\mathrm{H_{2}}$ (more reducing
-conditions) has a relatively more favorable mass-action effect on
-the formation of the proteins at the higher temperature sites.
-
-\subsection{Chemical affinities along a chemical gradient}
-
-Function to calculate $\log a_{\mathrm{H_{2}}}$ as a linear equation
-in $T$. This was used to fit the spatial distribution of proteins
-in the 2011 paper \citep{DS11} (also shown as Equation 2 in 2013
-\citep{DS13}).
-\begin{equation}
-\log a_{\mathrm{H_{2_{\left(aq\right)}}}}=-11+3/40\times T\mathrm{\left(^{\circ}C\right)}\label{eq:logaH2}
-\end{equation}
-
-<<logaH2>>=
-get.logaH2 <- function(T) -11 + T * 3/40
-@
-
-Calculate the residue-normalized chemical affinities of the formation
-reactions of the overall model proteins. First set activities of the
-proteins equal to unity (logarithm of activity equal to zero). Then
-calculate affinities per mole of protein for the temperature, pH and
-$\log a_{\mathrm{H_{2}}_{\left(aq\right)}}$ of each site. Use the
-lengths of the model proteins to calculate the affinities per residue.
-
-<<affinity>>=
-species(1:5, 0)
-a <- affinity(T=bison.T, pH=bison.pH, H2=get.logaH2(bison.T))
-a.res <- t(as.data.frame(a$values))/pl
-a.res
-@
-
-\href{https://doi.org/10.1371/journal.pone.0022782.t005}{DOI link}
-to published table.
-
-The affinities are expressed as dimensionless values, i.e. $\boldsymbol{A}/2.303RT$
-where $\boldsymbol{A}$, $R$ and $T$ stand for chemical affinity,
-gas constant, and temperature in Kelvin. The affinities are all negative,
-but show a progression from higher for protein (row) 1 at the conditions
-of sites (columns) 1 and 2 to higher for protein (row) 4 at the conditions
-of sites (columns) 3 to 5.
-
-<<affinitymax>>=
-apply(a.res, 2, which.max)
-@
-
-To go from residue-normalized affinities to the actual progression
-of relative stabilities of the proteins (i.e. size-adjusted affinities),
-we subtract the logarithm of the protein length from the per-residue
-affinities (also see description following Eq. 19 in Ref. \citep{DS11}).
-
-<<affinityadj>>=
-a.res <- a.res - log10(pl)
-apply(a.res, 2, which.max)
-@
-
-This 1-2-4 progression of stabilities is visualized below; for some
-groups of proteins, an update to the thermodynamic properties of the
-methionine sidechain causes proteins from site 3 also to become stable.
-
-\subsection{Relative stabilities along a chemical gradient}
-
-Set the temperature limits (\texttt{Tlim}) over which to perform the
-calculations.
-
-<<Tlim>>=
-Tlim <- c(50, 100)
-@
-
-Make two plots here. (1) The metastable equilibrium predominance diagram
-of the overall model proteins as a function of temperature and logarithm
-of activity of hydrogen. The stability fields for the proteins from
-the higher temperatures are at higher activities of hydrogen in the
-diagram. The dotted line passes through the stability fields of the
-model proteins at approximately the actual environmental temperatures.
-(2) Combine the gradients of temperature, pH and hydrogen activity
-to calculate the metastable equilibrium activities of the proteins.
-The total activity of residues is set by reference activities of the
-proteins equal to $10^{-3}$. In order to label the $x$-axis ``distance'',
-modify a couple of entries in the list returned by \texttt{affinity()}
-(\texttt{vars}, \texttt{vals}); otherwise the $x$-axis would represent
-temperature (the first variable in the argument list to \texttt{affinity()}).
-
-<<TlogaH2plot, fig.width=6, fig.height=3, message=FALSE, results='hide'>>=
-par(mfrow=c(1, 2))
-# first plot
-a <- affinity(T=Tlim, H2=c(-7, -4))
-diagram(a, fill=NULL, names=as.character(1:5), normalize=TRUE)
-lines(Tlim, get.logaH2(Tlim), lty=3)
-# second plot
-species(1:5, -3)
-xT <- Tfun(xpoints)
-xpH <- pHfun(xpoints)
-xH2 <- get.logaH2(xT)
-a <- affinity(T=xT, pH=xpH, H2=xH2)
-a$vars[1] <- "distance, m"
-a$vals[[1]] <- xpoints
-e <- equilibrate(a, normalize=TRUE)
-diagram(e, legend.x=NULL)
-legend("bottom", lty=1:5, legend=1:5, bty="n", cex=0.6)
-@
-
-\href{https://doi.org/10.1371/journal.pone.0022782.g005}{DOI link}
-to published figure (panels B and F).
-
-\clearpage
-
-Plot the equilibrium degrees of formation as a function of distance
-for different classes of proteins. Calculate the affinities for all
-of the proteins. For each group of three, make strip charts using
-the \texttt{strip()} function in CHNOSZ. The heights of the bars represent
-the relative abundances of the model proteins. The five steps of the
-color code go from site 1 (red) to site 5 (blue).
-
-<<stripplot, fig.width=6, fig.height=3.5, message=FALSE, results='hide'>>=
-loadclass <- function(class) {
- species(delete=TRUE)
- species(rep(class, each=5), rep(sitenames, length(class)))
-}
-xclasses <- c("overall", "transferase", "transport", "synthetase", "membrane", "permease")
-loadclass(xclasses)
-a <- affinity(T=xT, pH=xpH, H2=xH2)
-a$vars[1] <- "distance, m"
-a$vals[[1]] <- xpoints
-col <- c("red", "orange", "yellow", "green", "blue")
-par(mfrow=c(1, 2), mar=c(4, 4, 1, 1))
-for(i in 1:2) {
- ispecies <- lapply((1:3)+(i-1)*3, function(x) {1:5+(x-1)*5} )
- names(ispecies) <- xclasses[(1:3)+(i-1)*3]
- strip(a = a, ispecies = ispecies, col = col, xticks = distance, cex.names = 1)
-}
-@
-
-\href{https://doi.org/10.1371/journal.pone.0022782.g007}{DOI link}
-to published figure.
-
-\clearpage
-
-\section{Comparing old and new methionine sidechain parameters}
-
-Make some $T-\log a_{\mathrm{H_{2}}}$ metastable equilibrium predominance
-diagrams using different values for the thermodynamic properties of
-the methionine sidechain group {[}Met{]}. The first row shows the
-results using the old values (defined on page 1); for the second row
-(j=2), we reset the database to use the current (revised) parameters.
-
-<<methionine, fig.width=6, fig.height=4, message=FALSE, results='hide', cache=TRUE>>=
-par(mfrow=c(2, 3))
-for(j in 1:2) {
- # use old [Met] for first row and new [Met] for second row
- if(j==2) {
- reset()
- add.obigt("OldAA", c("[Gly]", "[UPBB]"))
- ip.annot <- add.protein(aa.annot)
- }
- # setup basis species and proteins
- setup.basis()
- # make the plots
- for(annot in c("overall", "transferase", "synthase")) {
- ip <- ip.annot[aa.annot$protein==annot]
- a <- affinity(T=c(50, 100), H2=c(-7, -4), iprotein=ip)
- diagram(a, fill=NULL, names=as.character(1:5), normalize=TRUE)
- # add logaH2-T line
- lines(par("usr")[1:2], get.logaH2(par("usr")[1:2]), lty=3)
- # add a title
- title(main=annot)
- }
-}
-@
-
-\href{https://doi.org/10.1371/journal.pone.0072395.g002}{DOI link}
-to published figure, which was made with higher resolution (256) and
-manual placement of some labels.
-
-\section{Relative abundance calculations}
-
-Function to return the fractional abundances based on BLAST counts,
-stored in the \texttt{ref} column of \texttt{aa.phyla}.
-
-<<alpha.blast>>=
-alpha.blast <- function() {
- out <- xtabs(ref ~ protein + organism, aa.phyla)
- # put it in correct order, then turn counts into fractions
- out <- out[c(1,5:2), c(1:7,9:11,8)]
- out <- out/rowSums(out)
- return(out)
-}
-@
-
-Function to calculate metastable equilibrium degrees of formation
-of proteins (normalized to residues) as a function of $\log a_{\mathrm{H_{2}}}$
-for a specified location. Note: in the call to \texttt{affinity()},
-increase the resolution from 101 to 1001 to reproduce the calculations
-in paper.
-
-<<alpha.equil>>=
-alpha.equil <- function(i=1) {
- # order the names and counts to go with the alphabetical phylum list
- iloc <- which(aa.phyla$protein==sitenames[i])
- iloc <- iloc[order(match(aa.phyla$organism[iloc], phyla.abc))]
- # set up basis species, with pH specific for this location
- setup.basis()
- basis("pH", bison.pH[i])
- # calculate metastable equilibrium activities of the residues
- a <- affinity(H2=c(-11, -1, 101), T=bison.T[i], iprotein=ip.phyla[iloc])
- e <- equilibrate(a, loga.balance=0, as.residue=TRUE)
- # remove the logarithms to get relative abundances
- a.residue <- 10^sapply(e$loga.equil, c)
- colnames(a.residue) <- aa.phyla$organism[iloc]
- # the BLAST profile
- a.blast <- alpha.blast()
- # calculate Gibbs energy of transformation (DGtr) and find optimal logaH2
- iblast <- match(colnames(a.residue), colnames(a.blast))
- r <- revisit(e, "DGtr", log10(a.blast[i, iblast]), plot.it=FALSE)
- # return the calculated activities, logaH2 range, DGtr values, and optimal logaH2
- return(list(alpha=a.residue, H2vals=a$vals[[1]], DGtr=r$H, logaH2.opt=r$xopt))
-}
-@
-
-Now we're ready to make a plot! We start by adding the proteins with
-amino acid composition that were read from the file. Then plot the
-degrees of formation as a function of $\log a_{\mathrm{H_{2}}}$ at
-sites 1, 3 and 5, but record the results (including calculated relative
-abundances and optimal values of $\log a_{\mathrm{H_{2}}}$) at all
-5 sites.
-
-<<alphaplot, fig.width=8, fig.height=6, tidy=FALSE, message=FALSE, results='hide', cache=TRUE>>=
-ip.phyla <- add.protein(aa.phyla)
-layout(matrix(1:6, ncol=3), heights=c(2, 1))
-equil.results <- list()
-for(i in 1:5) {
- # get the equilibrium degrees of formation and the optimal logaH2
- ae <- alpha.equil(i)
- equil.results[[i]] <- ae
- if(i %in% c(1, 3, 5)) {
- iphy <- match(colnames(ae$alpha), phyla.abc)
- # top row: equilibrium degrees of formation
- thermo.plot.new(xlim=range(ae$H2vals), ylim=c(0, 0.5), xlab=axis.label("H2"),
- ylab=expression(alpha[equil]), yline=2, cex.axis=1, mgp=c(1.8, 0.3, 0))
- for(j in 1:ncol(ae$alpha)) {
- lines(ae$H2vals, ae$alpha[, j], lty=phyla.lty[iphy[j]])
- ix <- seq(1, length(ae$H2vals), length.out=11)
- ix <- head(tail(ix, -1), -1)
- points(ae$H2vals[ix], ae$alpha[, j][ix], pch=iphy[j]-1)
- }
- title(main=paste("site", i))
- legend("topleft", pch=iphy-1, lty=phyla.lty[iphy], legend=phyla.abbrv[iphy], bg="white")
- # bottom row: Gibbs energy of transformation and position of minimum
- thermo.plot.new(xlim=range(ae$H2vals), ylim=c(0, 1/log(10)), xlab=axis.label("H2"),
- ylab=expr.property("DGtr/2.303RT"), yline=2, cex.axis=1, mgp=c(1.8, 0.3, 0))
- lines(ae$H2vals, ae$DGtr)
- abline(v=ae$logaH2.opt, lty=2)
- abline(v=get.logaH2(bison.T[i]), lty=3, lwd=1.5)
- if(i==1) legend("bottomleft", lty=c(3, 2), lwd=c(1.5, 1), bg="white",
- legend=c("Equation 2", "optimal"))
- }
-}
-@
-
-\href{https://doi.org/10.1371/journal.pone.0072395.g003}{DOI link}
-to published version of the figure.
-
-\clearpage
-
-\section{Activity of hydrogen comparison}
-
-Let's compare the computed activities of hydrogen with various redox
-indicators measured in the field.
-
-\subsection{Conversion of field measurements}
-
-\paragraph{ORP}
-
-Oxidation-reduction potential (ORP), measured in the field using a
-portable probe and pH/voltmeter, can be converted to Eh by adding
-the potential of the reference electrode, in this case silver-silver
-chloride (Ag/AgCl) in saturated KCl. As an approximation, the following
-equation from Ref. \citep{BPJ85} for Ag/AgCl (1M KCl) is used, with
-temperature in $^{\circ}$C and potential in volts.
-
-<<E.AgAgCl>>=
-E.AgAgCl <- function(T) {
- 0.23737 - 5.3783e-4 * T - 2.3728e-6 * T^2 - 2.2671e-9 * (T+273)
-}
-@
-
-Data from Bison Pool and another hot spring shown in Fig. 9 of Ref.
-\citep{DS11} are both included, but for simplicity here they are
-lumped together. The ORP values have units of mV.
-
-<<meters>>=
-T.ORP <- c(93.9, 87.7, 75.7, 70.1, 66.4, 66.2)
-pH.ORP <- c(8.28, 8.31, 7.82, 7.96, 8.76, 8.06)
-ORP <- c(-258, -227, -55, -58, -98, -41)
-@
-
-Convert ORP to effective values of $\log a_{\mathrm{H_{2}}_{\left(aq\right)}}$.
-First convert to Eh (in volts). Then convert Eh to pe. Then combine
-Eh and pH with the law of mass action for
-\[
-2e^{-}+2\mathrm{H^{+}}\rightleftharpoons\mathrm{H_{2}}_{\left(aq\right)}
-\]
-to calculate $\log a_{\mathrm{H_{2}}_{\left(aq\right)}}$. The ``law
-of mass action'' is the equality between the equilibrium constant
-($K$) and the activity product ($Q$) of the species in the reaction.
-
-<<ORP2Eh, message=FALSE, results='hide'>>=
-Eh <- ORP/1000 + E.AgAgCl(T.ORP)
-pe <- convert(Eh, "pe", T=convert(T.ORP, "K"))
-logK.ORP <- subcrt(c("e-", "H+", "H2"), c(-2, -2, 1), T=T.ORP)$out$logK
-logaH2.ORP <- logK.ORP - 2*pe - 2*pH.ORP
-@
-
-\paragraph{Sulfide/Sulfate}
-
-For sulfide/sulfate, assign activities that are equal to concentrations
-(in molal units) measured in the field season of 2005, and use the
-law of mass action for
-\[
-\mathrm{HS^{-}}+4\mathrm{H_{2}O}\rightleftharpoons\mathrm{SO_{4}}^{-2}+\mathrm{H^{+}}+4\mathrm{H_{2}}_{\left(aq\right)}
-\]
-to calculate an effective activity of hydrogen. Sulfide was determined
-spectrophotometrically at the hot spring, and sulfate was determined
-using ion chromatography on water samples returned from the field.
-
-<<sulfur, message=FALSE, results='hide'>>=
-loga.HS <- log10(c(4.77e-6, 2.03e-6, 3.12e-7, 4.68e-7, 2.18e-7))
-loga.SO4 <- log10(c(2.10e-4, 2.03e-4, 1.98e-4, 2.01e-4, 1.89e-4))
-logK.S <- subcrt(c("HS-", "H2O", "SO4-2", "H+", "H2"), c(-1, -4, 1, 1, 4), T=bison.T)$out$logK
-logaH2.S <- (logK.S + bison.pH - loga.SO4 + loga.HS) / 4
-@
-
-\paragraph{Dissolved Oxygen}
-
-Calculate the effective activity of hydrogen corresponding to the
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 539
More information about the CHNOSZ-commits
mailing list