HTqPCR.R
In HTqPCR: Automated analysis of high-throughput qPCR data

### R code from vignette source 'HTqPCR.Rnw'

###################################################
### code chunk number 1: Prepare parameters
###################################################
options(width=65)
set.seed(123)


###################################################
### code chunk number 2: Load package
###################################################
library("HTqPCR")


###################################################
### code chunk number 3: Extract R code (eval = FALSE)
###################################################
## all.R.commands <- system.file("doc", "HTqPCR.Rnw", package = "HTqPCR")
## Stangle(all.R.commands)


###################################################
### code chunk number 4: All functions
###################################################
ls("package:HTqPCR")


###################################################
### code chunk number 5: Load example data
###################################################
data(qPCRraw)
data(qPCRpros)
class(qPCRraw)


###################################################
### code chunk number 6: Information contained in qPCRsets
###################################################
slotNames(qPCRraw)
phenoData(qPCRraw)
pData(qPCRraw)
pData(qPCRraw) <- data.frame(Genotype=rep(c("A", "B"), each=3), Replicate=rep(1:3, 2))
pData(qPCRraw)
featureData(qPCRraw)
head(fData(qPCRraw))


###################################################
### code chunk number 7: Example input files
###################################################
path <- system.file("exData", package="HTqPCR")
head(read.delim(file.path(path, "files.txt")))


###################################################
### code chunk number 8: Read raw data
###################################################
files <- read.delim(file.path(path, "files.txt"))
raw <- readCtData(files=files$File, path=path)


###################################################
### code chunk number 9: Show qPCRset data object
###################################################
show(raw)


###################################################
### code chunk number 10: Ct overview ex 1
###################################################
g <- featureNames(raw)[1:10]
plotCtOverview(raw, genes=g, xlim=c(0,50), groups=files$Treatment, conf.int=TRUE,
ylim=c(0,55))


###################################################
### code chunk number 11: Ct overview ex 2
###################################################
plotCtOverview(raw, genes=g, xlim=c(0,50), groups=files$Treatment,
calibrator="Control")


###################################################
### code chunk number 12: HTqPCR.Rnw:242-245
###################################################
par(mfrow=c(2,1))
g <- featureNames(raw)[1:10]
plotCtOverview(raw, genes=g, xlim=c(0,50), groups=files$Treatment, conf.int=TRUE,
ylim=c(0,55))
plotCtOverview(raw, genes=g, xlim=c(0,50), groups=files$Treatment,
calibrator="Control")


###################################################
### code chunk number 13: Ct card ex 1
###################################################
plotCtCard(raw, col.range=c(10,35), well.size=2.6)


###################################################
### code chunk number 14: Ct card ex 2
###################################################
featureClass(raw) <- factor(c("Marker", "TF", "Kinase")[sample(c(1,1,2,2,1,3),
384, replace=TRUE)])
plotCtCard(raw, plot="class", well.size=2.6)


###################################################
### code chunk number 15: HTqPCR.Rnw:267-268
###################################################
plotCtCard(raw, col.range=c(10,35), well.size=2.6)


###################################################
### code chunk number 16: HTqPCR.Rnw:270-271
###################################################
featureClass(raw) <- factor(c("Marker", "TF", "Kinase")[sample(c(1,1,2,2,1,3),
384, replace=TRUE)])
plotCtCard(raw, plot="class", well.size=2.6)


###################################################
### code chunk number 17: Ct replicates
###################################################
plotCtReps(qPCRraw, card=2, percent=20)


###################################################
### code chunk number 18: HTqPCR.Rnw:288-289
###################################################
plotCtReps(qPCRraw, card=2, percent=20)


###################################################
### code chunk number 19: Ct variation ex 1
###################################################
raw.mix <- raw
#exprs(raw.mix)[,6] <- sample(exprs(raw[,6]))
plotCtVariation(raw.mix, variation="sd", log=TRUE, main="SD of replicated features", col="lightgrey")


###################################################
### code chunk number 20: Ct variation ex 2
###################################################
raw.variation <- plotCtVariation(raw.mix, type="detail", add.featurenames=TRUE, pch=" ", cex=1.2)


###################################################
### code chunk number 21: Ct variation ex 2 in detail
###################################################
names(raw.variation)
head(raw.variation[["Var"]][,1:4])
head(raw.variation[["Mean"]][,1:4])
apply(raw.variation[["Var"]][,3:7], 2, summary)
colSums(raw.variation[["Var"]][,3:7]>20)


###################################################
### code chunk number 22: HTqPCR.Rnw:328-329
###################################################
raw.mix <- raw
#exprs(raw.mix)[,6] <- sample(exprs(raw[,6]))
plotCtVariation(raw.mix, variation="sd", log=TRUE, main="SD of replicated features", col="lightgrey")


###################################################
### code chunk number 23: HTqPCR.Rnw:331-332
###################################################
raw.variation <- plotCtVariation(raw.mix, type="detail", add.featurenames=TRUE, pch=" ", cex=1.2)


###################################################
### code chunk number 24: Plot Ct categories ex 1
###################################################
raw.cat <- raw
plotCtCategory(raw.cat)


###################################################
### code chunk number 25: Plot Ct categories ex 2
###################################################
plotCtCategory(raw.cat, stratify="class")


###################################################
### code chunk number 26: HTqPCR.Rnw:373-376
###################################################
par(mfrow=c(2,1))
raw.cat <- raw
plotCtCategory(raw.cat)
plotCtCategory(raw.cat, stratify="class")


###################################################
### code chunk number 27: Plot Ct categories ex 3
###################################################
plotCtCategory(raw.cat, by.feature=TRUE, cexRow=0.1)


###################################################
### code chunk number 28: HTqPCR.Rnw:391-392
###################################################
plotCtCategory(raw.cat, by.feature=TRUE, cexRow=0.1)


###################################################
### code chunk number 29: Normalise data
###################################################
q.norm <- normalizeCtData(raw.cat, norm="quantile")
sr.norm <- normalizeCtData(raw.cat, norm="scale.rank")
nr.norm <- normalizeCtData(raw.cat, norm="norm.rank")
d.norm <- normalizeCtData(raw.cat, norm="deltaCt", deltaCt.genes=c("Gene1", "Gene60"))
g.norm <- normalizeCtData(raw.cat, norm="geometric.mean")


###################################################
### code chunk number 30: Normalisation comparison
###################################################
plot(exprs(raw), exprs(q.norm), pch=20, main="Quantile normalisation", col=rep(brewer.pal(6, "Spectral"), each=384))


###################################################
### code chunk number 31: HTqPCR.Rnw:444-460
###################################################
col <- rep(brewer.pal(6, "Spectral"), each=384)
col2 <- brewer.pal(5, "Dark2")
par(mfrow=c(3,2), mar=c(2,2,2,2))
# All methods individually
plot(exprs(raw), exprs(q.norm), pch=20, main="Quantile normalisation", col=col)
plot(exprs(raw), exprs(sr.norm), pch=20, main="Rank invariant scaling", col=col)
plot(exprs(raw), exprs(nr.norm), pch=20, main="Rank invariant normalisation", col=col)
plot(exprs(raw), exprs(d.norm), pch=20, main="deltaCt normalisation", col=col)
plot(exprs(raw), exprs(g.norm), pch=20, main="Geometric mean normalisation", col=col)
# Just a single sample, across methods
plot(exprs(raw)[,3], exprs(q.norm)[,3], pch=20, col=col2[1], main="Comparison of methods for sample 3", ylim=c(-10,40))
points(exprs(raw)[,3], exprs(sr.norm)[,3], pch=20, col=col2[2])
points(exprs(raw)[,3], exprs(nr.norm)[,3], pch=20, col=col2[3])
points(exprs(raw)[,3], exprs(d.norm)[,3], pch=20, col=col2[4])
points(exprs(raw)[,3], exprs(g.norm)[,3], pch=20, col=col2[5])
legend(8, 40, legend=c("Quantile", "Rank.invariant scaling", "Rank.invariant normalization", "deltaCt", "Geometric.mean"), col=col2, lwd=2, bty="n")


###################################################
### code chunk number 32: Ct correlations
###################################################
plotCtCor(raw, main="Ct correlation")


###################################################
### code chunk number 33: HTqPCR.Rnw:540-541
###################################################
plotCtCor(raw, main="Ct correlation")


###################################################
### code chunk number 34: Summary of Ct values
###################################################
summary(raw)


###################################################
### code chunk number 35: Ct density
###################################################
plotCtDensity(sr.norm)


###################################################
### code chunk number 36: Ct histogram
###################################################
plotCtHistogram(sr.norm)


###################################################
### code chunk number 37: HTqPCR.Rnw:571-574
###################################################
par(mfrow=c(1,2), mar=c(3,3,2,1))
plotCtDensity(sr.norm)
plotCtHistogram(sr.norm)


###################################################
### code chunk number 38: HTqPCR.Rnw:585-592
###################################################
par(mfrow=c(3,2), mar=c(2,2,2,1))
plotCtDensity(qPCRraw, main="Raw Ct values")
plotCtDensity(q.norm, main="quantile")
plotCtDensity(sr.norm, main="scale.rankinvariant")
plotCtDensity(nr.norm, main="norm.rankinvariant")
plotCtDensity(d.norm, main="deltaCt")
plotCtDensity(g.norm, main="geometric.mean")


###################################################
### code chunk number 39: Ct boxes
###################################################
plotCtBoxes(sr.norm, stratify="class")


###################################################
### code chunk number 40: HTqPCR.Rnw:607-608
###################################################
plotCtBoxes(sr.norm, stratify="class")


###################################################
### code chunk number 41: Ct scatter ex 1
###################################################
plotCtScatter(sr.norm, cards=c(1,2), col="type", diag=TRUE)


###################################################
### code chunk number 42: Ct scatter ex 2
###################################################
plotCtScatter(sr.norm, cards=c(1,4), col="class", diag=TRUE)


###################################################
### code chunk number 43: HTqPCR.Rnw:628-631
###################################################
par(mfrow=c(1,2), mar=c(3,3,2,1))
plotCtScatter(sr.norm, cards=c(1,2), col="type", diag=TRUE)
plotCtScatter(sr.norm, cards=c(1,4), col="class", diag=TRUE)


###################################################
### code chunk number 44: Ct pairs
###################################################
plotCtPairs(sr.norm, col="type", diag=TRUE)


###################################################
### code chunk number 45: HTqPCR.Rnw:649-650
###################################################
plotCtPairs(sr.norm, col="type", diag=TRUE)


###################################################
### code chunk number 46: Ct heatmap
###################################################
plotCtHeatmap(raw, gene.names="", dist="euclidean")


###################################################
### code chunk number 47: HTqPCR.Rnw:669-670
###################################################
plotCtHeatmap(raw, gene.names="", dist="euclidean")


###################################################
### code chunk number 48: CV across samples
###################################################
plotCVBoxes(qPCRraw, stratify="class")
plotCVBoxes(qPCRraw, stratify="type")


###################################################
### code chunk number 49: HTqPCR.Rnw:688-691
###################################################
par(mfrow=c(1,2), mar=c(2,2,2,1))
plotCVBoxes(qPCRraw, stratify="class")
plotCVBoxes(qPCRraw, stratify="type")


###################################################
### code chunk number 50: Cluster Ct
###################################################
clusterCt(sr.norm, type="samples")


###################################################
### code chunk number 51: HTqPCR.Rnw:717-718
###################################################
clusterCt(sr.norm, type="samples")


###################################################
### code chunk number 52: Plot subclusters
###################################################
cluster.list <- clusterCt(sr.norm, type="genes", n.cluster=6, cex=0.5)


###################################################
### code chunk number 53: HTqPCR.Rnw:737-738
###################################################
cluster.list <- clusterCt(sr.norm, type="genes", n.cluster=6, cex=0.5)


###################################################
### code chunk number 54: Principal components analysis
###################################################
plotCtPCA(qPCRraw)
plotCtPCA(qPCRraw, features=FALSE)


###################################################
### code chunk number 55: HTqPCR.Rnw:756-759
###################################################
par(mfrow=c(1,2), mar=c(2,2,2,1))
plotCtPCA(qPCRraw)
plotCtPCA(qPCRraw, features=FALSE)


###################################################
### code chunk number 56: Object history
###################################################
getCtHistory(sr.norm)


###################################################
### code chunk number 57: Multiple samples per card
###################################################
# Example with 2 or 4 samples per 384 well card.
sample2.order <- rep(c("subSampleA", "subSampleB"), each=192)
sample4.order <- rep(c("subA", "subB", "subC", "subD"), each=96)
# Splitting the data into all individual samples
qPCRnew2 <- changeCtLayout(sr.norm, sample.order=sample2.order)
show(qPCRnew2)
qPCRnew4 <- changeCtLayout(sr.norm, sample.order=sample4.order)
show(qPCRnew4)


###################################################
### code chunk number 58: Card history
###################################################
getCtHistory(qPCRnew4)


###################################################
### code chunk number 59: Example SDS data
###################################################
path <- system.file("exData", package="HTqPCR")
cat(paste(readLines(file.path(path, "SDS_sample.txt"), n=19), "\n"))


###################################################
### code chunk number 60: Example SDS data 2
###################################################
readLines(file.path(path, "SDS_sample.txt"), n=20)


###################################################
### code chunk number 61: SDS format
###################################################
path <- system.file("exData", package = "HTqPCR")
raw <- readCtData(files="SDS_sample.txt", path=path, format="SDS")
show(raw)


###################################################
### code chunk number 62: SDS format
###################################################
path <- system.file("exData", package = "HTqPCR")
raw <- readCtData(files="LightCycler_sample.txt", path=path, format="LightCycler")
show(raw)


###################################################
### code chunk number 63: CFX format
###################################################
path <- system.file("exData", package = "HTqPCR")
raw <- readCtData(files="CFX_sample.txt", path=path, format="CFX", n.features=330)
show(raw)


###################################################
### code chunk number 64: BioMark format
###################################################
exPath <- system.file("exData", package="HTqPCR")
raw1 <- readCtData(files="BioMark_sample.csv", path=exPath, format="BioMark", n.features=48, n.data=48)
dim(raw1)
raw2 <- readCtData(files="BioMark_sample.csv", path=exPath, format="BioMark", n.features=48*48, n.data=1)
dim(raw2)


###################################################
### code chunk number 65: Ct microfluidic array
###################################################
plotCtArray(raw1)


###################################################
### code chunk number 66: HTqPCR.Rnw:1078-1079
###################################################
plotCtArray(raw1)


###################################################
### code chunk number 67: OpenArray format
###################################################
exPath <- system.file("exData", package="HTqPCR")
raw1 <- readCtData(files="OpenArray_sample.csv", path=exPath, format="OpenArray", n.features=846, n.data=6)
dim(raw1)
raw2 <- readCtData(files="OpenArray_sample.csv", path=exPath, format="OpenArray", n.features=846*6, n.data=1)
dim(raw2)


###################################################
### code chunk number 68: Check HTqPCR news
###################################################
news(Version>1.7, package="HTqPCR")


###################################################
### code chunk number 69: sessionInfo
###################################################
toLatex(sessionInfo())