Next: transfer Up: CHNOSZ examples Previous: get.protein

revisit

revist>   ## Don't show:
revist> data(thermo)
thermo: loaded 1997 aqueous, 3089 total species to thermo$obigt
thermo: loaded 5264 proteins to thermo$ECO
thermo: loaded 6717 proteins to thermo$SGD
thermo: loaded 4155 localizations and 3570 abundances to thermo$yeastgfp

revist> ## End Don't show
revist>   ## carboxylases from different organisms
revist>   # ribulose phosphate carboxylase/oxygenase
revist>   rubisco <- c("Q6JAI0_9RHOD","A5CKC7_9CHRO","RBL_SYNJA","A1E8R4_9CHLO",
revist+     "A8C9T6_9MYCO","A3EQE1_9BACT","A6YF84_9PROT", "RBL_BRAJA",
revist+     "A1RZJ5_THEPD","RBL_METJA","A3DND9_STAMF","RBL_PYRHO")

revist>   # acetyl-coenzyme A carboxylase
revist>   accoaco <- c("Q05KD0_HYDTH","Q9F7M8_PRB01","A6CDM2_9PLAN","A0GZU2_9CHLR",
revist+     "ACCA_DEIRA","A1VC70_DESVV","A7WGI1_9AQUI","Q2JSS7_SYNJA",
revist+     "A4AGS7_9ACTN","ACCA_AQUAE","ACCA_CAUCR","A6VIX9_METM7")

revist>   # calculate affinities as a function of T with
revist>   # buffered logfH2 and CO2
revist>   # adjust the glutathione buffer
revist>   # (more oxidizing than default)
revist>   mod.buffer("GSH-GSSG",c("GSH","GSSG"),logact=c(-3,-7))
mod.buffer: changed state and/or logact of GSH GSSG in GSH-GSSG buffer.

revist>   # add a CO2 gas saturation buffer
revist>   mod.buffer("CO2","CO2","gas",-1)
mod.buffer: added CO2.

revist>   basis(c("CO2","H2O","NH4+","SO4-2","H2","H+"),
revist+     c("CO2",0,-4,-4,"GSH-GSSG",-7))
      C H N O S  Z ispecies   logact state
CO2   1 0 0 2 0  0       69      CO2    aq
H2O   0 2 0 1 0  0        1        0   liq
NH4+  0 4 1 0 0  1       18       -4    aq
SO4-2 0 0 0 4 1 -2       24       -4    aq
H2    0 2 0 0 0  0       66 GSH-GSSG    aq
H+    0 1 0 0 0  1        3       -7    aq

revist>   species(c(rubisco,accoaco))
protein: found Q6JAI0_9RHOD (C2029H3169N551O575S16, 405 residues)
protein: found A5CKC7_9CHRO (C2345H3606N646O687S19, 470 residues)
protein: found RBL_SYNJA (C2366H3652N650O683S21, 474 residues)
protein: found A1E8R4_9CHLO (C1325H2091N383O379S14, 269 residues)
protein: found A8C9T6_9MYCO (C2344H3614N640O690S19, 476 residues)
protein: found A3EQE1_9BACT (C1901H3060N530O543S23, 392 residues)
protein: found A6YF84_9PROT (C1061H1670N312O299S11, 213 residues)
protein: found RBL_BRAJA (C2403H3735N663O702S21, 486 residues)
protein: found A1RZJ5_THEPD (C2215H3468N598O630S13, 443 residues)
protein: found RBL_METJA (C2151H3412N574O627S16, 425 residues)
protein: found A3DND9_STAMF (C2264H3529N603O633S13, 443 residues)
protein: found RBL_PYRHO (C2179H3409N587O615S18, 430 residues)
protein: found Q05KD0_HYDTH (C2218H3568N612O641S19, 444 residues)
protein: found Q9F7M8_PRB01 (C2222H3505N607O675S19, 447 residues)
protein: found A6CDM2_9PLAN (C1611H2578N436O477S13, 322 residues)
protein: found A0GZU2_9CHLR (C2831H4515N821O827S17, 590 residues)
protein: found ACCA_DEIRA (C1494H2419N453O455S8, 316 residues)
protein: found A1VC70_DESVV (C3761H5876N1084O1085S29, 751 residues)
protein: found A7WGI1_9AQUI (C2360H3732N622O702S13, 472 residues)
protein: found Q2JSS7_SYNJA (C1608H2590N464O478S9, 327 residues)
protein: found A4AGS7_9ACTN (C3108H4983N889O981S19, 675 residues)
protein: found ACCA_AQUAE (C1635H2616N440O457S14, 323 residues)
protein: found ACCA_CAUCR (C1517H2440N450O464S11, 320 residues)
protein: found A6VIX9_METM7 (C2443H3887N659O729S17, 494 residues)

revist>   a <- affinity(T=c(0,160))
energy.args: pressure is Psat
energy.args: variable 1 is T at 128 increments from 273.15 to 433.15
affinity: loading buffer species
affinity: loading ionizable protein groups
subcrt: 50 species at 128 values of T and P (wet)
buffer: ( GSH-GSSG ) for activity of H2 (active), CO2 (conserved).
buffer: ( CO2 ) for activity of CO2 (active).

revist>   # create speciation diagram with colors
revist>   col <- c(rep("red",12),rep("blue",12))

revist>   d <- diagram(a,residue=TRUE,color=col,ylim=c(-6,-1),legend.x=NULL)
diagram: immobile component is protein backbone group
diagram: conservation coefficients are 405 470 474 269 476 392 213 486 443 425 443 430 444 447 322 590 316 751 472 327 675 323 320 494
diagram: using residue equivalents
diagram: log total activity of PBB (from species) is 1.017326

revist>   legend("topleft",col=c("red","blue"),lty=1,
revist+     legend=c("ribulose phosphate carboxylase",
revist+     "acetyl-coenzyme A carboxylase"))

revist>   title(main=paste("Calculated relative abundances of",
revist+     "carboxylases from different organisms",sep="\n"))

$\begin{figure}\par \includegraphics{pictures/revisit1} \par \par \end{figure}$

revist>   # ... now on to a species richness diagram
revist>   # all the proteins, then rubisco and accoaco
revist>   draw.diversity(d,"richness",logactmin=-3.6)

revist>   draw.diversity(d,"richness",logactmin=-3.6,
revist+     ispecies=1:12,col="red",add=TRUE)

revist>   draw.diversity(d,"richness",logactmin=-3.6,
revist+     ispecies=13:24,col="blue",add=TRUE)

revist>   legend("bottomleft",col=c("red","blue","black"),lty=1,
revist+     legend=c("ribulose phosphate carboxylase",
revist+     "acetyl-coenzyme A carboxylase","all"))

revist>   title(main=paste("Carboxylases with activities",
revist+     "greater than 10^(-3.6)",sep="\n"))

$\begin{figure}\par \includegraphics{pictures/revisit2} \par \par \end{figure}$

revist>   ## continuing from above ... make a rank-abundance
revist>   ## diagram and fit with a lognormal distribution
revist>   if(require(vegan)) {
revist+     basis("H2",-4)
revist+     a <- affinity()
revist+     logact <- diagram(a,residue=TRUE,do.plot=FALSE)$logact
revist+     act <- 10^as.numeric(logact)
revist+     # we use family=Gamma because our species have activities
revist+     # (i.e., proportional abundances) and not integer counts
revist+     mod <- rad.lognormal(act,family=Gamma)
revist+     plot(mod,main=paste("Relative abundances of carboxylases",
revist+       "fit with lognormal distribution",sep="\n"))
revist+     # calculate Shannon diversity index
revist+     # using revisit (CHNOSZ)
revist+     H1 <- revisit(act,"shannon",as.is=TRUE)
revist+     legend("topright",legend=paste("H =",round(H1,2)),pch=1)
revist+     # using diversity (vegan)
revist+     H2 <- diversity(matrix(act,nrow=1))
revist+     stopifnot(isTRUE(all.equal(H1,H2)))
revist+   }
Loading required package: vegan
This is vegan 1.15-4
affinity: temperature is 25 C
energy.args: pressure is Psat
affinity: loading buffer species
affinity: loading ionizable protein groups
subcrt: 48 species at 298.15 K and 1 bar (wet)
buffer: ( CO2 ) for activity of CO2 (active).
diagram: immobile component is protein backbone group
diagram: conservation coefficients are 405 470 474 269 476 392 213 486 443 425 443 430 444 447 322 590 316 751 472 327 675 323 320 494
diagram: using residue equivalents
diagram: log total activity of PBB (from species) is 1.017326

$\begin{figure}\par \includegraphics{pictures/revisit3} \par \par \end{figure}$

revist>   ### using grep.file, read.fasta, add.protein
revist>   # calculations for Pelagibacter ubique
revist>   f <- system.file("HTCC1062.faa",package="CHNOSZ")

revist>   # line numbers of all entries in the file
revist>   j <- grep.file(f)  # length = 1354

revist>   # locate entries whose names contain DNA
revist>   j <- grep.file(f,"hypothetical")

revist>   # get the amino acid compositions of these protein
revist>   p <- read.fasta(f,j)

revist>   # add these proteins to CHNOSZ's inventory
revist>   i <- add.protein(p)
add.protein: added 199 proteins

revist>   # set up a thermodynamic system
revist>   basis("CHNOS+")
    C H N O S Z ispecies logact state
CO2 1 0 0 2 0 0       69     -3    aq
H2O 0 2 0 1 0 0        1      0   liq
NH3 0 3 1 0 0 0       68     -4    aq
H2S 0 2 0 0 1 0       70     -7    aq
O2  0 0 0 2 0 0     2852    -80   gas
H+  0 1 0 0 0 1        3     -7    aq

revist>   # calculate affinities in logfO2-pH space
revist>   a <- affinity(H2O=c(-10,-2),O2=c(-82,-76),iprotein=i)
affinity: temperature is 25 C
energy.args: pressure is Psat
energy.args: variable 1 is H2O at 128 increments from -10 to -2
energy.args: variable 2 is O2 at 128 increments from -82 to -76
protein: found H2O_RESIDUE (H2O, 0 residues)
protein: found Ala_RESIDUE (C3H5NO, 1 residues)
protein: found Cys_RESIDUE (C3H5NOS, 1 residues)
protein: found Asp_RESIDUE (C4H5NO3, 1 residues)
protein: found Glu_RESIDUE (C5H7NO3, 1 residues)
protein: found Phe_RESIDUE (C9H9NO, 1 residues)
protein: found Gly_RESIDUE (C2H3NO, 1 residues)
protein: found His_RESIDUE (C6H7N3O, 1 residues)
protein: found Ile_RESIDUE (C6H11NO, 1 residues)
protein: found Lys_RESIDUE (C6H12N2O, 1 residues)
protein: found Leu_RESIDUE (C6H11NO, 1 residues)
protein: found Met_RESIDUE (C5H9NOS, 1 residues)
protein: found Asn_RESIDUE (C4H6N2O2, 1 residues)
protein: found Pro_RESIDUE (C5H7NO, 1 residues)
protein: found Gln_RESIDUE (C5H8N2O2, 1 residues)
protein: found Arg_RESIDUE (C6H12N4O, 1 residues)
protein: found Ser_RESIDUE (C3H5NO2, 1 residues)
protein: found Thr_RESIDUE (C4H7NO2, 1 residues)
protein: found Val_RESIDUE (C5H9NO, 1 residues)
protein: found Trp_RESIDUE (C11H10N2O, 1 residues)
protein: found Tyr_RESIDUE (C9H9NO2, 1 residues)
affinity: loading ionizable protein groups
subcrt: 44 species at 298.15 K and 1 bar (wet)
50..100..150..

revist>   # calculate the logarithms of activities
revist>   d <- diagram(a,do.plot=FALSE,mam=FALSE)
diagram: immobile component is protein backbone group
diagram: using residue equivalents
diagram: log total activity of PBB (from species) is 1.648994

revist>   # show the protein richness
revist>   draw.diversity(d,"richness",logactmin="")

revist>   mtitle(c("Richness of hypothetical proteins in",
revist+     expression(italic("Pelagibacter ubique"))))

$\begin{figure}\par \includegraphics{pictures/revisit4} \par \par \end{figure}$

revist>   # show the coefficient of variation of activities
revist>   draw.diversity(d,"cv")

revist>   mtitle(c("Coefficient of variation of hypothetical",
revist+     expression("protein activities in"~italic("P. ubique"))))

$\begin{figure}\par \includegraphics{pictures/revisit5} \par \par \end{figure}$

Next: transfer Up: CHNOSZ examples Previous: get.protein