HPDsummary {MCMC.qpcr}

Summarizes and plots results of mcmc.qpcr function series.

Description

Calculates abundances of each gene across factor combinations; calculates pairwise differences between all factor combinations and their significances for each gene; plots results as bar or line graphs with credible intervals (ggplot2) NOTE: only works for experiments involving a single multi-level fixed factor or two fully crossed multi-level fixed factors.

Usage

HPDsummary(model, data, xgroup=NULL,genes = NA, relative = FALSE, 
log.base = 2, summ.plot = TRUE, ptype="z", ...)

Arguments

model

Model generated by mcmc.qpcr(),mcmc.qpcr.lognormal() or mcmc.qpcr.classic()

data

Dataset used to build the model (returned by cq2counts() or cq2log())

xgroup

The factor to form the x-axis of the plot. By default the first factor in the model will be used.

genes

A vector of gene names to summarize and plot. If left unspecified, all genes will be summarized.

relative

Whether to plot absolute transcript abundances (relative = FALSE) or fold- changes relative to the sample that is considered to be "global control" (relative = TRUE). The "global control" is the combination of factors that served as a reference during model fitting, either because it is alphanumerically first (that happens by default) or because it has been explicitly designated as such using relevel() function (see tutorial).

log.base

Base of the logarithm to use.

summ.plot

By default, the function generates a summary plot, which is a line-points-95% credible intervals plot of log(absolute abundances) with 'relative=FALSE' and a more typical bar graph of log(fold change relative to the control), again with 95% credible intervals, with 'relative=TRUE'. Specify 'summ.plot=FALSE' if you don't want the summary plot.

ptype

Which type of p-values to use. By default p-values based on the Bayesian z-score are used. Specify 'ptype="mcmc"' to output more conventional p-values based on MCMC sampling (these will be limited on the lower end by the size of MCMC sample).

...

Additional options for summaryPlot() function. Among those, 'x.order' can be a vector specifying the order of factor levels on the x-axis.

Value

A list of three items:

summary

Summary table containing calculated abundances, their SD and 95% credible limits

geneWise

A series of matrices listing pairwise differences between factor combinations (upper triangle) and corresponding p-values (lower triangle)

ggPlot

the ggplot2 object for plotting. See http://docs.ggplot2.org/0.9.2.1/theme.html for ways to modify it, such as add text, rotate labels, change fonts, etc.

Author(s)

Mikhail V. Matz, University of Texas at Austin <[e-mail]>

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

See Also

See function summaryPlot() for plotting the summary table in other ways.

Examples

data(beckham.data)
data(beckham.eff)

# analysing the first 5 genes 
# (to try it with all 10 genes, change the line below to gcol=4:13)
gcol=4:8 
ccol=1:3 # columns containing experimental conditions

# recalculating into molecule counts, reformatting
qs=cq2counts(data=beckham.data,genecols=gcol,
condcols=ccol,effic=beckham.eff,Cq1=37)

# creating a single factor, 'treatment.time', out of 'tr' and 'time'
qs$treatment.time=as.factor(paste(qs$tr,qs$time,sep="."))

# fitting a naive model
naive=mcmc.qpcr(
	fixed="treatment.time",
	data=qs,
	nitt=3000,burnin=2000 # remove this line in actual analysis!
)

#summary plot of inferred abundances
# s1=HPDsummary(model=naive,data=qs)

#summary plot of fold-changes relative to the global control
s0=HPDsummary(model=naive,data=qs,relative=TRUE)

# pairwise differences and their significances for each gene:
s0$geneWise