Last data update: 2014.03.03

R: For two-way designs, plots mcmc.qpcr model predictions gene...
trellisByGeneR Documentation

For two-way designs, plots mcmc.qpcr model predictions gene by gene

Description

For each gene, plots model-predicted values and 95% credible intervals.

Usage

trellisByGene(modelSummary,xFactor,groupFactor,
nrow=2,lineWidth=0.4,whiskerWidth=0.2,pointSize=2.5,
facetScales="free_y",ylab="log(abundance)",
legendPos="bottom",posDodge=0.3)

Arguments

modelSummary

two-way design model summary produced by HPDsummary()

xFactor

factor to form the x-axis

groupFactor

factor to form separate lines on the plot

nrow

number of rows in the resulting trellis plot

lineWidth

line width, passed as 'lwd' to geom_errorbar function (ggplot2)

whiskerWidth

width of the line denoting 95% CI margin, passed as 'width' to geom_errorbar function (ggplot2)

pointSize

passed as 'size' to geom_point function of ggplot2

facetScales

passed as 'scales' to facet_wrap function of ggplot2

ylab

y-axis label

legendPos

passed as 'legend.position' to theme function of ggplot2

posDodge

position dodge, increase for more jitter

Value

A ggplot2 type object

Author(s)

Mikhal V. Matz, UT Austin, matz@utexas.edu

References

Matz MV, Wright RM, Scott JG (2013) No Control Genes Required: Bayesian Analysis of qRT-PCR Data. PLoS ONE 8(8): e71448. doi:10.1371/journal.pone.0071448

Examples


# loading Cq data and amplification efficiencies
data(coral.stress) 
data(amp.eff) 

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=coral.stress,genecols=genecolumns,
condcols=conditions,effic=amp.eff,Cq1=37) 

# fitting the 2-way model
mm=mcmc.qpcr(
	fixed="condition+timepoint+condition:timepoint",
	data=dd,
	nitt=4000 # remark this line to analyze real data!
)

# summarizing results
ss=HPDsummary(mm,data=dd,summ.plot=FALSE)

# plotting predicted means and 95% CIs gene by gene
trellisByGene(ss,xFactor="condition",groupFactor="timepoint")

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(MCMC.qpcr)
Loading required package: MCMCglmm
Loading required package: Matrix
Loading required package: coda
Loading required package: ape
Loading required package: ggplot2
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/MCMC.qpcr/trellisByGene.Rd_%03d_medium.png", width=480, height=480)
> ### Name: trellisByGene
> ### Title: For two-way designs, plots mcmc.qpcr model predictions gene by
> ###   gene
> ### Aliases: trellisByGene
> 
> ### ** Examples
> 
> 
> # loading Cq data and amplification efficiencies
> data(coral.stress) 
> data(amp.eff) 
> 
> 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=coral.stress,genecols=genecolumns,
+ condcols=conditions,effic=amp.eff,Cq1=37) 
> 
> # fitting the 2-way model
> mm=mcmc.qpcr(
+ 	fixed="condition+timepoint+condition:timepoint",
+ 	data=dd,
+ 	nitt=4000 # remark this line to analyze real data!
+ )
$PRIOR
$PRIOR$R
$PRIOR$R$V
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1

$PRIOR$R$nu
[1] 3.002


$PRIOR$G
$PRIOR$G$G1
$PRIOR$G$G1$V
[1] 1

$PRIOR$G$G1$nu
[1] 0




$FIXED
[1] "count~0+gene++gene:condition+gene:timepoint+gene:condition:timepoint"

$RANDOM
[1] "~sample"


                       MCMC iteration = 0

 Acceptance ratio for liability set 1 = 0.000188

 Acceptance ratio for liability set 2 = 0.000587

 Acceptance ratio for liability set 3 = 0.000281

 Acceptance ratio for liability set 4 = 0.000317

                       MCMC iteration = 1000

 Acceptance ratio for liability set 1 = 0.123172

 Acceptance ratio for liability set 2 = 0.320222

 Acceptance ratio for liability set 3 = 0.158875

 Acceptance ratio for liability set 4 = 0.307175

                       MCMC iteration = 2000

 Acceptance ratio for liability set 1 = 0.189891

 Acceptance ratio for liability set 2 = 0.334127

 Acceptance ratio for liability set 3 = 0.223250

 Acceptance ratio for liability set 4 = 0.333968

                       MCMC iteration = 3000

 Acceptance ratio for liability set 1 = 0.220109

 Acceptance ratio for liability set 2 = 0.331111

 Acceptance ratio for liability set 3 = 0.251422

 Acceptance ratio for liability set 4 = 0.340333

                       MCMC iteration = 4000

 Acceptance ratio for liability set 1 = 0.239656

 Acceptance ratio for liability set 2 = 0.327810

 Acceptance ratio for liability set 3 = 0.280687

 Acceptance ratio for liability set 4 = 0.335587
> 
> # summarizing results
> ss=HPDsummary(mm,data=dd,summ.plot=FALSE)
> 
> # plotting predicted means and 95% CIs gene by gene
> trellisByGene(ss,xFactor="condition",groupFactor="timepoint")
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>