### R code from vignette source 'formation.Rnw'
### Encoding: ISO8859-1

###################################################
### code chunk number 1: libraryCHNOSZ
###################################################

### loading CHNOSZ

library(CHNOSZ)

# end Sweave chunk



###################################################
### code chunk number 2: setup_four
###################################################

setup_four <- function() {

  prot <<- c("LYSC","CYC","RNAS1","MYG")

  org <<- c("CHICK","BOVIN","BOVIN","PHYCA")

  proteins <<- po <- paste(prot,org,sep="_")

  basis("CHNOS+")

  species(po)

  pl <- protein.length(po)

  loga.prot <- log10(1/pl)

  species(1:4,loga.prot)

  mymar <<- c(3.5,4,0.5,2)

}

setup_four()



###################################################
### code chunk number 3: table (eval = FALSE)
###################################################
## 
## # get values for table in presentation
## 
## # properties of formation reactions of residues
## 
## prot <- c("LYSC","CYC","RNAS1","MYG")
## 
## org <- c("CHICK","BOVIN","BOVIN","PHYCA")
## 
## po <- paste(prot,org,sep="_")
## 
## basis("CHNOS+")
## 
## #basis("CHNOS")
## 
## species(po)
## 
## pl <- protein.length(po)
## 
## # equal-activity reference state
## 
## # (unit activity of residues)
## 
## loga.prot <- log10(1/pl)
## 
## species(1:4,loga.prot)
## 
## A.in <- affinity()
## 
## # we are dealing with the ionized proteins
## 
## ionize()
## 
## A <- affinity()
## 
## # the standard state properties of reactions
## 
## logK <- affinity(property="logK")
## 
## logK <- as.numeric(ionize(A$values,logK$values))
## 
## # the activity products
## 
## logQ.act <- affinity(property="logQ")
## 
## logQ.act <- as.numeric(ionize(A$values,logQ.act$values))
## 
## # done with equal-activity reference state
## 
##   # first get to per residue
## 
## A.act <- as.numeric(A.in$values)[1:4]/pl
## 
##   # then account for the protein-residue activity changes
## 
## A.act <- -log10(pl) - (pl-1)/pl*loga.prot + A.act
## 
## logK <- logK[1:4]/pl
## 
## logQ.act <- logQ.act[1:4]/pl
## 
## logQ.act <- log10(pl) + (pl-1)/pl*loga.prot + logQ.act
## 
## ionize(FALSE)
## 
## # go to equal-affinity reference state
## 
## loga.prot <- as.numeric(diagram(A.in,do.plot=FALSE)$logact)
## 
## loga.residue <- log10(pl*10^loga.prot)
## 
## species(1:4,loga.prot)
## 
## A.aff <- affinity()
## 
## # again deal with the ionized proteins
## 
## ionize()
## 
## A <- affinity()
## 
## logQ.aff <- affinity(property="logQ")
## 
## logQ.aff <- as.numeric(ionize(A$values,logQ.aff$values))
## 
## # done with equal-affinity reference state
## 
## A.aff <- as.numeric(A.aff$values)[1:4]/pl
## 
## A.aff <- -log10(pl) - (pl-1)/pl*loga.prot + A.aff
## 
## logQ.aff <- logQ.aff[1:4]/pl
## 
## logQ.aff <- log10(pl) + (pl-1)/pl*loga.prot + logQ.aff
## 
## ionize(FALSE)
## 
## # summarize the results
## 
## out <- data.frame(logK=logK,logQ.act=logQ.act,A.act=A.act,loga.prot=loga.prot,loga.residue=loga.residue,logQ.aff=logQ.aff,A.aff=A.aff)
## 
## rownames(out) <- po
## 


###################################################
### code chunk number 4: Gresidue
###################################################
### plot amino acid Gibbs energies
g <- colnames(thermo$group)
i <- info(g)
G <- info(i)$G
# add protein backbone to sidechain groups to get residues
G <- G[1:20] + G[22]
Gf.aa <- G/1000
i <- i[1:20]
f <- info(i)$formula
M <- element(f)$mass
aa <- aminoacids()
mymar <- c(3.5,4,0.5,2)
par(cex=1.5,mar=mymar,las=1,mgp=c(2.5,1,0))
ihp <- which(aa %in% c("A","G","I","L","M","F","P","W","Y","V"))
col <- rep("blue",length(M))
col[ihp] <- "red"
plot(M,Gf.aa,pch=aa,xlab="Mass, g/mol",ylab="",col=col)
thermo.axis(axis.label("DG0f","k"),side=2,line=3.2)


###################################################
### code chunk number 5: Gprotein
###################################################
### plot protein Gibbs energies
ip <- info(proteins)
Gf.prot <- info(ip)$G/1000
# highest values of unfolding Gibbs energies reported by Pfeil, 1998
# (Protein Stability and Folding, Springer) convert kJ/mol to kcal/mol
Gunfold <- convert(c(73.6,64.4,79.1,62.8),"cal")
Gfold <- -Gunfold
par(cex=1.5,mar=mymar,las=1,mgp=c(2,0.7,0))
height <- matrix(c(Gf.prot,Gfold),ncol=4,byrow=TRUE)
barplot(height,ylim=c(pretty(Gf.prot)[1],0),ylab=NULL,names=prot,col=c("grey","green4"),border=NA)
thermo.axis(axis.label("DG0f","k"),side=2,line=3.2)


###################################################
### code chunk number 6: 80O2_four
###################################################
# just at logfO2 = -80
par(cex=1.5)
a <- affinity()
d <- diagram(a,do.plot=FALSE)
thermo.plot.new(xlim=c(-85,-65),ylim=c(-8,0),xlab=axis.label("O2"),ylab=as.expression(quote(log~italic(a))),mar=mymar,yline=par("mgp")[1]+1)
myloga <- as.numeric(d$logact)
points(rep(-80,4),myloga,pch=19)
# put some space around the text labels
text(rep(-79.5,3),myloga[c(1,2,4)],prot[c(1,2,4)],adj=0)
text(rep(-80.5,1),myloga[3],prot[3],adj=1)


###################################################
### code chunk number 7: O2_four
###################################################
par(cex=1.5,lwd=1.5)
a <- affinity(O2=c(-85,-65))
d <- diagram(a,ylim=c(-8,0),legend.x=NULL,mar=mymar)
text(-82,-1.3,"MYG")
text(-73.5,-1.2,"CYC")
text(-68,-1.2,"RNAS1")
text(-67.5,-2.5,"LYSC")
points(rep(-80,4),myloga,pch=19,cex=0.7)


###################################################
### code chunk number 8: setupEFTu
###################################################
# set up for EF-Tu study
# read the amino acid compositions of the proteins
f <- system.file("extdata/fasta/EF-Tu.aln",package="CHNOSZ")
# we only use six of the proteins
iip <- c(4:6,1:3)
p <- read.fasta(f,pnff=TRUE)[iip,]
# add the proteins to this session
ip <- add.protein(p)
basis("CHNOS+")
species(paste(p$protein,p$organism,sep="_"))


###################################################
### code chunk number 9: O2_EFTu
###################################################
a <- affinity(O2=c(-80,-66))
col <- c("blue","red","red","red","red","blue")
lty <- c(1,1,1,2,2,2)
# use a small top margin since the titles will be on the slides
mar <- c(3.5,4,0.5,2)
d <- diagram(a,col=col,names=p$organism,legend.x=NULL,mar=mar,cex=1.5,lwd=1.5,lty=lty)
text(-76.5,-3.4,p$organism[1])
text(-78,-3.05,p$organism[2])
lines(c(-77.5,-76.5),c(-3.09,-3.19))
text(-68,-3.08,p$organism[3])
text(-76,-2.63,p$organism[4])
text(-78,-2.89,p$organism[5])
text(-70,-2.8,p$organism[6])


###################################################
### code chunk number 10: O2-T_EFTu
###################################################
a <- affinity(T=c(0,120),O2=c(-80,-50))
diagram(a,names=p$organism,mar=mar,cex=1.5,color=NULL)


###################################################
### code chunk number 11: T-pH_EFTu
###################################################
basis("O2",-60)
a <- affinity(pH=c(0,14),T=c(20,120))
diagram(a,names=p$organism,mar=mymar,cex=1.5,color=NULL)


###################################################
### code chunk number 12: O2-pH_EFTu
###################################################
a <- affinity(pH=c(0,14),O2=c(-80,-60),T=25)
diagram(a,names=p$organism,mar=mymar,cex=1.5,color=NULL)


###################################################
### code chunk number 13: Zprotein_ionized
###################################################
### plot charges of CYC_BOVIN
par(cex=1.5,mar=mymar,las=1,mgp=c(2,0.7,0))
setup_four()
  # add the ionizable groups
  ionize()
  # get the affinities and charges (along equal temperature increments)
  T <- c(25,150,6); pH <- c(0,14,50)
  x <- affinity(pH=pH,T=T)
  z <- ionize(x$values,x$values)
  # plot charges at the temperatures we're interested in
  plot(z[[2]][,1],x=seq(0,14,length.out=50),type="l",xlab="pH",
    ylab="net charge (Z)")  # 25  deg C
  lines(z[[2]][,4],x=seq(0,14,length.out=50),col="red")    # 100 deg C
  lines(z[[2]][,6],x=seq(0,14,length.out=50),col="orange") # 150 deg C
  text(x=c(12,6.8),y=c(-15,-17),labels=paste("T=",c(25,150),sep=""))
  ionize(FALSE)


###################################################
### code chunk number 14: Gprotein_ionized
###################################################
### plot Gibbs energies of ionized CYC_BOVIN
par(cex=1.5,mar=mymar,las=1,mgp=c(2,0.7,0))
  # add the ionizable groups
  ionize()
  # get the affinities and charges (along equal temperature increments)
  T <- c(25,150,6); pH <- c(0,14,50)
  a <- affinity(pH=pH,T=T)
  g <- affinity(pH=pH,T=T,property="G.species")
  myg <- ionize(a$values,g$values)
  # plot energies at the temperatures we're interested in
  plot(myg[[2]][,1]/1000000,x=seq(0,14,length.out=50),type="l",xlab="pH",
    ylab="",ylim=c(-4.8,-3.2))  # 25  deg C
  # add an axis label
  mtext("Gibbs energy (Mcal/mol)",2,line=2.5,cex=par("cex"),las=0)
  lines(myg[[2]][,4]/1000000,x=seq(0,14,length.out=50),col="red")    # 100 deg C
  lines(myg[[2]][,6]/1000000,x=seq(0,14,length.out=50),col="orange") # 150 deg C
  text(x=c(2,2,2),y=c(-4.6,-4.2,-3.9),labels=paste(c(150,100,25)," degC",sep=""))
  ionize(FALSE)


###################################################
### code chunk number 15: Eh-pH_four
###################################################
par(cex=1.5,lwd=1.5)
setup_four()
# we change the basis species to look at Eh instead of logfO2
basis(c("CO2","H2O","NH3","H2S","e-","H+"))
a <- affinity(pH=c(0,14),Eh=c(-0.6,0.6))
diagram(a,names=prot,mar=mymar,color=NULL)
water.lines()


###################################################
### code chunk number 16: Eh-pH_O2
###################################################
par(cex=1.5,lwd=1.5)
# for some reason when compiling the vignette
# as this point the only species loaded is H2O
# but we do want the proteins!
  setup_four()
  # we change the basis species to look at Eh instead of logfO2
  basis(c("CO2","H2O","NH3","H2S","e-","H+"))
  a <- affinity(pH=c(0,14),Eh=c(-0.6,0.6))
diagram(a,names=prot,col="darkgrey",col.names="darkgrey",mar=mymar,color=NULL)
# now overlay the logfO2 contours
basis(c("O2","H+","e-"))
species("H2O")
a <- affinity(pH=c(0,14),Eh=c(-1,1))
diagram(a,what="O2",add=TRUE,cex=1.5,lty=3)
water.lines()


