Imperial College/Courses/Fall2009/Synthetic Biology (MRes class)/'R' Tutorial/Practical: Difference between revisions
| Line 169: | Line 169: | ||
* Modify the code so that the initials values used by the non-linear regression can be read from the parameter file you created in Part 2a. | * Modify the code so that the initials values used by the non-linear regression can be read from the parameter file you created in Part 2a. | ||
* Plot on the same graph both fits (linear regression and non-linear regression) | * Plot on the same graph both fits (linear regression and non-linear regression) | ||
==Part 3: Gene Expression Modelling== | |||
In this section, we will use 'R' to simulate a simple constitutive gene expression model. | |||
The model is given by a system of 2 differential equations to account for the expression of mRNA molecules and Protein molecules: | |||
... | |||
model | |||
... | |||
The script below uses this model to simulate mRNA and Protein time series. | |||
<syntax type='R'> | |||
# | |||
# Gene Expression ODE model | |||
# | |||
params <- c(k1=0.1, d1=0.05, k2=0.1, d2=0.0025) | |||
times <- seq(from=1, to=3600, by=10) | |||
my.atol <- c(1e-6, 1e-6) | |||
geneExpressionModel <-function(t, y, p) | |||
{ | |||
dm <- p["k1"] - p["d1"]*y["m"] | |||
dp <- y["m"]*p["k2"] - p["d2"]*y["p"] | |||
res<-c(dm, dp) | |||
list(res) | |||
} | |||
require(odesolve) | |||
results <- lsoda(c(m=0,p=0),times,geneExpressionModel, params, rtol=1e-4, atol= my.atol) | |||
print(results) | |||
plot(results[,1], results[,2], xlim=c(0,3800), ylim=c(0,100)) | |||
lines(results[,1], results[,3]) | |||
</syntax> | |||
===Part 3: Questions=== | |||
Revision as of 05:53, 9 October 2009
<html> <body> <!-- Start of StatCounter Code --> <script type="text/javascript"> var sc_project=3315864; var sc_invisible=0; var sc_partition=36; var sc_security="8bb2efcd"; </script>
<script type="text/javascript" src="http://www.statcounter.com/counter/counter_xhtml.js"></script><noscript><div class="statcounter"><a class="statcounter" href="http://www.statcounter.com/"><img class="statcounter" src="http://c37.statcounter.com/3315864/0/8bb2efcd/0/" alt="blog stats" /></a></div></noscript> <!-- End of StatCounter Code -->
</body> </html>
Introduction to 'R'
Part 1: Random Walks
Below is a basic 'R' script to simulate a random walk along the X axis. Random walks are important to model as they relate for example to the phenomenon of diffusion.
- Copy and save the script into a file.
- Run this script within 'R'
<syntax type='R'>
- Random Walk 1D / 1 Trace
-
- number of steps
nbSteps <- 1000
- sampling and path construction
- First method: Looping
- draw from uniform distribtion between [0,1]
x <- runif(n = nbSteps)
- initialise path variable
xSteps <- rep(0, nbSteps) #as opposed to x1Steps <-0
- build independent steps series
for(i in 1:nbSteps){ if(x[i]>0.5){ xSteps[i] <- 1} else{ xSteps[i] <- -1 } }
- build path from independent steps
xPath <- cumsum(xSteps)
- Second method: vectorisation (faster)
- x <- runif(n = nbSteps)
- xSteps <- ifelse(x>0.5, 1, -1)
- xPath <- cumsum(xSteps)
- plot path
plot(xPath, 1:nbSteps, type='l', col='red', lwd =2,
main = 'Random Walk 1D',
xlab = "X coordinate",
ylab = "Steps")
</syntax>
Questions
- Modify this script so that you can overlay as many independent paths as you want (Tip = use the function capability from 'R').
- Write a function to only return the last X coordinate after nbSteps.
- Use the previous function to build the distribution of X coordinates from 1000 independent simulations, after nbSteps.
- Extend the first script to be able to simulate 2D random walks.
Part 2: Enzymatic reaction analysis
In this section, we will be looking at enzymatic reaction data, using a Michaelis-Menten model. The purpose of the 'R' script is to automatise enzyme characterisation from experimental data.
Part 2a: Linear regression analysis
- Copy this script and run it.
<syntax type='R'>
- Enzymatic analysis: Linear regression analysis
concentration <- c(0.3330, 0.1670, 0.0833, 0.0416, 0.0208, 0.0104, 0.0052) velocity <- c(3.636, 3.636, 3.236, 2.666, 2.114, 1.466, 0.866)
df <- data.frame(concentration, velocity) plot(df$concentration, df$velocity)
df$ytrans <- 1/df$velocity df$xtrans <- 1/df$concentration
par(ask=T) plot(df$xtrans, df$ytrans)
lmfit <- lm(ytrans~xtrans, data=df)
summary(lmfit)
coefficients <- coef(lmfit) print(coefficients)
par(ask=T) plot(df$xtrans, df$ytrans) abline(coefficients)
Vmax <- 1/coef(lmfit)[1] Km <- Vmax*coef(lmfit)[2]
print(Vmax) print(Km)
</syntax>
Part 2a:Questions
- Using the help functions in 'R', add comments to the script to detail each command. Use also 'plot' arguments to make plot labels more explicit.
- Modify the script so that experimental data can be read from a file, and analytical results can be exported into a file.
Part 2b: Non-linear regression
In this section, we will use a non-linear regression method to estimate the Michaelis-Menten parameters from the data. The non-linear regression method explored here is the least-square method.
<syntax type'R'>
- Enzymatic analysis: Non-linear regression (Least-square optimisation)
concentration <- c(0.3330, 0.1670, 0.0833, 0.0416, 0.0208, 0.0104, 0.0052) velocity <- c(3.636, 3.636, 3.236, 2.666, 2.114, 1.466, 0.866)
df <- data.frame(concentration, velocity) plot(df$concentration, df$velocity)
nlsfit <- nls(velocity~Vmax*concentration/(Km+concentration),data=df, start=list(Km=0.0166, Vmax=3.852))
summary(nlsfit)
par(ask=T) plot(df$concentration, df$velocity, xlim=c(0,0.4), ylim=c(0,4))
x <- seq(0, 0.4, length=100) y2 <- (coef(nlsfit)["Vmax"]*x)/(coef(nlsfit)["Km"]+x) lines(x, y2)
</syntax>
Part 2:Questions
- As before, using the 'R' help comment this code.
- Modify the code so that the initials values used by the non-linear regression can be read from the parameter file you created in Part 2a.
- Plot on the same graph both fits (linear regression and non-linear regression)
Part 3: Gene Expression Modelling
In this section, we will use 'R' to simulate a simple constitutive gene expression model. The model is given by a system of 2 differential equations to account for the expression of mRNA molecules and Protein molecules:
... model ...
The script below uses this model to simulate mRNA and Protein time series.
<syntax type='R'>
- Gene Expression ODE model
params <- c(k1=0.1, d1=0.05, k2=0.1, d2=0.0025)
times <- seq(from=1, to=3600, by=10)
my.atol <- c(1e-6, 1e-6)
geneExpressionModel <-function(t, y, p) { dm <- p["k1"] - p["d1"]*y["m"] dp <- y["m"]*p["k2"] - p["d2"]*y["p"] res<-c(dm, dp)
list(res)
}
require(odesolve)
results <- lsoda(c(m=0,p=0),times,geneExpressionModel, params, rtol=1e-4, atol= my.atol)
print(results)
plot(results[,1], results[,2], xlim=c(0,3800), ylim=c(0,100)) lines(results[,1], results[,3])
</syntax>
