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


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


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


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


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


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


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


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


###################################################
### chunk number 9: 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))


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


###################################################
### chunk number 11: 
###################################################
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")


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


###################################################
### chunk number 13: 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)


###################################################
### chunk number 14: 
###################################################
plotCtCard(raw, col.range=c(10,35), well.size=2.6)


###################################################
### chunk number 15: 
###################################################
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)


###################################################
### chunk number 16: Set Ct categories
###################################################
raw.cat <- setCategory(raw, groups=files$Treatment, quantile=0.8)


###################################################
### chunk number 17: Plot Ct categories ex 1
###################################################
plotCtCategory(raw.cat)


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


###################################################
### chunk number 19: 
###################################################
par(mfrow=c(2,1))
plotCtCategory(raw.cat)
plotCtCategory(raw.cat, stratify="class")


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


###################################################
### chunk number 21: 
###################################################
plotCtCategory(raw.cat, by.feature=TRUE, cexRow=0.1)


###################################################
### chunk number 22: 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"))


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


###################################################
### chunk number 24: 
###################################################
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)
# Just a single sample, across methods
plot(exprs(raw)[,1], exprs(q.norm)[,1], pch=20, col=col2[1], main="Comparison of methods for sample 1", ylim=c(-10,40))
points(exprs(raw)[,1], exprs(sr.norm)[,1], pch=20, col=col2[2])
points(exprs(raw)[,1], exprs(nr.norm)[,1], pch=20, col=col2[3])
points(exprs(raw)[,1], exprs(d.norm)[,1], pch=20, col=col2[4])
legend(8, 40, legend=c("Quantile", "Rank.invariant scaling", "Rank.invariant normalization", "deltaCt"), col=col2, lwd=2, bty="n")


###################################################
### chunk number 25: Subset data
###################################################
nr.norm[1:10,]
nr.norm[,c(1,3,5)]


###################################################
### chunk number 26: Filter data 1
###################################################
qFilt   <- filterCtData(nr.norm, remove.type="Endogenous Control")
qFilt   <- filterCtData(nr.norm, remove.name=c("Gene1", "Gene20", "Gene30"))
qFilt   <- filterCtData(nr.norm, remove.class="Kinase")
qFilt   <- filterCtData(nr.norm, remove.type=c("Endogenous Control"), remove.name=c("Gene1", "Gene20", "Gene30"))


###################################################
### chunk number 27: Filter data 2
###################################################
qFilt   <- filterCtData(nr.norm, remove.category="Undetermined")
qFilt   <- filterCtData(nr.norm, remove.category="Undetermined", n.category=5)


###################################################
### chunk number 28: IQR plot
###################################################
iqr.values <- apply(exprs(nr.norm), 1, IQR)
hist(iqr.values, n=20, main="", xlab="IQR across samples")
abline(v=1.5, col=2)


###################################################
### chunk number 29: Filter data 3
###################################################
qFilt   <- filterCtData(nr.norm, remove.IQR=1.5)


###################################################
### chunk number 30: 
###################################################
iqr.values <- apply(exprs(nr.norm), 1, IQR)
hist(iqr.values, n=20, main="", xlab="IQR across samples")
abline(v=1.5, col=2)


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


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


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


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


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


###################################################
### chunk number 36: 
###################################################
par(mfrow=c(1,2), mar=c(3,3,2,1))
plotCtDensity(sr.norm)
plotCtHistogram(sr.norm)


###################################################
### chunk number 37: 
###################################################
par(mfrow=c(2,2), mar=c(2,2,2,1))
plotCtDensity(q.norm, main="quantile")
plotCtDensity(sr.norm, main="scale.rankinvariant")
plotCtDensity(nr.norm, main="norm.rankinvariant")
plotCtDensity(d.norm, main="deltaCt")


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


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


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


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


###################################################
### chunk number 42: 
###################################################
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)


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


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


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


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


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


###################################################
### chunk number 48: 
###################################################
plotCtReps(qPCRraw, card=2, percent=20)


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


###################################################
### chunk number 50: 
###################################################
par(mfrow=c(1,2), mar=c(2,2,2,1))
plotCVBoxes(qPCRraw, stratify="class")
plotCVBoxes(qPCRraw, stratify="type")


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


###################################################
### chunk number 52: 
###################################################
clusterCt(sr.norm, type="samples")


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


###################################################
### chunk number 54: Show subcluster
###################################################
c6 <- cluster.list[[6]]
print(c6)
show(sr.norm[c6,])


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


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


###################################################
### chunk number 57: 
###################################################
par(mfrow=c(1,2), mar=c(2,2,2,1))
plotCtPCA(qPCRraw)
plotCtPCA(qPCRraw, features=FALSE)


###################################################
### chunk number 58: Object history
###################################################
getCtHistory(sr.norm)
getCtHistory(qFilt)


###################################################
### chunk number 59: Perform standard t-test
###################################################
qDE.ttest <- ttestCtData(sr.norm[,1:4], groups=files$Treatment[1:4], 
calibrator="Control")
qDE.ttest[1:10,]


###################################################
### chunk number 60: Perform Mann-Whitney test
###################################################
qDE.mwtest <- mannwhitneyCtData(sr.norm[,1:4], groups=files$Treatment[1:4], 
calibrator="Control")
qDE.mwtest[1:10,]


###################################################
### chunk number 61: Perform limma test
###################################################
# Preparing experiment design
design <- model.matrix(~0+files$Treatment)
colnames(design) <- c("Control", "LongStarve", "Starve")
print(design)
contrasts <- makeContrasts(LongStarve-Control, LongStarve-Starve, 
Starve-Control, (Starve+LongStarve)/2-Control, levels=design)
colnames(contrasts) <- c("LS-C", "LS-S", "S-C", "bothS-C")
print(contrasts)
# Reorder data to get the genes in consecutive rows
sr.norm2    <- sr.norm[order(featureNames(sr.norm)),] 
qDE.limma       <- limmaCtData(sr.norm2, design=design, contrasts=contrasts, 
ndups=2, spacing=1)


###################################################
### chunk number 62: limma test output
###################################################
class(qDE.limma)
names(qDE.limma)
qDE.limma[["LS-C"]][1:10,]


###################################################
### chunk number 63: limma summary output
###################################################
qDE.limma[["Summary"]][21:30,]


###################################################
### chunk number 64: Relative quantification ex 1
###################################################
plotCtRQ(qDE.ttest, genes=1:15)


###################################################
### chunk number 65: Relative quantification ex 2
###################################################
plotCtRQ(qDE.limma, p.val=0.085, transform="log10", col="#9E0142")


###################################################
### chunk number 66: 
###################################################
plotCtRQ(qDE.ttest, genes=1:15)


###################################################
### chunk number 67: 
###################################################
plotCtRQ(qDE.limma, p.val=0.085, transform="log10", col="#9E0142")


###################################################
### chunk number 68: Significant Ct
###################################################
plotCtSignificance(qDE.limma, q=sr.norm, 
groups=files$Treatment, target="LongStarve", 
calibrator="Control", genes=featureNames(sr.norm)[11:20], un.col="#3288BD", jitter=0.2)


###################################################
### chunk number 69: 
###################################################
plotCtSignificance(qDE.limma, q=sr.norm, 
groups=files$Treatment, target="LongStarve", 
calibrator="Control", genes=featureNames(sr.norm)[11:20], un.col="#3288BD", jitter=0.2)


###################################################
### chunk number 70: Heatmap significant Ct
###################################################
heatmapSig(qDE.limma, dist="euclidean")


###################################################
### chunk number 71: 
###################################################
heatmapSig(qDE.limma, dist="euclidean")


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


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


###################################################
### chunk number 74: 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)


###################################################
### chunk number 75: Card history
###################################################
getCtHistory(qPCRnew4)


###################################################
### chunk number 76: Combine qPCRset objects
###################################################
q.comb	<- cbind(q.norm[,1:3], sr.norm[,4], nr.norm[,c(1,5,6)])
q.comb
q.comb2	<- rbind(q.norm, sr.norm[1:4,], nr.norm)
q.comb2


###################################################
### chunk number 77: Combined card history
###################################################
getCtHistory(q.comb)


###################################################
### chunk number 78: sessionInfo
###################################################
toLatex(sessionInfo())