NormqPCR.R
In NormqPCR: Functions for normalisation of RT-qPCR data

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

###################################################
### code chunk number 1: read.qPCR.tech.reps
###################################################
library(ReadqPCR) # load the ReadqPCR library
library(NormqPCR)
path <- system.file("exData", package = "NormqPCR")
qPCR.example.techReps <- file.path(path, "qPCR.techReps.txt")
qPCRBatch.qPCR.techReps <- read.qPCR(qPCR.example.techReps)
rownames(exprs(qPCRBatch.qPCR.techReps))[1:8]


###################################################
### code chunk number 2: combine read.qPCR.tech.reps
###################################################
combinedTechReps <- combineTechReps(qPCRBatch.qPCR.techReps)
combinedTechReps


###################################################
### code chunk number 3: taqman read
###################################################
path <- system.file("exData", package = "NormqPCR")
taqman.example <- file.path(path, "/example.txt")
qPCRBatch.taqman <- read.taqman(taqman.example)


###################################################
### code chunk number 4: taqman detector example
###################################################
exprs(qPCRBatch.taqman)["Ccl20.Rn00570287_m1",]


###################################################
### code chunk number 5: replace above cutoff
###################################################
qPCRBatch.taqman.replaced <- replaceAboveCutOff(qPCRBatch.taqman, 
                                                newVal = NA, cutOff = 35)
exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]


###################################################
### code chunk number 6: replace NAs with 40
###################################################
qPCRBatch.taqman.replaced <- replaceNAs(qPCRBatch.taqman, newNA = 40)
exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]


###################################################
### code chunk number 7: construct contrast matrix
###################################################
sampleNames(qPCRBatch.taqman)
a <- c(0,0,1,1,0,0,1,1) # one for each sample type, with 1 representing
b <- c(1,1,0,0,1,1,0,0) # position of sample type in samplenames vector 
contM <- cbind(a,b)
colnames(contM) <- c("case","control") # set the names of each sample type
rownames(contM) <- sampleNames(qPCRBatch.taqman) # set row names
contM
sMaxM <- t(as.matrix(c(3,3))) # now make the contrast matrix
colnames(sMaxM) <- c("case","control") # make sure these line up with samples
sMaxM


###################################################
### code chunk number 8: replace 3 or more NAs with all NAs
###################################################
qPCRBatch.taqman.replaced <- makeAllNewVal(qPCRBatch.taqman, contM, 
                                           sMaxM, newVal=NA)


###################################################
### code chunk number 9: ccl20 is now all NAs
###################################################
exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]


###################################################
### code chunk number 10: NormqPCR
###################################################
options(width = 68)
data(geNorm)
str(exprs(geNorm.qPCRBatch))


###################################################
### code chunk number 11: geNorm
###################################################
tissue <- as.factor(c(rep("BM", 9),  rep("FIB", 20), rep("LEU", 13),
                    rep("NB", 34), rep("POOL", 9)))
res.BM <- selectHKs(geNorm.qPCRBatch[,tissue == "BM"], method = "geNorm", 
                    Symbols = featureNames(geNorm.qPCRBatch), 
                    minNrHK = 2, log = FALSE)
res.POOL <- selectHKs(geNorm.qPCRBatch[,tissue == "POOL"], 
                      method = "geNorm", 
                      Symbols = featureNames(geNorm.qPCRBatch), 
                      minNrHK = 2, trace = FALSE, log = FALSE)
res.FIB <- selectHKs(geNorm.qPCRBatch[,tissue == "FIB"], 
                     method = "geNorm", 
                     Symbols = featureNames(geNorm.qPCRBatch), 
                     minNrHK = 2, trace = FALSE, log = FALSE)
res.LEU <- selectHKs(geNorm.qPCRBatch[,tissue == "LEU"], 
                     method = "geNorm", 
                     Symbols = featureNames(geNorm.qPCRBatch), 
                     minNrHK = 2, trace = FALSE, log = FALSE)
res.NB <- selectHKs(geNorm.qPCRBatch[,tissue == "NB"], 
                    method = "geNorm", 
                    Symbols = featureNames(geNorm.qPCRBatch), 
                    minNrHK = 2, trace = FALSE, log = FALSE)


###################################################
### code chunk number 12: table3
###################################################
ranks <- data.frame(c(1, 1:9), res.BM$ranking, res.POOL$ranking, 
                    res.FIB$ranking, res.LEU$ranking, 
                    res.NB$ranking)
names(ranks) <- c("rank", "BM", "POOL", "FIB", "LEU", "NB")
ranks


###################################################
### code chunk number 13: fig2
###################################################
library(RColorBrewer)
mypalette <- brewer.pal(5, "Set1")
matplot(cbind(res.BM$meanM, res.POOL$meanM, res.FIB$meanM, 
              res.LEU$meanM, res.NB$meanM), type = "b", 
        ylab = "Average expression stability M", 
        xlab = "Number of remaining control genes", 
        axes = FALSE, pch = 19, col = mypalette, 
        ylim = c(0.2, 1.22), lty = 1, lwd = 2, 
        main = "Figure 2 in Vandesompele et al. (2002)")
axis(1, at = 1:9, labels = as.character(10:2))
axis(2, at = seq(0.2, 1.2, by = 0.2), labels = seq(0.2, 1.2, by = 0.2))
box()
abline(h = seq(0.2, 1.2, by = 0.2), lty = 2, lwd = 1, col = "grey")
legend("topright", legend = c("BM", "POOL", "FIB", "LEU", "NB"), 
       fill = mypalette)


###################################################
### code chunk number 14: fig3a
###################################################
mypalette <- brewer.pal(8, "YlGnBu")
barplot(cbind(res.POOL$variation, res.LEU$variation, res.NB$variation, 
              res.FIB$variation, res.BM$variation), beside = TRUE, 
        col = mypalette, space = c(0, 2), 
        names.arg = c("POOL", "LEU", "NB", "FIB", "BM"),
        ylab = "Pairwise variation V",
        main = "Figure 3(a) in Vandesompele et al. (2002)")
legend("topright", legend = c("V9/10", "V8/9", "V7/8", "V6/7", 
                              "V5/6", "V4/5", "V3/4", "V2/3"), 
       fill = mypalette, ncol = 2)
abline(h = seq(0.05, 0.25, by = 0.05), lty = 2, col = "grey")
abline(h = 0.15, lty = 1, col = "black")


###################################################
### code chunk number 15: NormFinder
###################################################
data(Colon)
Colon
Class <- pData(Colon)[,"Classification"]
res.Colon <- stabMeasureRho(Colon, group = Class, log = FALSE)
sort(res.Colon) # see Table 3 in Andersen et al (2004)

data(Bladder)
Bladder
grade <- pData(Bladder)[,"Grade"]
res.Bladder <- stabMeasureRho(Bladder, group = grade,
                              log = FALSE)
sort(res.Bladder)


###################################################
### code chunk number 16: NormFinder1
###################################################
selectHKs(Colon, log = FALSE, trace = FALSE, 
          Symbols = featureNames(Colon))$ranking
selectHKs(Bladder, log = FALSE, trace = FALSE, 
          Symbols = featureNames(Bladder))$ranking


###################################################
### code chunk number 17: NormFinder2
###################################################
Class <- pData(Colon)[,"Classification"]
selectHKs(Colon, group = Class, log = FALSE, trace = TRUE,
          Symbols = featureNames(Colon), minNrHKs = 12,
          method = "NormFinder")$ranking


###################################################
### code chunk number 18: NormFinder3
###################################################
grade <- pData(Bladder)[,"Grade"]
selectHKs(Bladder, group = grade, log = FALSE, trace = FALSE, 
          Symbols = featureNames(Bladder), minNrHKs = 13,
          method = "NormFinder")$ranking


###################################################
### code chunk number 19: taqman read dCq
###################################################
path <- system.file("exData", package = "NormqPCR")
taqman.example <- file.path(path, "example.txt")
qPCR.example <- file.path(path, "qPCR.example.txt")
qPCRBatch.taqman <- read.taqman(taqman.example)


###################################################
### code chunk number 20: dCq
###################################################
hkgs<-"Actb-Rn00667869_m1"
qPCRBatch.norm <- deltaCq(qPCRBatch =  qPCRBatch.taqman, hkgs = hkgs, calc="arith")
head(exprs(qPCRBatch.norm))


###################################################
### code chunk number 21: dCq many genes
###################################################
hkgs<-c("Actb-Rn00667869_m1", "B2m-Rn00560865_m1", "Gapdh-Rn99999916_s1")
qPCRBatch.norm <- deltaCq(qPCRBatch =  qPCRBatch.taqman, hkgs = hkgs, calc="arith")
head(exprs(qPCRBatch.norm))


###################################################
### code chunk number 22: taqman read
###################################################
path <- system.file("exData", package = "NormqPCR")
taqman.example <- file.path(path, "example.txt")
qPCR.example <- file.path(path, "qPCR.example.txt")
qPCRBatch.taqman <- read.taqman(taqman.example)


###################################################
### code chunk number 23: contrast
###################################################
contM <- cbind(c(0,0,1,1,0,0,1,1),c(1,1,0,0,1,1,0,0))
colnames(contM) <- c("interestingPhenotype","wildTypePhenotype")
rownames(contM) <- sampleNames(qPCRBatch.taqman)
contM


###################################################
### code chunk number 24: ddCq
###################################################
hkg <- "Actb-Rn00667869_m1"
ddCq.taqman <- deltaDeltaCq(qPCRBatch = qPCRBatch.taqman, maxNACase=1, maxNAControl=1, 
                            hkg=hkg, contrastM=contM, case="interestingPhenotype", 
                            control="wildTypePhenotype", statCalc="geom", hkgCalc="arith")
head(ddCq.taqman)


###################################################
### code chunk number 25: ddCq Avg
###################################################
hkg <- "Actb-Rn00667869_m1"
ddCqAvg.taqman <- deltaDeltaCq(qPCRBatch = qPCRBatch.taqman, maxNACase=1, maxNAControl=1, 
                               hkg=hkg, contrastM=contM, case="interestingPhenotype", 
                               control="wildTypePhenotype", paired=FALSE, statCalc="geom", 
                               hkgCalc="arith")
head(ddCqAvg.taqman)


###################################################
### code chunk number 26: taqman gM
###################################################
qPCRBatch.taqman <- read.taqman(taqman.example)
contM <- cbind(c(0,0,1,1,0,0,1,1),c(1,1,0,0,1,1,0,0))
colnames(contM) <- c("interestingPhenotype","wildTypePhenotype")
rownames(contM) <- sampleNames(qPCRBatch.taqman)
hkgs<-c("Actb-Rn00667869_m1", "B2m-Rn00560865_m1", "Gapdh-Rn99999916_s1")
ddCq.gM.taqman <- deltaDeltaCq(qPCRBatch = qPCRBatch.taqman, maxNACase=1, maxNAControl=1, 
                               hkgs=hkgs, contrastM=contM, case="interestingPhenotype", 
                               control="wildTypePhenotype", statCalc="arith", hkgCalc="arith")
head(ddCq.gM.taqman)


###################################################
### code chunk number 27: taqman gM Avg
###################################################
qPCRBatch.taqman <- read.taqman(taqman.example)
contM <- cbind(c(0,0,1,1,0,0,1,1),c(1,1,0,0,1,1,0,0))
colnames(contM) <- c("interestingPhenotype","wildTypePhenotype")
rownames(contM) <- sampleNames(qPCRBatch.taqman)
hkgs<-c("Actb-Rn00667869_m1", "B2m-Rn00560865_m1", "Gapdh-Rn99999916_s1")
ddAvgCq.gM.taqman <-deltaDeltaCq(qPCRBatch = qPCRBatch.taqman, maxNACase=1, maxNAControl=1, 
                                 hkgs=hkgs, contrastM=contM, case="interestingPhenotype", 
                                 control="wildTypePhenotype", paired=FALSE, statCalc="arith", 
                                 hkgCalc="arith")
head(ddAvgCq.gM.taqman)


###################################################
### code chunk number 28: ComputeNRQs1
###################################################
path <- system.file("exData", package = "ReadqPCR")
qPCR.example <- file.path(path, "qPCR.example.txt")
Cq.data <- read.qPCR(qPCR.example)


###################################################
### code chunk number 29: ComputeNRQs2
###################################################
Cq.data1 <- combineTechRepsWithSD(Cq.data)


###################################################
### code chunk number 30: ComputeNRQs3
###################################################
Effs <- file.path(path, "Efficiencies.txt")
Cq.effs <- read.table(file = Effs, row.names = 1, header = TRUE)
rownames(Cq.effs) <- featureNames(Cq.data1)
effs(Cq.data1) <- as.matrix(Cq.effs[,"efficiency",drop = FALSE])
se.effs(Cq.data1) <- as.matrix(Cq.effs[,"SD.efficiency",drop = FALSE])


###################################################
### code chunk number 31: ComputeNRQs4
###################################################
res <- ComputeNRQs(Cq.data1, hkgs = c("gene_az", "gene_gx"))
## NRQs
exprs(res)
## SD of NRQs
se.exprs(res)