# NOT RUN {
# loading Cq data and amplification efficiencies
data(coral.stress)
data(amp.eff)
# extracting a subset of data
cs.short=subset(coral.stress, timepoint=="one")
genecolumns=c(5,6,16,17) # specifying columns corresponding to genes of interest
conditions=c(1:4) # specifying columns containing factors
# calculating molecule counts and reformatting:
dd=cq2counts(data=cs.short,genecols=genecolumns,
condcols=conditions,effic=amp.eff,Cq1=37)
# fitting the model
mm=mcmc.qpcr(
fixed="condition",
data=dd,
controls=c("nd5","rpl11"),
nitt=3000,burnin=2000 # remove this line when analyzing real data!
)
# plotting abundances of individual genes across all conditions
# step 1: defining conditions
cds=list(
control=list(factors=0), # gene-specific intercept
stress=list(factors=c(0,"conditionheat")) # multiple effects will be summed up
)
# step 2: plotting gene after gene on the same panel
HPDplotBygene(model=mm,gene="actin",conditions=cds,col="cyan3",
pch=17,jitter=-0.1,ylim=c(-3.5,15),pval="z")
HPDplotBygene(model=mm,gene="hsp16",conditions=cds,
newplot=FALSE,col="coral",pch=19,jitter=0.1,pval="z")
# step 3: adding legend
legend(0.5,10,"actin",lty=1,col="cyan3",pch=17,bty="n")
legend(0.5,7,"hsp16",lty=1,col="coral",pch=19,bty="n")
# }
Run the code above in your browser using DataLab