### 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()