2014年9月24日 星期三

R script (grofit)

#before everything starts, type this script/

library(grofit)


#open 2 R console in terminal
open -n -a R.app

#print working directory
getwd()



Workflow

#Define option (e.g. interactive = FALSE means manually choose each growth curve

#
grofit.control(neg.nan.act = FALSE, clean.bootstrap = TRUE,
               suppress.messages = FALSE, fit.opt = "b",
               log.x.gc = FALSE, log.y.gc = FALSE,
               interactive = TRUE, nboot.gc = 0,
               smooth.gc= NULL, model.type=c("logistic",
               "richards","gompertz", "gompertz.exp"),
               have.atleast = 6, parameter = 9, smooth.dr = NULL,
               log.x.dr = FALSE, log.y.dr = FALSE, nboot.dr = 0)

MyOpt1 <- grofit.control(fit.opt="b",interactive = FALSE,log.y.gc=TRUE,smooth.gc=0.5)

#use specific model

MyOpt2 <- grofit.control(fit.opt="m",model.type=c("gompertz"),interactive = FALSE,log.y.gc=TRUE,smooth.gc=0.5)

#read data from csv file

data<-read.table ("/Users/andesli/Desktop/20140926_Mutantbglucan_merged_log.csv",header=FALSE,sep=",")

#Set time

timepoints<-1:44
OR
#Set time from 0.5 to 22, interval = 0.5
timepoint2<-seq(0.5,22, by=0.5)

#Input time into a matrix form

time<-matrix(rep(timepoints,86),nrow=86,ncol=44,byrow=TRUE)

#Fit data with grofit model and FALSE mean turn off dose response curve function/
TestRun1 <- grofit(time, data, FALSE, MyOpt1)
#plot(TestRun1) is not working
OR
TestRun1<-gcFit(time,data,control=MyOpt1)

Plot Graph

#set a new pop up window

#make an empty chart
plot(1,1,xlim=c(0,42),ylim=c(0,1),type="n",ann=FALSE)
title(main="Mutant_1",xlab="Time (hours)",ylab="log(1+OD at 600nm)")

dev.new(width=8,height=3)
par(mfrow=c(1,3),mar=c(5.1,4.1,3.1,1.1))

#define dot color

colData<-c("black","cyan","magenta","green","blue","orange","grey","red","darkgreen","pink")

#pch = dot shape, slope = TURE-> tangent, add = whether plot in the existing window
#opt=m => plot model, s => spline

plot(TestRun1, add=FALSE, opt="m", pch=1:10, cex=1,colModel=1,colData=colData, slope=TRUE)

#plot specific curve in large data

plot(TestRun1$gcFittedModel[[no. of row], add=FALSE, opt="m", pch=1, cex=1,colModel=1,colData=1, slope=TRUE)
OR
#for data cannot fit model
plot(TestRun$gcFittedSpline[[no. of row]], add=FALSE, opt="s", pch=1, cex=1,colSpline=1,colData=1, slope=FALSE,raw=TRUE)

# call out legend

legend("topleft",c(TestRun_log$gcFittedModel[[no. of row]]$gcID))
legend("topright",c("RCM","beta glucan"),pch=c(1,2),col=c(1,2))

#plot single curve with best fitted curve and tagent

plot(TestRun_log$gcFittedModel[[no. of row]]$raw.time,TestRun_log$gcFittedModel[[no. of row]]$raw.data,xlab="Time (hours)",ylab="log(1+OD (nm600))")
plot(plot(TestRun_log$gcFittedModel[[no. of row]],add=TURE)

#plot multiple curve at the same time and save

for (i in 1:10)
{
pdf(paste("/Users/andesli/Desktop/Mutant_",b,".pdf")
+ plot(1,1,xlim=c(0,42),ylim=c(0,1),type="n",ann=FALSE)
+ title(main=paste("Mutant ",i),xlab="Time (hours)",ylab="log(1+OD at 600nm)")
+ legend("topright",c("RCM","beta glucan"),pch=c(1,2),col=c(1,2))
+ plot(bifido_merged$gcFittedModel[[i+111]],add=TRUE,pch=1,cex=1,colModel=1,colData=1,slope=TRUE,raw=TRUE)
+ plot(bifido_merged$gcFittedModel[[i]],add=TRUE,pch=2,cex=1,colModel=2,colData=2,slope=TRUE,raw=TRUE)
dev.off()
}



#Export data to a .txt

sink(file = "/Users/andesli/Desktop/TestRun1.txt")
TestRun1
sink()

#Run R in terminal with script, output in same directory

R CMD batch "filename"

Reference link
#Generic plot function for gcFit objects
http://rgm3.lab.nig.ac.jp/RGM/R_rdfile?f=grofit/man/plot.gcFit.Rd&d=R_CC

#pch shape chart
http://www.harding.edu/fmccown/r/

#Scatterplots tutorial
http://www.statmethods.net/graphs/scatterplot.html


沒有留言:

張貼留言