[CHNOSZ-commits] r164 - in pkg/CHNOSZ: . R data demo inst man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Feb 19 13:59:01 CET 2017
Author: jedick
Date: 2017-02-19 13:59:00 +0100 (Sun, 19 Feb 2017)
New Revision: 164
Added:
pkg/CHNOSZ/demo/Shh.R
Removed:
pkg/CHNOSZ/inst/doc/
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/examples.R
pkg/CHNOSZ/data/protein.csv
pkg/CHNOSZ/demo/00Index
pkg/CHNOSZ/demo/protbuff.R
pkg/CHNOSZ/demo/revisit.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/examples.Rd
pkg/CHNOSZ/man/protein.Rd
pkg/CHNOSZ/tests/testthat/test-util.program.R
pkg/CHNOSZ/vignettes/anintro.Rmd
pkg/CHNOSZ/vignettes/vig.bib
Log:
anintro.Rmd: add an affinity baseline
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/DESCRIPTION 2017-02-19 12:59:00 UTC (rev 164)
@@ -1,6 +1,6 @@
Date: 2017-02-19
Package: CHNOSZ
-Version: 1.0.8-53
+Version: 1.0.8-54
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-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/R/diagram.R 2017-02-19 12:59:00 UTC (rev 164)
@@ -174,6 +174,7 @@
col.names <- rep(col.names, length.out=ngroups)
## make up some names for lines/fields if they are missing
+ is.pname <- FALSE
if(missing(names)) {
# properties of basis species or reactions?
if(eout$property %in% c("G.basis", "logact.basis")) names <- rownames(eout$basis)
@@ -185,6 +186,7 @@
else names <- as.character(eout$species$name)
# remove non-unique organism or protein names
if(all(grepl("_", names))) {
+ is.pname <- TRUE
# everything before the underscore (the protein)
pname <- gsub("_.*$", "", names)
# everything after the underscore (the organism)
@@ -201,13 +203,10 @@
}
## apply formatting to chemical formulas 20170204
- if(format.names) {
+ if(all(grepl("_", names))) is.pname <- TRUE
+ if(format.names & !is.pname) {
exprnames <- as.expression(names)
- for(i in seq_along(exprnames)) {
- # can the name be parsed as a chemical formula?
- mtry <- suppressWarnings(try(makeup(exprnames[[i]]), TRUE))
- if(!identical(class(mtry), "try-error")) exprnames[[i]] <- expr.species(exprnames[[i]])
- }
+ for(i in seq_along(exprnames)) exprnames[[i]] <- expr.species(exprnames[[i]])
names <- exprnames
}
Modified: pkg/CHNOSZ/R/examples.R
===================================================================
--- pkg/CHNOSZ/R/examples.R 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/R/examples.R 2017-02-19 12:59:00 UTC (rev 164)
@@ -31,7 +31,7 @@
demos <- function(which=c("sources", "protein.equil", "add.obigt", "affinity", "NaCl", "density",
"ORP", "revisit", "findit", "ionize", "buffer", "protbuff", "yeastgfp", "mosaic",
- "copper", "solubility", "wjd", "dehydration", "bugstab", "activity_ratios"), to.file=FALSE) {
+ "copper", "solubility", "wjd", "dehydration", "bugstab", "Shh", "activity_ratios"), to.file=FALSE) {
# run one or more demos from CHNOSZ with ask=FALSE, and return the value of the last one
for(i in 1:length(which)) {
# say something so the user sees where we are
Modified: pkg/CHNOSZ/data/protein.csv
===================================================================
--- pkg/CHNOSZ/data/protein.csv 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/data/protein.csv 2017-02-19 12:59:00 UTC (rev 164)
@@ -452,3 +452,20 @@
RBL1,HALNC,BBA+03,O85040,1,45,6,32,31,23,46,14,25,25,35,11,16,21,12,28,16,30,31,10,16
RBL2,HALNC,BBA+03,Q9ZHZ4,1,47,2,33,26,23,51,14,31,28,30,19,14,23,14,20,20,20,18,5,21
RBS,HALNC,BBA+03,P45686,1,11,1,3,8,4,5,2,5,6,7,5,7,5,8,4,6,2,8,3,10
+SHH,HUMAN,UniProt,Q15465.N,1,14,3,10,14,5,16,6,7,15,11,3,8,7,3,13,12,8,9,3,7
+OLIG2,HUMAN,UniProt,Q13516,1,51,5,10,10,5,35,18,7,14,29,11,4,31,5,13,50,11,10,1,3
+NKX22,HUMAN,UniProt,O95096,1,29,1,14,18,7,18,6,4,16,23,4,8,25,19,17,24,15,11,4,10
+FOXA2,HUMAN,UniProt,Q9Y261,1,55,4,9,18,9,53,23,8,19,31,30,19,42,24,13,55,12,8,4,21
+IRX3,HUMAN,UniProt,P78415,1,80,2,23,43,11,47,9,8,18,51,4,13,62,13,25,44,17,13,5,13
+PAX6,HUMAN,UniProt,P26367,1,19,6,14,21,9,32,6,19,13,25,14,21,37,30,33,51,31,23,6,12
+NKX62,HUMAN,UniProt,Q9C056,1,35,1,15,11,9,35,5,5,21,30,5,7,26,9,16,16,9,12,4,6
+DBX1,HUMAN,UniProt,A6NMT0,1,23,1,12,22,17,27,6,8,18,31,6,5,46,14,24,46,17,11,3,6
+DBX2,HUMAN,UniProt,Q6ZNG2,1,40,8,10,19,11,22,3,8,17,35,4,10,39,17,24,37,10,14,6,5
+NKX61,HUMAN,UniProt,P78426,1,52,1,11,17,8,29,8,6,21,35,8,6,44,12,15,60,16,7,4,7
+PAX7,HUMAN,UniProt,P23759,1,37,8,23,29,14,43,14,22,25,34,14,14,50,27,35,55,29,26,5,16
+GLI1,HUMAN,UniProt,P08151,1,65,31,39,63,24,117,40,17,33,89,21,35,149,48,66,128,58,44,5,34
+GLI2,HUMAN,UniProt,P10070,1,148,27,64,71,28,140,68,33,48,135,40,50,174,94,80,188,81,75,4,38
+GLI3,HUMAN,UniProt,P10071,1,112,29,71,75,34,122,67,52,58,121,42,63,152,88,83,215,90,63,4,39
+PTC1,HUMAN,UniProt,Q13635,1,116,27,63,78,64,101,42,62,50,153,31,49,102,57,79,105,83,108,21,56
+SMO,HUMAN,UniProt,Q99835,1,70,26,28,36,37,57,15,33,29,71,14,25,61,26,47,52,43,53,19,18
+GLI3R,HUMAN,UniProt,P10071,1,10,0,6,4,5,5,12,7,0,7,5,3,26,0,7,17,6,3,0,8
Modified: pkg/CHNOSZ/demo/00Index
===================================================================
--- pkg/CHNOSZ/demo/00Index 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/demo/00Index 2017-02-19 12:59:00 UTC (rev 164)
@@ -17,4 +17,5 @@
wjd Gibbs energy minimization: prebiological atmospheres and cell periphery of yeast
dehydration log K of dehydration reactions; SVG file contains tooltips and links
bugstab formation potential of microbial proteins in colorectal cancer
+Shh affinities of transcription factors relative to Sonic hedgehog
activity_ratios mineral stability plots with activity ratios on the axes
Added: pkg/CHNOSZ/demo/Shh.R
===================================================================
--- pkg/CHNOSZ/demo/Shh.R (rev 0)
+++ pkg/CHNOSZ/demo/Shh.R 2017-02-19 12:59:00 UTC (rev 164)
@@ -0,0 +1,197 @@
+# Compare affinities of Sonic hedgehog and transcription factors involved in dorsal-ventral patterning
+# (Dick, 2015. Chemical integration of proteins in signaling and development. http://dx.doi.org/10.1101/015826)
+
+# UniProt names of the proteins
+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
+fill <- c(SHH = "#c8c7c8", OLIG2 = "#f9a330", NKX22 = "#ef2b2c", FOXA2 = "#6ab0b9",
+ NKX61 = "#b76775", DBX2 = "#35bcba", PAX7 = "#f6ef42",
+ PAX6 = "#4d2a59", IRX3 = "#63c54e", NKX62 = "#f24e33", DBX1 = "#d4e94e",
+ PTC1 = "#c7b0ee", SMO = "#8fd4ef", GLI1 = "#fcdc7e", GLI2 = "#c7e3b0", GLI3 = "#fcdc7e",
+ GLI3R = "#f1b1ae")
+
+# 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",
+ PTC1 = "Ptch1", SMO = "Smo", GLI1 = "Gli1A", GLI2 = "Gli2A", GLI3 = "Gli3A",
+ GLI3R = "Gli3R")
+
+# protein indices of Shh and the transcription factors
+ip <- match(pname, thermo$protein$protein)
+aa <- thermo$protein[ip, ]
+
+# set up basis species
+basis("CHNOS")
+basis("NH3", -7)
+
+# save as PDF file?
+pdf <- FALSE
+# draw interpretive legend?
+interp <- TRUE
+
+# 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
+pl <- protein.length(ip)
+names(A$values) <- pname
+for(i in 1:length(A$values)) A$values[[i]] <- A$values[[i]] / pl[i]
+A.SSH <- A$values$SHH
+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
+if(interp) {
+ 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")
+ 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))
+}
+for(i in 1:length(pl)) {
+ lty <- 3
+ lwd <- 1
+ # highlight SHH with solid line
+ if(pname[i]=="SHH") {
+ lty <- 1
+ lwd <- 3
+ }
+ lines(H2O, A$values[[i]], lty=lty, lwd=lwd)
+}
+# 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)
+# 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)
+# Nkx2.2
+iNKX22 <- A$NKX22 > A$FOXA2 & A$NKX22 > A$SHH
+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)
+# 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)
+# Dbx2
+iDBX2 <- A$DBX2 > A$NKX61 & A$DBX2 < A$IRX3
+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)
+# Nkx6.2
+iNKX62 <- A$NKX62 > A$IRX3 & A$NKX62 > A$DBX1
+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)
+# Shh
+#iSHH <- A$SHH > A$DBX1
+#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)
+} else {
+ text(-1.5, 0.5, "Nkx2.2")
+ text(-0.1, 0.3, "Pax6")
+ text(-0.2, 3.8, "Olig2")
+ text(-0.75, 2.45, "Irx3")
+ text(-1.13, 1.3, "Nkx6.1")
+ text(-1.5, 1, "Dbx2")
+ text(-2.4, 1.3, "Nkx6.2")
+ text(-3.8, 0.3, "Dbx1")
+}
+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?
+if(interp) {
+ # 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
+ 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
+ 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
+ xNKX22 <- H2O[max(which(A$NKX22 > A$SHH))]
+ 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)
+ 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
+ xtop <- H2O[2]
+ ytop <- 4.7
+ wtop <- 5
+ # 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
+ 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
+ iIRX3 <- min(which(A$IRX3 > A$OLIG2))
+ 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)
+ 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
+ arrows(-2.7, 2, -2.7, 2.5, 0.1)
+ 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.37, -1.55, "FP")
+ text(-1.6, -1.55, "p3")
+ text(-0.53, -1.55, "pMN")
+ text(0.2, -1.55, "p2")
+ text(0.2, 5.25, "pMN")
+ 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)
+} else {
+ # 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)
+}
+# all done!
+if(pdf) dev.off()
Modified: pkg/CHNOSZ/demo/protbuff.R
===================================================================
--- pkg/CHNOSZ/demo/protbuff.R 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/demo/protbuff.R 2017-02-19 12:59:00 UTC (rev 164)
@@ -14,9 +14,9 @@
a <- affinity(pH=c(4, 10, 300), T=c(40, 60, 300))
e <- equilibrate(a, normalize=TRUE)
fill <- ZC.col(ZC(protein.formula(species()$name)))
-diagram(e, fill=fill, format.names=FALSE)
+diagram(e, fill=fill)
title(main="Thiol peroxidases from bacteria")
-legend("topleft", describe.basis(thermo$basis[-6,]), bg="grey80", box.lwd=0)
+legend("topleft", describe.basis(thermo$basis[-6,]), bg="slategray1", box.lwd=0)
## Buffer + ionization: relative stabilities
## of E. coli sigma factors on a T-pH diagram
@@ -36,10 +36,10 @@
species(paste(proteins, "ECOLI", sep="_"))
a <- affinity(pH=c(5, 10, 300), T=c(10, 40, 300))
fill <- ZC.col(ZC(protein.formula(species()$name)))
-diagram(a, normalize=FALSE, fill=fill, format.names=FALSE)
+diagram(a, normalize=FALSE, fill=fill)
title(main=expression("Sigma factors in"~italic("E. coli")))
ptext <- c(describe.property("T", 25),
describe.basis(ibasis=c(2, 6), oneline=TRUE))
-legend("topleft", legend=c("preset conditions:", ptext), bg="grey80", box.lwd=0)
+legend("topleft", legend=c("preset conditions:", ptext), bg="slategray1", box.lwd=0)
btext <- describe.basis(ibasis=c(1, 3, 4, 5), oneline=TRUE)
-legend("bottomright", legend=c("buffered conditions:", btext), bg="grey80", box.lwd=0)
+legend("bottomright", legend=c("buffered conditions:", btext), bg="slategray1", box.lwd=0)
Modified: pkg/CHNOSZ/demo/revisit.R
===================================================================
--- pkg/CHNOSZ/demo/revisit.R 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/demo/revisit.R 2017-02-19 12:59:00 UTC (rev 164)
@@ -29,7 +29,7 @@
# calculate affinities in logfO2-logaH2O space
a <- affinity(O2=c(-85, -65), H2O=c(-5, 5), iprotein=ip)
# show the predominances
-diagram(a, normalize=TRUE, fill="heat", format.names=FALSE)
+diagram(a, normalize=TRUE, fill="heat")
# calculate the equilibrium activities
e <- equilibrate(a, loga.balance=0, normalize=TRUE)
# show the coefficient of variation
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/inst/NEWS 2017-02-19 12:59:00 UTC (rev 164)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.0.8-53 (2017-02-19)
+CHANGES IN CHNOSZ 1.0.8-54 (2017-02-19)
---------------------------------------
DOCUMENTATION:
@@ -11,7 +11,8 @@
- New demos: bugstab.R (potential diagrams for microbial proteins in
colorectal cancer), activity_ratios.R (mineral stability diagrams
- with activity ratios on the axes).
+ with activity ratios on the axes), Shh.R (affinities of
+ transcription factors relative to Sonic hedgehog).
- Converted demos: move examples from help pages to protein.equil.R,
add.obigt.R, and affinity.R.
Modified: pkg/CHNOSZ/man/examples.Rd
===================================================================
--- pkg/CHNOSZ/man/examples.Rd 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/man/examples.Rd 2017-02-19 12:59:00 UTC (rev 164)
@@ -16,7 +16,7 @@
examples(do.png = FALSE)
demos(which = c("sources", "protein.equil", "add.obigt", "affinity", "NaCl", "density",
"ORP", "revisit", "findit", "ionize", "buffer", "protbuff", "yeastgfp", "mosaic",
- "copper", "solubility", "wjd", "dehydration", "bugstab", "activity_ratios"),
+ "copper", "solubility", "wjd", "dehydration", "bugstab", "Shh", "activity_ratios"),
to.file=FALSE)
}
@@ -51,6 +51,7 @@
\code{wjd} \tab Gibbs energy minimization: prebiological atmospheres and cell periphery of yeast \cr
\code{dehydration} \tab log K of dehydration reactions; SVG file contains tooltips and links \cr
\code{bugstab} \tab formation potential of microbial proteins in colorectal cancer (Dick, 2016) \cr
+ \code{Shh} \tab affinities of transcription factors relative to Sonic hedgehog \cr
\code{activity_ratios} \tab mineral stability plots with activity ratios on the axes \cr
}
Modified: pkg/CHNOSZ/man/protein.Rd
===================================================================
--- pkg/CHNOSZ/man/protein.Rd 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/man/protein.Rd 2017-02-19 12:59:00 UTC (rev 164)
@@ -87,7 +87,7 @@
col <- rep("red", length(prot))
col[ZC > 0] <- "blue"
e <- equilibrate(a, normalize=TRUE)
-d <- diagram(e, col=col, format.names=FALSE)
+d <- diagram(e, col=col)
title(main="Bovine proteins, GSH/GSSG redox buffer")
}
Modified: pkg/CHNOSZ/tests/testthat/test-util.program.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-util.program.R 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/tests/testthat/test-util.program.R 2017-02-19 12:59:00 UTC (rev 164)
@@ -17,10 +17,10 @@
# CHNOSZ no longer has a large FASTA file to test this with 20170205
#ff <- system.file("extdata/fasta/HTCC1062.faa.xz", package="CHNOSZ")
#expect_message(aa <- read.fasta(ff), "read.fasta running 1354 calculations")
+ basis("CHNOS")
ip <- 1:nrow(thermo$protein)
- basis("CHNOS")
- expect_message(a <- affinity(iprotein=rep(ip, 3)), "affinity running 1359 calculations")
- expect_message(e <- equilibrate(a, normalize=TRUE), "equil.boltzmann running 1359 calculations")
+ expect_message(a <- affinity(iprotein=rep(ip, 3)), "affinity running .* calculations")
+ expect_message(e <- equilibrate(a, normalize=TRUE), "equil.boltzmann running .* calculations")
# test reaction method
species(c("CO2", "acetic acid"))
a <- affinity(O2=c(-90, -60, 1000))
Modified: pkg/CHNOSZ/vignettes/anintro.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/anintro.Rmd 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/vignettes/anintro.Rmd 2017-02-19 12:59:00 UTC (rev 164)
@@ -1226,9 +1226,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, col = 1:5, lwd = 2)
+diagram(e, ylim = c(-5, -2), col = 1:5, lwd = 2)
e <- equilibrate(a, normalize = TRUE)
-diagram(e, ylim = c(-5, -2.5), format.names = FALSE, col = 1:5, lwd = 2)
+diagram(e, ylim = c(-5, -2.5), col = 1:5, lwd = 2)
abline(h = logabundance, lty = 1:5, col = 1:5)
revisit(e, "DGinf", logabundance)
```
@@ -1253,6 +1253,75 @@
The minimum free energy difference occurs near `r logfO2` = -78.
This agrees with the assessment shown in Figure 4 of @Dic09 (note that the old parameters for the methionine sidechain group were used in that study).
+## An affinity baseline
+
+Because affinities of proteins often vary strongly with oxygen fugacity and other variables, it can be helpful to express the values as differences from a baseline.
+The following example compares the affinities of transcription factors involved in embryonic dorsal-ventral patterning with that of Sonic hedgehog (Shh) as a function of `r logfO2` and log*a*<sub>`r h2o`</sub> [@Dic15].
+We first list the UniProt names of Shh and 10 transcription factors, and get the `iprotein` indices (rownumbers of `thermo$protein`):
+```{r Shh_pname}
+pname <- c("SHH", "OLIG2", "NKX22", "FOXA2", "IRX3",
+ "PAX6", "NKX62", "DBX1", "DBX2", "NKX61", "PAX7")
+ip <- match(pname, thermo$protein$protein)
+```
+
+Next, set up the basis species:
+```{r Shh_basis, results="hide"}
+basis("CHNOS")
+basis("NH3", -7)
+```
+
+Now calculate the affinities of formation of the proteins from the basis species.
+The values chosen for `r logfO2` and log*a*<sub>`r h2o`</sub> covary, so we are using the transect mode of <span style="color:green">`affinity()`</span>:
+```{r Shh_affinity, message=FALSE}
+O2 <- seq(-70, -100, length.out = 10)
+H2O <- seq(0.5, -4.5, length.out = 10)
+a <- affinity(H2O = H2O, O2 = O2, iprotein = ip)
+```
+
+We would like to compare the affinities of the proteins on a per-residue scale.
+R's `lapply()` could be used here, but to show the operation more clearly we use a `for()` loop:
+```{r Shh_residue}
+pl <- protein.length(ip)
+for(i in seq_along(a$values)) a$values[[i]] <- a$values[[i]] / pl[i]
+```
+
+Then, we calculate the relative affinities with Shh as the baseline:
+```{r Shh_minusShh}
+a.Shh <- a$values[[1]]
+for(i in 1:length(a$values)) a$values[[i]] <- a$values[[i]] - a.Shh
+```
+
+Next we use <span style="color:green">`diagram()`</span> to plot the affinities.
+We set `balance = 1` to plot the values as they are---without that, <span style="color:green">`diagram()`</span> divides the values by protein length, which we have already done!
+```{marginfigure}
+See <span style="color:blue">`demo(Shh)`</span> for a plot with more interpretive labels and comments.
+```
+For this plot, we highlight and label the proteins with the highest relative affinity at some combination of `logfO2` and log*a*<sub>`r h2o`</sub> along the transect.
+Those proteins are Olig2, Irx3, Nkx6.2, Dbx1, and Shh (numbers 2, 5, 7, 8, 1 in the set we have identified).
+The last few lines are used to set up the second (upper) *x* axis, using a label generated with <span style="color:green">`axis.label()`</span>:
+
+```{r Shh_diagram, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, dpi=dpi, out.width="100%", echo=FALSE, results="hide", message=FALSE, fig.cap="Per-residue affinities of formation of transcription factors relative to Shh.", cache=TRUE, pngquant=pngquant, timeit=timeit}
+# line type, width, and color
+twc <- lapply(c(3, 1, 1), rep, length(pname))
+ihigh <- c(2, 5, 7, 8, 1)
+twc[[1]][ihigh] <- 1
+twc[[2]][ihigh] <- 3
+col <- c("#f9a330", "#63c54e", "#f24e33", "#d4e94e", "#0f0f0f")
+twc[[3]][ihigh] <- col
+diagram(a, balance = 1, ylim = c(-0.5, 4.5), xlim = c(0.5, -4.5),
+ lty = twc[[1]], lwd=twc[[2]], col = twc[[3]], names = NA)
+legend("topright", legend = c("Olig2", "Irx3", "Nkx6.2", "Dbx1", "Shh"),
+ lwd = 2, col = col)
+par(usr = c(-70, -100, -0.5, 4.5), tcl = -0.3)
+axis(3, at = seq(-70, -100, by = -5))
+mtext(axis.label("O2"), line = 1.2)
+```
+```{r Shh_diagram, eval=FALSE}
+```
+
+Among these proteins, Olig2 and Shh have the greatest affinities for formation at the extremes of oxidation and hydration state.
+This thermodynamic description suggests links between the chemical compositions of the proteins and their order of appearance along with chemical changes in developing embryos (Dick, 2015).
+
## Adding proteins
In the Rubisco example above, we saw the use of <span style="color:green">`read.fasta()`</span> to read amino acid sequences from a FASTA file.
@@ -1391,7 +1460,7 @@
res <- 300
a <- affinity(T = c(T, res), H2 = c(-8, -3, res), iprotein = ip)
fill <- ZC.col(ZC(protein.formula(ip)))
-diagram(a, normalize = TRUE, fill = fill, names = 1:5, format.names = FALSE)
+diagram(a, normalize = TRUE, fill = fill, names = 1:5)
T <- c(93.3, 79.4, 67.5, 65.3, 57.1)
logaH2 <- c(-3.38, -4.14, -5.66, -7.47, -10.02)
lines(T, logaH2, lty = 2, lwd = 2)
@@ -1442,7 +1511,7 @@
```
Use the following commands to set the basis species, add the variously protonated ATP species, calculate the affinities of the formation reactions, equilibrate the system, 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 M and *T* = 100 °C:
+This is similar to Figure 1.3 of Alberty (2003), but is calculated for *I* = 0 M 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.
Modified: pkg/CHNOSZ/vignettes/vig.bib
===================================================================
--- pkg/CHNOSZ/vignettes/vig.bib 2017-02-19 05:30:45 UTC (rev 163)
+++ pkg/CHNOSZ/vignettes/vig.bib 2017-02-19 12:59:00 UTC (rev 164)
@@ -96,6 +96,14 @@
doi = {10.1098/rsif.2013.1095},
}
+ at Article{Dic15,
+ author = {Dick, Jeffrey M.},
+ journal = {bioRxiv},
+ title = {{C}hemical integration of proteins in signaling and development},
+ year = {2015},
+ doi = {10.1101/015826},
+}
+
@Article{DLH06,
author = {Dick, Jeffrey M. and LaRowe, Douglas E. and Helgeson, Harold C.},
journal = {Biogeosciences},
More information about the CHNOSZ-commits
mailing list