User:Vincent Rouilly/Computational Biology With R

Computational Biology with R

Books

• R Programming for Bioinformatics, Robert Gentleman, Fred Hutchinson, 2008
• Introduction to Scientific Programming and Simulation Using R, Owen Jones et al, 2009
• More R Books (from R-Project.org)

Software resources

• StatET with Eclipse, develop in R within a great IDE.
• Sweave (literate programming / reproducible science))

Tutorials

RSBML lib

• Install library
  • > source("http://bioconductor.org/biocLite.R")
  • > biocLite("rsbml")
• Load library
  • > library("rsbml")

• Load SBML model
• Simulate a SBML model
  • > dom <- rsbml_read("/my_sbml_file.xml") ## e.g. dom <- rsbml_read(system.file("sbml", "GlycolysisLayout.xml", package = "rsbml"))
  • > mod <- model(dom)
  • > sub <- new("SOSSubject", mod)
  • > exp <- new("SOSExperiment", subject=sub)
  • > simulate(exp) (Warning: it does seem to work for me at the moment)

Running SBML models in R

1. Install SBMLR package
   1. Start R, and enter:
   2. > source("http://bioconductor.org/biocLite.R")
   3. > biocLite("SBMLR")
2. Load library 'SBMLR' in R by typing:
   1. > library("SBMLR")
3. Download sample SBML file describing a A-->B reaction: media:here
4. Read SBML file into R
   1. > myModel <- readSBML("AtoB.xml")
5. Simulate model between [0,100] with 100 points
   1. > results=simulate(myModel,seq(0,100,1))
6. Plot results
   1. > attach(results)
   2. > par(mfrow=c(2,1))
   3. > plot(time,s1,type="l")
   4. > plot(time,s2,type="l")
   5. > detach(results)

Perturbation analysis on compound concentration

1. Simulate first section between [0,50], where s1=100 at t=0 (following SBML description)
   1. > myModel=readSBML("AtoB.xml")
   2. > result_1 = simulate(myModel,seq(0,50,1))
2. Simulate second section between [50,100], where s1=100 at t=50
   1. > myModel$species$s1$ic=50
   2. > results_2 = simulate(myModel,seq(50,100,1))
   3. > results=data.frame(rbind(results_1, results_2))
   4. > attach(results)
   5. > par(mfrow=c(2,1))
   6. > plot(time,IMP,type="l")
   7. > plot(time,HX,type="l")
   8. > par(mfrow=c(1,1))
   9. > detach(results)