[Qpcr-commits] r100 - in pkg/NormqPCR: . R inst/doc man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Apr 15 02:09:26 CEST 2010
Author: jperkins
Date: 2010-04-15 02:09:25 +0200 (Thu, 15 Apr 2010)
New Revision: 100
Removed:
pkg/NormqPCR/R/produceHkgsDF.R
pkg/NormqPCR/man/produceHkgsDF.Rd
Modified:
pkg/NormqPCR/NAMESPACE
pkg/NormqPCR/R/deltaDeltaAvgCt.R
pkg/NormqPCR/R/deltaDeltaCt.R
pkg/NormqPCR/inst/doc/NormqPCR.Rnw
Log:
cleaned up package, added to vignette, cleaned up the ddCt methods and reintroduced ddCt by geomMean
Modified: pkg/NormqPCR/NAMESPACE
===================================================================
--- pkg/NormqPCR/NAMESPACE 2010-03-17 19:05:12 UTC (rev 99)
+++ pkg/NormqPCR/NAMESPACE 2010-04-15 00:09:25 UTC (rev 100)
@@ -1,2 +1,2 @@
importClasses(qPCRSet)
-export(geomMean,stabMeasureM,stabMeasureRho,selectHKs,deltaCt,replaceNAs,deltaDeltaCt,deltaDeltaAvgCt,produceHkgsDF,makeAllNAs,combineTechReps,replaceAboveCutOff, makeAllNAs, replaceNAs)
+export(geomMean,stabMeasureM,stabMeasureRho,selectHKs,deltaCt,replaceNAs,deltaDeltaCt,deltaDeltaAvgCt,makeAllNAs,combineTechReps,replaceAboveCutOff, makeAllNAs, replaceNAs, gM_ddCt)
Modified: pkg/NormqPCR/R/deltaDeltaAvgCt.R
===================================================================
--- pkg/NormqPCR/R/deltaDeltaAvgCt.R 2010-03-17 19:05:12 UTC (rev 99)
+++ pkg/NormqPCR/R/deltaDeltaAvgCt.R 2010-04-15 00:09:25 UTC (rev 100)
@@ -38,9 +38,15 @@
for (detector in featureNames(qPCRBatch)) {
VCase <- caseM[detector,]
VControl <- controlM[detector,]
- if(is.na(VCase)) VCase <- rep(NA, length(VCase))
+ if(! FALSE %in% is.na(VCase)) {
+ warning("No Detector for Case")
+ VCase <- rep(NA, length(VCase))
+ }
else meanCase <- geomMean(VCase, na.rm=TRUE)
- if(is.na(VControl)) VControl <- rep(NA, length(VControl))
+ if(! FALSE %in% is.na(VControl)) {
+ warning("No Detector for Control")
+ VControl <- rep(NA, length(VControl))
+ }
else meanControl <- geomMean(VControl, na.rm=TRUE)
sdCase <- sd(VCase, na.rm=TRUE)
sdControl <- sd(VControl, na.rm=TRUE)
Modified: pkg/NormqPCR/R/deltaDeltaCt.R
===================================================================
--- pkg/NormqPCR/R/deltaDeltaCt.R 2010-03-17 19:05:12 UTC (rev 99)
+++ pkg/NormqPCR/R/deltaDeltaCt.R 2010-04-15 00:09:25 UTC (rev 100)
@@ -26,13 +26,12 @@
for (detector in featureNames(qPCRBatch)) {
VCase <- caseM[detector,]
VControl <- controlM[detector,]
-warning("length")
#stop("length of case is:",VCase,"_",length(VCase))
if(length(VCase) == 1) {
- warning("Only one Detector for Case")
+ warning("Only one Detector for Control")
dCtCase <- VCase
sdCase <- NA
- } else if(is.na(VCase)) {
+ } else if(! FALSE %in% is.na(VCase)) {
warning("No Detector for Case")
dCtCase <- rep(NA, length = VCase)
dCtControl <- NA
@@ -45,7 +44,7 @@
warning("Only one Detector for Control")
dCtControl <- VControl
sdControl <- NA
- } else if(is.na(VControl)) {
+ } else if(! FALSE %in% is.na(VControl)) {
warning("No Detector for Control")
dCtControl <- rep(NA, length = VControl)
sdControl <- NA
Deleted: pkg/NormqPCR/R/produceHkgsDF.R
===================================================================
--- pkg/NormqPCR/R/produceHkgsDF.R 2010-03-17 19:05:12 UTC (rev 99)
+++ pkg/NormqPCR/R/produceHkgsDF.R 2010-04-15 00:09:25 UTC (rev 100)
@@ -1,125 +0,0 @@
-# nSet = "data.frame",
-# hkgs = "character"))
-
-###################################
-# normalise qPCRSet with 1 HKG #
-###################################
-
-produceHkgsDF <- function(qPCRSet, hkgs, design, cutoff = 40, verbose = FALSE){ # takes qPCRSet, vector of housekeeping genes and a design vector to say what comparisons to make
- hkgs <- make.names(hkgs)
- ##########
- # Use design 'matrix' to work out which is case and which is control
-cat("HHHH\n")
- logicase <- design == "case"
- logicontrol <- design == "control"
- lenCase <- sum(logicase==TRUE)
- lenControl <- sum(logicontrol==TRUE)
- maxNAinCaseHkg <- ceiling(lenCase/2) # max number of NA values allowed for case
- maxNAinControlHkg <- ceiling(lenControl/2) # max number of NA values allowed for control
- DiscardNACase <- floor(lenCase/4) # less or equal to this number of NAs we discard
- DiscardNAControl <- floor(lenControl/4) # max number of NA values allowed for control
- ##########
-cat("HHHH\n")
- normSet <- exprs(qPCRSet) # turn exprs component into a data frame
-cat("HHHH\n")
- normSet[normSet > cutoff] <- NA
- tabFormat <- vector()
- laterColNames <- vector()
-cat("HHHH\n")
- for(hkg in hkgs) { # For each nominal housekeeping gene
- if(hkg %in% featureNames(qPCRSet) == FALSE) stop (hkg," not found in file. Ensure entered housekeeping genes appear in the file")
- # first see if it's suitable
- hkgCts <- as.numeric(exprs(qPCRSet[hkg,])) # get the Ct values for given hkg
- hkgCtsCase <- hkgCts[logicase] # seperate for case and control
- hkgCtsControl <- hkgCts[logicontrol]
- # if unsuitable, stop the loop
- if(sum(is.na(hkgCtsCase)) > maxNAinCaseHkg) stop (hkg, " is unsuitable as a housekeeping gene because a value was obtained for it less than", maxNAinCaseHkg, "times out of", lenControl, ".")
- if(sum(is.na(hkgCtsControl)) > maxNAinControlHkg) stop (hkg, " is unsuitable as a housekeeping gene because a value was obtained for it less than", maxNAinControlHkg, "times out of ", lenControl, ".")
-cat("HHHHHHHHHH\n")
-
- hkg <- gsub("-.+$","",hkg) # regexp to remove any rubbish from end of control gene spec
- laterColNames <- c(laterColNames, paste(hkg, "Control_mean", sep = "_"), paste(hkg, "Control_Sds", sep = "_"), paste(hkg, "Case_mean", sep = "_"), paste(hkg, "Case_Sds", sep = "_"), paste(hkg, "ddCt", sep = "_"), paste(hkg, "2^DDCt", sep = "_"), paste(hkg, "2^DDCt min", sep = "_"), paste(hkg, "2^DDCt max", sep = "_")) # get the column names for the tSet matrix
- }
-
- for(hkg in hkgs) { # for each housekeeping gene
- if(verbose) cat("For HKG: ", hkg, "\n")
- if(hkg %in% featureNames(qPCRSet) == FALSE) stop (hkg," not found in file. Ensure entered housekeeping genes appear in the file")
- hkgCts <- as.numeric(exprs(qPCRSet[hkg, ]))
- hkgCtsCase <- hkgCts[logicase]
- hkgCtsControl <- hkgCts[logicontrol]
- ##########
- # first see if it's suitable
- hkgCts <- as.numeric(exprs(qPCRSet[hkg, ])) # get the Ct values for given hkg
- hkgCtsCase <- hkgCts[logicase] # seperate for case and control
- hkgCtsControl <- hkgCts[logicontrol]
- # if unsuitable, stop the loop
- if(sum(is.na(hkgCtsCase)) > maxNAinCaseHkg) stop (hkg, " is unsuitable as a housekeeping gene because a value was only obtained for it once or less")
- if(sum(is.na(hkgCtsControl)) > maxNAinControlHkg) stop (hkg, " is unsuitable as a housekeeping gene because a value was only obtained for it once or less")
- ##########
- # initialise vectors for different stats for the different detectors for this hkg and either case/control
- noOfFeatures <- length(featureNames(qPCRSet))
- deltaCtsMeanCase <- vector(length=noOfFeatures)
- deltaCtsSdCase <- vector(length=noOfFeatures)
- deltaCtsMeanControl <- vector(length=noOfFeatures)
- deltaCtsSdControl <- vector(length=noOfFeatures)
- deltadeltaCt <- vector(length=noOfFeatures)
- twoDDCt <- vector(length=noOfFeatures)
- minBound <- vector(length=noOfFeatures)
- maxBound <- vector(length=noOfFeatures)
- ##########
- for (detector in featureNames(qPCRSet)) {
- i_detector <- 1
- Cts <- as.numeric(exprs(qPCRSet[detector,])) # the raw values for the detector
- Cts[Cts > cutoff] <- NA
- CtsControl <- Cts[logicontrol]
- CtsCase <- Cts[logicase]
- if(sum(is.na(CtsControl)) > DiscardNAControl) { # if up to a quarter NAs
- normControl <- mean(CtsControl - hkgCtsControl, na.rm=T) # work out stats using existing values
- v_dCtsMeanControl <- mean(CtsControl - hkgCtsControl, na.rm=T)
- v_dCtsSdsControl <- sd(CtsControl - hkgCtsControl, na.rm=T)
- }
- else { # otherwise make all the NAs cutoff value and continue
- CtsControl[CtsControl == NA] <- cutoff
- hkgCtsControl[hkgCtsControl == NA] <- cutoff
- normControl <- mean(CtsControl - hkgCtsControl)
- v_dCtsMeanControl <- mean(CtsControl - hkgCtsControl)
- v_dCtsSdsControl <- sd(CtsControl - hkgCtsControl)
- }
- if(sum(is.na(CtsCase)) > DiscardNACase) { # if up to a quarter NAs
- normCase <- mean(CtsCase - hkgCtsCase,na.rm=T) # work out stats using existing values
- v_dCtsMeanCase <- mean(CtsCase - hkgCtsCase,na.rm=T)
- v_dCtsSdsCase <- sd(CtsCase - hkgCtsCase,na.rm=T)
- }
- else {
- CtsCase[CtsControl == NA] <- cutoff
- hkgCtsCase[hkgCtsControl == NA] <- cutoff
- normCase <- mean(CtsCase - hkgCtsCase)
- v_dCtsMeanCase <- mean(CtsCase - hkgCtsCase)
- v_dCtsSdsCase <- sd(CtsCase - hkgCtsCase)
- }
- ######## Now add these values to the vectors
- deltaCtsMeanCase[i_detector] <- v_dCtsMeanCase
- deltaCtsSdCase[i_detector] <- v_dCtsSdsCase
- deltaCtsMeanControl[i_detector] <- v_dCtsMeanControl
- deltaCtsSdControl[i_detector] <- v_dCtsSdsControl
-
- ddct <- normCase - normControl # log ratio
- twoDDCt[i_detector] <- 2^-ddct # fold change
-
- if(typeof(v_dCtsSdsControl) == "double") {
- maxBound[i_detector] <- 2^-(ddct - v_dCtsSdsControl)
- minBound[i_detector] <- 2^-(ddct + v_dCtsSdsControl)
- }
- else {
- minBound[i_detector] <- NA
- maxBound[i_detector] <- NA
- }
- deltadeltaCt[i_detector] <- ddct # add to vector
- i_detector <- i_detector + 1
- }
- normSet <- data.frame(normSet,deltaCtsMeanControl,deltaCtsSdControl,deltaCtsMeanCase,deltaCtsSdCase,deltadeltaCt,twoDDCt,minBound,maxBound) # add the vectors together
- }
- normSetColNames <- c(sampleNames(qPCRSet),laterColNames) # make column names
- names(normSet) <- normSetColNames
- return(normSet)
-}
Modified: pkg/NormqPCR/inst/doc/NormqPCR.Rnw
===================================================================
--- pkg/NormqPCR/inst/doc/NormqPCR.Rnw 2010-03-17 19:05:12 UTC (rev 99)
+++ pkg/NormqPCR/inst/doc/NormqPCR.Rnw 2010-04-15 00:09:25 UTC (rev 100)
@@ -1,6 +1,6 @@
%\VignetteIndexEntry{NormqPCR}
%\VignetteDepends{stats,RColorBrewer,Biobase,methods,ReadqPCR}
-%\VignetteKeywords{real-time, quantitative, PCR, housekeeper, reference gene, geNorm, NormFinder}
+%\VignetteKeywords{real-time, quantitative, PCR, housekeeper, reference gene,geNorm, NormFinder}
%\VignettePackage{NormqPCR}
%
\documentclass[11pt]{article}
@@ -11,13 +11,14 @@
pdftitle={NormqPCR: Functions for normalisation of RT-qPCR data},%
pdfauthor={Matthias Kohl and James Perkins},%
pdfsubject={NormqPCR},%
-pdfkeywords={real-time, quantitative, PCR, housekeeper, reference gene, geNorm, NormFinder},%
+pdfkeywords={real-time, quantitative, PCR, housekeeper, reference gene, geNorm,
+NormFinder},%
pagebackref,bookmarks,colorlinks,linkcolor=darkblue,citecolor=darkblue,%
pagecolor=darkblue,raiselinks,plainpages,pdftex]{hyperref}
%
\markboth{\sl Package ``{\tt NormqPCR}''}{\sl Package ``{\tt NormqPCR}''}
%
-% -------------------------------------------------------------------------------
+%------------------------------------------------------------------------------
\newcommand{\code}[1]{{\tt #1}}
\newcommand{\pkg}[1]{{\tt "#1"}}
\newcommand{\myinfig}[2]{%
@@ -28,10 +29,12 @@
\end{center}
% \end{figure}
}
-% -------------------------------------------------------------------------------
+%------------------------------------------------------------------------------
%
-% -------------------------------------------------------------------------------
+%------------------------------------------------------------------------------
+
\begin{document}
+
\SweaveOpts{keep.source = TRUE, eval = TRUE, include = FALSE}
%-------------------------------------------------------------------------------
\title{NormqPCR: Functions for normalisation of RT-qPCR data}
@@ -46,15 +49,25 @@
%-------------------------------------------------------------------------------
The package \pkg{NormqPCR} provides methods for the normalization of
real-time quantitative RT-PCR data. In this vignette we describe and
-demonstrate the available functions. Firstly we describe how the user can deal with undertermined values, and choose their own threshold for this. Then we show how the user may combine technical replicates, and deal with NA values. The rest of the vignette is split into two distinct sections, the first giving details of different methods to select the best houskeeping gene/genes for normalisation, and the second showing how to use the selected housekeeping genes to produce $2^{-\Delta Ct}$ normalised estimators and $2^{-\Delta \Delta Ct}$ estimators of differential expression.
+demonstrate the available functions. Firstly we show
+how the user may combine technical replicates, and deal with undetermined values
+and values above a user-chosen threshold. The rest of the vignette is split into
+two distinct sections, the first giving details of different methods to select
+the best houskeeping gene/genes for normalisation, and the second showing how to
+use the selected housekeeping genes to produce $2^{-\Delta Ct}$ normalised
+estimators and $2^{-\Delta \Delta Ct}$ estimators of differential expression.
%-------------------------------------------------------------------------------
-\section{Combining technical replicates together}
+\section{Combining technical replicates}
%-------------------------------------------------------------------------------
-When a raw data file is read in using read.qPCR contains technical replicates, they are dealt with by concatenating the suffix \_TechRep.n to the detector name, where
-n in {1, 2...N } is the number of the replication in the total number of replicates, N, based
+When a raw data file is read in using read.qPCR contains technical replicates,
+they are dealt with by concatenating the suffix \_TechRep.n to the detector
+name, where
+n in {1, 2...N } is the number of the replication in the total number of
+replicates, N, based
on order of appearence in the qPCR data file.
-So if we read in a file with technical replicates, we can see that the detector/feature names are thus suffixed:
+So if we read in a file with technical replicates, we can see that the
+detector/feature names are thus suffixed:
<<read.qPCR.tech.reps>>=
library(ReadqPCR) # load the ReadqPCR library
@@ -65,7 +78,10 @@
rownames(exprs(qPCRBatch.qPCR.techReps))[1:8]
@
-It is likely that before continuing with the analysis, the user would wish to average the technical replicates by finding their geometric mean. This can be achieved using the combineTechReps function, which will produce a new qPCRBatch object, with all tech reps reduced to one reading:
+It is likely that before continuing with the analysis, the user would wish to
+average the technical replicates by finding their geometric mean. This can be
+achieved using the combineTechReps function, which will produce a new qPCRBatch
+object, with all tech reps reduced to one reading:
<<combine read.qPCR.tech.reps>>=
combinedTechReps <- combineTechReps(qPCRBatch.qPCR.techReps)
@@ -76,65 +92,95 @@
\section{Dealing with undetermined values}
%-------------------------------------------------------------------------------
-When an RT-qPCR experiment does not produce a reading after a certain number of cycles (cycle threshold), the reading is given as undetermined. These are represented in qPCRBatch objects as "NA". Different users may have different ideas about how many cycles they wish to allow before declaring a detector as not present in the sample. There are two methods for the user to decide what to do with numbers above a given cycle threshold:
+When an RT-qPCR experiment does not produce a reading after a certain number of
+cycles (cycle threshold), the reading is given as undetermined. These are
+represented in qPCRBatch objects as "NA". Different users may have different
+ideas about how many cycles they wish to allow before declaring a detector as
+not present in the sample. There are two methods for the user to decide what to
+do with numbers above a given cycle threshold:
-First the user might decide that anything above 38 cycles means there is nothing present in their sample, instead of the standard 40 used by the taqman software. They can replace the value of all readings above 38 as NA using the following:
+First the user might decide that anything above 38 cycles means there is nothing
+present in their sample, instead of the standard 40 used by the taqman software.
+They can replace the value of all readings above 38 as NA using the following:
-Firstly read in the taqman example file from \pkg{ReadqPCR} which has 96 detectors, with 4 replicates for mia (case) and 4 non-mia (control):
+Firstly read in the taqman example file from \pkg{ReadqPCR} which has 96
+detectors, with 4 replicates for mia (case) and 4 non-mia (control):
<<taqman read>>=
path <- system.file("exData", package = "ReadqPCR")
taqman.example <- paste(path, "/example.txt", sep="")
qPCRBatch.taqman <- read.taqman(taqman.example)
@
-We can see that for the detector: Ccl20.Rn00570287\_m1 have these readings for the different samples:
+We can see that for the detector: Ccl20.Rn00570287\_m1 have these readings for
+the different samples:
<<taqman detector example>>=
exprs(qPCRBatch.taqman)["Ccl20.Rn00570287_m1",]
@
-Now use the \code{replaceAboveCutOff} method in order to replace anything above 35 with NA:
+Now use the \code{replaceAboveCutOff} method in order to replace anything above
+35 with NA:
<<replace above cutoff>>=
-qPCRBatch.taqman.replaced <- replaceAboveCutOff(qPCRBatch.taqman, newVal = NA, cutOff = 35)
+qPCRBatch.taqman.replaced <- replaceAboveCutOff(qPCRBatch.taqman, newVal = NA,
+cutOff = 35)
exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]
@
-It may also be the case that the user wants to get rid of all NA values, and replace them with an arbitrary number. This can be done using the \code{replaceNAs} method. So if the user wanted to replace all NAs with 40, it can be done as follows:
+It may also be the case that the user wants to get rid of all NA values, and
+replace them with an arbitrary number. This can be done using the
+\code{replaceNAs} method. So if the user wanted to replace all NAs with 40, it
+can be done as follows:
<<replace NAs with 40>>=
qPCRBatch.taqman.replaced <- replaceNAs(qPCRBatch.taqman, newNA = 40)
exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]
@
-In addition, the situation sometimes arises where some readings for a given detector are above a given cycle threshold, but some others are not. The user may decide for example that if a given number of readings are NAs, then all of the readings for this detector should be NAs. This is important because otherwise an unusual reading for one detector might lead to an inaccurate estimate for the expression of a given gene.
+In addition, the situation sometimes arises where some readings for a given
+detector are above a given cycle threshold, but some others are not. The user
+may decide for example that if a given number of readings are NAs, then all of
+the readings for this detector should be NAs. This is important because
+otherwise an unusual reading for one detector might lead to an inaccurate
+estimate for the expression of a given gene.
-This process will necessarily be separate for the different sample types, since you might expect a given gene to show expression in one sample type compared to another. Therefore it is necessary to designate the replicates per sample type using a contrast matrix. It is also necessary to make a sampleMaxMatrix which gives a maximum number of NAs allowed for each sample type.
+This process will necessarily be separate for the different sample types, since
+you might expect a given gene to show expression in one sample type compared to
+another. Therefore it is necessary to designate the replicates per sample type
+using a contrast matrix. It is also necessary to make a sampleMaxMatrix which
+gives a maximum number of NAs allowed for each sample type.
-So in the example file above we two sample types, with 4 biological replicates for each. the contrastMatrix and sampleMaxMatrix might be contructed like this:
+So in the example file above we two sample types, with 4 biological replicates
+for each. the contrastMatrix and sampleMaxMatrix might be contructed like this:
<<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) # the position of the sample type in the samplenames vector
+b <- c(1,1,0,0,1,1,0,0) # position of sample type in the samplenames vector
contM <- cbind(a,b)
colnames(contM) <- c("case","control") # then give the names of each sample type
rownames(contM) <- sampleNames(qPCRBatch.taqman) # and the rows of the matrix
contM
sMaxM <- t(as.matrix(c(3,3))) # now make the contrast matrix
-colnames(sMaxM) <- c("case","control") # make sure these line up with samplenames
+colnames(sMaxM) <- c("case","control") # make sure these line up with samples
sMaxM
@
More details on the contrast matrix can be found in the limma manual.
-For example, if the user decides that if at least 3 out of 4 readings are NAs for a given detector, then all readings should be NA, they will can do the following, using the \code{makeAllNAs} method:
+For example, if the user decides that if at least 3 out of 4 readings are NAs
+for a given detector, then all readings should be NA, they will can do the
+following, using the \code{makeAllNAs} method:
<<replace 3 or more NAs with all NAs>>=
qPCRBatch.taqman.replaced <- makeAllNAs(qPCRBatch.taqman, contM, sMaxM)
@
-Here you can see for the Ccl20.Rn00570287\_m1 detector, the control values have been made all NA, wheras before 3 were NA and one was 35. However the case values have been kept, since they were all below the NA threshold. It is important to filter the data in this way to ensure the correct calculations are made downstream.
+Here you can see for the Ccl20.Rn00570287\_m1 detector, the control values have
+been made all NA, wheras before 3 were NA and one was 35. However the case
+values have been kept, since they were all below the NA threshold. It is
+important to filter the data in this way to ensure the correct calculations are
+made downstream.
<< ccl20 is now all NAs >>=
exprs(qPCRBatch.taqman.replaced)["Ccl20.Rn00570287_m1",]
@@ -143,11 +189,16 @@
%-------------------------------------------------------------------------------
\section{Selection of most stable reference/housekeeping genes}
%-------------------------------------------------------------------------------
+This section contains two subsections containing different methods for the
+selection of appropriate housekeeping genes.
+
%-------------------------------------------------------------------------------
\subsection{geNorm}
%-------------------------------------------------------------------------------
-We describe the selection of the best (most stable) reference/housekeeping genes
-using the method of Vandesompele et al (2002)~\cite{geNorm} (in the sequel: Vand02)
+We describe the selection of the best (most stable) reference/housekeeping
+genes
+using the method of Vandesompele et al (2002)~\cite{geNorm} (in the sequel:
+Vand02)
which is called {\it geNorm}. We first load the package and the data
<<NormqPCR>>=
options(width = 68)
@@ -217,7 +268,8 @@
@
\myinfig{1}{NormqPCR-fig2.pdf}
\par
-Second, we plot the pairwise variation for each cell type (cf. Figure~3~(a) in Vand02)
+Second, we plot the pairwise variation for each cell type (cf. Figure~3~(a) in
+Vand02)
<<fig3a, fig = TRUE>>=
mypalette <- brewer.pal(8, "YlGnBu")
barplot(cbind(res.POOL$variation, res.LEU$variation, res.NB$variation,
@@ -241,14 +293,16 @@
\subsection{NormFinder}
%-------------------------------------------------------------------------------
The second method for selection reference/housekeeping genes implemented in
-package is the method derived by \cite{NormFinder} (in the sequel: And04) called
+package is the method derived by \cite{NormFinder} (in the sequel: And04)
+called
{\it NormFinder}.\\
The ranking contained in Table~3 of And04 can be obtained via
<<NormFinder>>=
Colon.example <- paste(path, "/qPCR.colon.txt", sep="")
Colon.qPCRBatch <- read.qPCR(Colon.example)
str(exprs(Colon.qPCRBatch))
-Colon.qPCRBatch[["Group"]] <- c(rep("Normal",10),rep("Dukes A",10),rep("Dukes B",10),rep("Dukes C",10))
+Colon.qPCRBatch[["Group"]] <- c(rep("Normal",10),rep("Dukes A",10),rep("Dukes
+B",10),rep("Dukes C",10))
pData(Colon.qPCRBatch)
res.Colon <- stabMeasureRho(Colon.qPCRBatch,
log = FALSE)
@@ -273,8 +327,10 @@
Symbols = featureNames(Bladder.qPCRBatch))$ranking
@
As we are often interested in more than one reference/housekeeping gene we also
-implemented a step-wise procedure of the NormFinder algorithm explained in Section
-``Average control gene'' in the supplementary information of And04. This procedure
+implemented a step-wise procedure of the NormFinder algorithm explained in
+Section
+``Average control gene'' in the supplementary information of And04. This
+procedure
is available via function \code{selectHKs}.
<<NormFinder2>>=
selectHKs(Colon.qPCRBatch,
@@ -289,22 +345,37 @@
%-------------------------------------------------------------------------------
\section{Normalization by means of reference/housekeeping genes}
%-------------------------------------------------------------------------------
-
\subsection{$2^{-\Delta \Delta Ct}$ method using a single housekeeper}
+%-------------------------------------------------------------------------------
+It is possible to use the $2^{-\Delta \Delta Ct}$ method for calculating the
+relative gene expression.
+Both the same well and the separate well methods as detailed in \cite{ddCt} can
+be used for this purpose.
+They have been named \code{deltaDeltaCt} and \code{deltaDeltaAvgCt}
+respectively, with \code{deltaDeltaAvgCt} so named since it calculates the
+average standard deviation between case and control as s = $\sqrt{s_{1}^{2} +
+s_{2}^{2}}$.
+This approach is not recommended when the housekeeper and genes to be compared
+are from the same sample, as is the case when using the taqman cards, but is
+included for compleness and for situations where readings for the housekeeper
+might be taken from a separate biological replicate for example in a {\it post
+hoc} manner due to the originally designated housekeeping genes not performing
+well, or for when NormqPCR is used for more traditional qPCR where the products
+undergo amplifications from separate wells.
-It is possible to use the $2^{-\Delta \Delta Ct}$ method for calculating the relative gene expression.
-Both the same well and the separate well methods as detailed in \cite{ddCt} can be used for this purpose.
-They have been named \code{deltaDeltaCt} and \code{deltaDeltaAvgCt} respectively, with \code{deltaDeltaAvgCt} so named since it calculates the average standard deviation between case and control as s = $\sqrt{s_{1}^{2} + s_{2}^{2}}$.
-This approach is not recommended when the housekeeper and genes to be compared are from the same sample, as is the case when using the taqman cards, but is included for compleness and for situations where readings for the housekeeper might be taken from a separate biological replicate for example in a {\it post hoc} manner due to the originally designated housekeeping genes not performing well, or for when NormqPCR is used for amplifications from separate wells.
-
-for the example taqman dataset from \pkg{ReadqPCR} we must first read in the data:
+for the example taqman dataset from \pkg{ReadqPCR} we must first read in the
+data:
<<taqman read>>=
path <- system.file("exData", package = "ReadqPCR")
taqman.example <- paste(path, "/example.txt", sep="")
qPCRBatch.taqman <- read.taqman(taqman.example)
@
-\code{deltaDeltaCt} and \code{deltaDeltaAvgCt} also require a contrast matrix. This is to contain columns which will be used to specify the samples representing \code{case} and \code{control} which are to be compared, in a similar way to the \pkg{limma} package. these columns should contain 1s or 0s which refer to the samples in either category:
+\code{deltaDeltaCt} and \code{deltaDeltaAvgCt} also require a contrast matrix.
+This is to contain columns which will be used to specify the samples
+representing \code{case} and \code{control} which are to be compared, in a
+similar way to the \pkg{limma} package. these columns should contain 1s or 0s
+which refer to the samples in either category:
<< contrast >>=
contM <- cbind(c(0,0,0,0,1,1,1,1),c(1,1,1,1,0,0,0,0))
@@ -313,35 +384,54 @@
contM
@
-We can now normalise each sample by a given housekeeping gene and then look at the ratio of expression between the case and control samples. Results show the difference between the mean value for case and control following subtraction of the housekeeping genes (first column) followed by the range of values that correspond to 1 s.d. either side of the mean value, as detailed in \cite{ddCt}
+We can now normalise each sample by a given housekeeping gene and then look at
+the ratio of expression between the case and control samples. Results show the
+difference between the mean value for case and control following subtraction of
+the housekeeping genes (first column) followed by the range of values that
+correspond to 1 s.d. either side of the mean value, as detailed in \cite{ddCt}
<< ddCt >>=
hkg <- "Actb-Rn00667869_m1"
-ddCt.taqman <- deltaDeltaCt(qPCRBatch.taqman,1,1,hkg, contM, "interestingPhenotype","wildTypePhenotype")
+ddCt.taqman <- deltaDeltaCt(qPCRBatch.taqman,1,1,hkg, contM,
+"interestingPhenotype","wildTypePhenotype")
head(ddCt.taqman)
@
-We can also average the taqman data using the separate samples/wells method \code{deltaDeltaAvgCt}.
-Note how although the values are the same the ranges (+/- 1 s.d.) are generally a little higher.
-This is because the averages for case or control are calculated, and the average values for the respective housekeepers are calculated independently and then subtracted. Therefore the pairing of housekeeper with the detector value within the same sample is lost.
+We can also average the taqman data using the separate samples/wells method
+\code{deltaDeltaAvgCt}.
+Note how although the values are the same the ranges (+/- 1 s.d.) are generally
+a little higher.
+This is because the averages for case or control are calculated, and the average
+values for the respective housekeepers are calculated independently and then
+subtracted. Therefore the pairing of housekeeper with the detector value within
+the same sample is lost, increasing the variance.
<< ddCt Avg >>=
hkg <- "Actb-Rn00667869_m1"
-ddCtAvg.taqman <- deltaDeltaAvgCt(qPCRBatch.taqman,1,1,hkg, contM, "interestingPhenotype","wildTypePhenotype")
+ddCtAvg.taqman <- deltaDeltaAvgCt(qPCRBatch.taqman,1,1,hkg, contM,
+"interestingPhenotype","wildTypePhenotype")
head(ddCtAvg.taqman)
@
+%-------------------------------------------------------------------------------
+\subsection{$2^{\Delta \Delta Ct}$ method) using geometric mean of a number of
+housekeeping genes}
+%-------------------------------------------------------------------------------
+If the user wishes to normalise by more than one housekeeping gene, for example
+if they have found a more than one housekeeping gene using the NormFinder/geNorm
+algorithms described above, they can. This is implemented by calculating the
+geometric
+mean of these values to form a "pseudo-housekeeper" from which the other
+values are subtracted. We illustrate this using the data from the geNorm dataset
+above.
-\subsection{$2^{\Delta \Delta Ct}$ method) using geometric mean of a number of housekeeping genes}
-If the user wishes to normalise by more than one housekeeping gene, for example if they have found a more than one housekeeping gene using the NormFinder/geNorm algorithms described above, and wish to combine them they can use the geometric mean of these values as a pseudo-housekeeper from which to subtract the other values. We illustrate this using the data from the geNorm dataset above.
-
-
%-------------------------------------------------------------------------------
\begin{thebibliography}{1}
\bibitem{NormFinder}
Claus Lindbjerg Andersen, Jens Ledet Jensen and Torben Falck Orntoft (2004).
-\newblock Normalization of Real-Time Quantitative Reverse Transcription-PCR Data:
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/qpcr -r 100
More information about the Qpcr-commits
mailing list