[Qpcr-commits] r102 - in pkg: NormqPCR NormqPCR/R NormqPCR/inst/doc NormqPCR/man QCqPCR/R ReadqPCR/R
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun May 16 19:23:06 CEST 2010
Author: jperkins
Date: 2010-05-16 19:23:06 +0200 (Sun, 16 May 2010)
New Revision: 102
Modified:
pkg/NormqPCR/NAMESPACE
pkg/NormqPCR/R/deltaDeltaCt.R
pkg/NormqPCR/inst/doc/NormqPCR.Rnw
pkg/NormqPCR/man/deltaDeltaCt.Rd
pkg/QCqPCR/R/qPCRPairs.R
pkg/ReadqPCR/R/qPCRbatch.R
pkg/ReadqPCR/R/taqmanbatch.R
Log:
few changes..nearly there
Modified: pkg/NormqPCR/NAMESPACE
===================================================================
--- pkg/NormqPCR/NAMESPACE 2010-05-13 21:31:12 UTC (rev 101)
+++ pkg/NormqPCR/NAMESPACE 2010-05-16 17:23:06 UTC (rev 102)
@@ -1,2 +1,2 @@
importClasses(qPCRSet)
-export(geomMean,stabMeasureM,stabMeasureRho,selectHKs,deltaCt,replaceNAs,deltaDeltaCt,deltaDeltaAvgCt,makeAllNAs,combineTechReps,replaceAboveCutOff, makeAllNAs, replaceNAs, gM_ddCt, gM_ddAvgCt)
+export(geomMean,stabMeasureM,stabMeasureRho,selectHKs,deltaCt,replaceNAs,deltaDeltaCt,deltaDeltaAvgCt,makeAllNAs,combineTechReps,replaceAboveCutOff, makeAllNAs, replaceNAs, gM_ddAvgCt)
Modified: pkg/NormqPCR/R/deltaDeltaCt.R
===================================================================
--- pkg/NormqPCR/R/deltaDeltaCt.R 2010-05-13 21:31:12 UTC (rev 101)
+++ pkg/NormqPCR/R/deltaDeltaCt.R 2010-05-16 17:23:06 UTC (rev 102)
@@ -1,26 +1,50 @@
setGeneric("deltaDeltaCt",
- function(qPCRBatch, maxNACase=0, maxNAControl=0, hkg, contrastM, case, control, paired=TRUE)
+ function(qPCRBatch, maxNACase=0, maxNAControl=0, hkgs, contrastM, case, control, paired=TRUE, combineHkg=FALSE)
standardGeneric("deltaDeltaCt")
)
setMethod("deltaDeltaCt", signature = "qPCRBatch", definition =
- function(qPCRBatch, maxNACase, maxNAControl, hkg, contrastM, case, control, paired) {
- hkg <- make.names(hkg)
- if(! hkg %in% featureNames(qPCRBatch)) stop("invalid housekeeping gene")
- if(sum(is.na(hkg)) > 0) warning(hkg, " May be a bad housekeeping gene to normalise with since it did not produce a reading ", sum(is.na(hkg)), "times out of", length(hkg), ". deltaDeltaAvgCt() function might give more robust results")
-
+ function(qPCRBatch, maxNACase, maxNAControl, hkgs, contrastM, case, control, paired, combineHkg) {
+ hkgs <- make.names(hkgs)
+ if(combineHkg == TRUE) {
+ if(length(hkgs) == 1) stop("Not enough hkgs given")
+ }
+# for(hkg in hkgs) {
+# if(! hkg %in% featureNames(qPCRBatch)) stop("invalid housekeeping gene")
+ if(FALSE %in% (hkgs %in% featureNames(qPCRBatch))) stop("invalid housekeeping gene given")
+ for(hkg in hkgs){
+ if(sum(is.na(hkg)) > 0) warning(hkg, " May be a bad housekeeping gene to normalise with since it did not produce a reading ", sum(is.na(hkg)), "times out of", length(hkg))
+ }
case <- row.names(contrastM)[contrastM[,case] == 1]
control <- row.names(contrastM)[contrastM[,control] == 1]
expM <- exprs(qPCRBatch)
caseM <- expM[,case]
controlM <- expM[,control]
+
+# hkgMCase <- caseM[hkgs, ]
+# hkgMControl <- controlM[hkgs, ]
+ if(combineHkg == TRUE) {
+ hkgMCase <- caseM[hkgs, ]
+ hkgMControl <- controlM[hkgs, ]
+ hkgVCase <- apply(hkgMCase, 2, geomMean, na.rm=TRUE)
+ hkgVControl <- apply(hkgMControl, 2, geomMean, na.rm=TRUE)
+ } else {
+ hkg <- hkgs[1]
+ }
+
hkgVCase <- caseM[hkg, ]
hkgVControl <- controlM[hkg, ]
- if(! FALSE %in% is.na(hkgVCase) || ! FALSE %in% is.na(hkgVControl)) stop("Need at least 1 non NA for the housekeeper")
sdHkgCase <- sd(hkgVCase, na.rm=TRUE)
+ sdHkgControl <- sd(hkgVControl, na.rm=TRUE)
+ if(! FALSE %in% is.na(hkgVCase) || ! FALSE %in% is.na(hkgVControl)) stop("Need at least 1 non NA for the housekeeper")
+
ddCts <- vector(length=length(featureNames(qPCRBatch)))
+ dCtCases <- vector(length=length(featureNames(qPCRBatch)))
+ dCtControls <- vector(length=length(featureNames(qPCRBatch)))
+ sdCtCases <- vector(length=length(featureNames(qPCRBatch)))
+ sdCtControls <- vector(length=length(featureNames(qPCRBatch)))
minddCts <- vector(length=length(ddCts))
maxddCts <- vector(length=length(ddCts))
@@ -54,6 +78,11 @@
dCtControl <- rep(NA, length = VControl)
} else {
dCtControl <- geomMean(VControl, na.rm=TRUE) - geomMean(hkgVControl, na.rm=TRUE)
+ if (paired == TRUE) {
+ sdControl <- sd(VControl - hkgVControl, na.rm=TRUE)
+ } else {
+ sdControl <- sqrt(sd(VControl, na.rm=TRUE)^2 + sdHkgControl^2)
+ }
}
if(sum(is.na(VCase)) > maxNACase) {
dCtCase <- NA
@@ -83,8 +112,12 @@
ddCts[i] <- 2^-ddCt
}
}
+ dCtCases[i] <- 2^-dCtCase
+ sdCtCases[i] <- sdCase
+ dCtControls[i] <- 2^-dCtControl
+ sdCtControls[i] <- sdControl
i <- i+1
}
- return(cbind(featureNames(qPCRBatch),ddCts,minddCts,maxddCts))
+ return(cbind(featureNames(qPCRBatch),dCtCases,sdCtCases,dCtControls,sdCtControls,ddCts,minddCts,maxddCts))
}
)
Modified: pkg/NormqPCR/inst/doc/NormqPCR.Rnw
===================================================================
--- pkg/NormqPCR/inst/doc/NormqPCR.Rnw 2010-05-13 21:31:12 UTC (rev 101)
+++ pkg/NormqPCR/inst/doc/NormqPCR.Rnw 2010-05-16 17:23:06 UTC (rev 102)
@@ -39,8 +39,8 @@
%-------------------------------------------------------------------------------
\title{NormqPCR: Functions for normalisation of RT-qPCR data}
%-------------------------------------------------------------------------------
-\author{Matthias Kohl\\
-University of Bayreuth (Germany)\medskip\\
+\author{James Perkins and Matthias Kohl\\
+University College London (UK) / University of Bayreuth (Germany)\medskip\\
}
\maketitle
\tableofcontents
@@ -397,13 +397,8 @@
@
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.
+\code{deltaDeltaAvgCt}. Here sd is calculated separately for the sample of interest and the sample used for callibration (i.e. case versus control) and combined Therefore the pairing of housekeeper with the detector value within
+the same sample is lost, in some cases causing increased variance.
<< ddCt Avg >>=
hkg <- "Actb-Rn00667869_m1"
@@ -428,8 +423,8 @@
colnames(contM) <- c("interestingPhenotype","wildTypePhenotype")
rownames(contM) <- sampleNames(qPCRBatch.taqman)
hkgs<-c("Actb-Rn00667869_m1", "B2m-Rn00560865_m1", "Gapdh-Rn99999916_s1")
-ddCt.gM.taqman <- gM_ddCt(qPCRBatch = qPCRBatch.taqman, maxNACase=1, maxNAControl=1, hkgs=hkgs, contrastM=contM, case="interestingPhenotype", control="wildTypePhenotype")
-head(ddCt.gM.taqman)
+#ddCt.gM.taqman <- gM_ddCt(qPCRBatch = qPCRBatch.taqman, maxNACase=1, maxNAControl=1, hkgs=hkgs, contrastM=contM, case="interestingPhenotype", control="wildTypePhenotype")
+#head(ddCt.gM.taqman)
@
There is also the option of using the geometric mean method using shared variance between the samples being compared, similar
Modified: pkg/NormqPCR/man/deltaDeltaCt.Rd
===================================================================
--- pkg/NormqPCR/man/deltaDeltaCt.Rd 2010-05-13 21:31:12 UTC (rev 101)
+++ pkg/NormqPCR/man/deltaDeltaCt.Rd 2010-05-16 17:23:06 UTC (rev 102)
@@ -8,7 +8,7 @@
Suitable when housekeeping genes are from same wells/sample as the other detectors
}
\usage{
-deltaDeltaCt(qPCRBatch, maxNACase=0, maxNAControl=0, hkg, contrastM, case, control)
+deltaDeltaCt(qPCRBatch, maxNACase=0, maxNAControl=0, hkgs, contrastM, case, control, paired=TRUE, combineHkg=FALSE)
}
\arguments{
\item{qPCRBatch}{ Expression set containing qPCR data.
@@ -17,7 +17,7 @@
}
\item{maxNAControl}{ Maximum number of NA values allowed before a detector's reading is discarded for samples designated as control
}
- \item{hkg}{ String containing the name of th name of the housekeeping gene which will be used to normalise the rest of the genes.
+ \item{hkgs}{ String containing the name of th name of the housekeeping gene which will be used to normalise the rest of the genes.
}
\item{contrastM}{ A binary matrix which designates case and control samples
}
@@ -25,12 +25,16 @@
}
\item{control}{ The name of the column in contrastM that corresponds to the control samples
}
+ \item{paired}{ Logical - if TRUE the detectors and housekeepers in the same sample will be paired for calculating standard deviation, effectively meaning we will be calculating standard deviation of the differences. If FALSE, there will be no pairing, and standard deviation will be pooled between the detector and housekeepers
}
+ \item{combineHkgs}{ Logical - if TRUE, then as long as more than one housekeeper given for argument hkgs, it will combine the housekeepers by finding the geometric mean. Housekeepers can be found using geNorm or NormFinder algorithms
+}
+}
\details{
Takes expression set of qPCR values and normalises them using different housekeeping genes. Returns seperate sets of values for each housekeeping gene given.
}
\value{
- matrix with columns containing the detector ids, the 2^delta delta Ct values and the minimum and maximum values (ie the values that are 1 sd away )
+ matrix with columns containing the detector ids, 2^delta Ct values for the sample of interest and the callibrator sample, alongside their respective standard deviations, the 2^delta delta Ct values and the minimum and maximum values (ie the values that are 1 sd away )
}
%\references{ }
\author{ James Perkins \email{jperkins at biochem.ucl.ac.uk}}
Modified: pkg/QCqPCR/R/qPCRPairs.R
===================================================================
--- pkg/QCqPCR/R/qPCRPairs.R 2010-05-13 21:31:12 UTC (rev 101)
+++ pkg/QCqPCR/R/qPCRPairs.R 2010-05-16 17:23:06 UTC (rev 102)
@@ -5,8 +5,7 @@
setMethod("qPCRPairs", signature = "qPCRBatch", definition =
function(qPCRBatch, plotType, writeToFile, pairsToPlot) {
if(plotType == "Sample") {
-#cat("NOTHEREIHOPE\n")
- if(pairsToPlot == "All") {
+ if(pairsToPlot == "All") {
pairsToPlot <- combn(sampleNames(qPCRBatch),2)
}
else {
@@ -21,9 +20,6 @@
pairsToPlot <- combn(sampleNames(qPCRBatch),2)
}
else {
-
- }
-
plateVec <- as.vector(gsub("-.*", "", orderMat))
wellVec <- as.numeric(gsub(".*-", "", orderMat))
plotMat <- matrix(ncol = length(levels(as.factor(plateVec))), nrow = max(wellVec))
@@ -37,7 +33,7 @@
apply(pairsToPlot, 2, .plotPairs, plotMat, writeToFile)
}
else stop("incorrect plotType argument given")
- }
+}
)
.plotPairs <- function(samples, plotMat, writeToFile) # plots graph between the 2 samples
{
Modified: pkg/ReadqPCR/R/qPCRbatch.R
===================================================================
--- pkg/ReadqPCR/R/qPCRbatch.R 2010-05-13 21:31:12 UTC (rev 101)
+++ pkg/ReadqPCR/R/qPCRbatch.R 2010-05-16 17:23:06 UTC (rev 102)
@@ -11,7 +11,7 @@
n <- length(colnames(exprs))
if (dim(pdata)[1] != n) { # so if we don't have a row for each sample in the pData matrix
warning("Incompatible phenoData object. Created a new one using sample name data derived from raw data.\n")
- samplenames <- sub("^/?([^/]*/)*", "", colnames(exprs), extended = TRUE)
+ samplenames <- sub("^/?([^/]*/)*", "", colnames(exprs))
pdata <- data.frame(sample = 1:length(samplenames), row.names = samplenames)
phenoData <- new("AnnotatedDataFrame", data = pdata,
varMetadata = data.frame(labelDescription = "arbitrary numbering",
@@ -33,7 +33,7 @@
raw.data <- read.table(filename, header=TRUE)
if(is.null(raw.data$Well) || is.null(raw.data$PlateID)) {
noWellData <- TRUE
- if (verbose) cat("No Well and/or Plate info found, skipping this part")
+ if (verbose) cat("No Well and/or Plate info found, skipping this part", "\n")
}
else {
raw.data$PlateID <- paste(raw.data$PlateID, as.character(raw.data$Well), sep= "-")
Modified: pkg/ReadqPCR/R/taqmanbatch.R
===================================================================
--- pkg/ReadqPCR/R/taqmanbatch.R 2010-05-13 21:31:12 UTC (rev 101)
+++ pkg/ReadqPCR/R/taqmanbatch.R 2010-05-16 17:23:06 UTC (rev 102)
@@ -40,7 +40,7 @@
n <- length(colnames(exprs))
if (dim(pdata)[1] != n) { # so if we don't have a row for each sample in the pData matrix
warning("Incompatible phenoData object. Created a new one using sample name data derived from raw data.\n")
- samplenames <- sub("^/?([^/]*/)*", "", colnames(exprs), extended = TRUE)
+ samplenames <- sub("^/?([^/]*/)*", "", colnames(exprs))
pdata <- data.frame(sample = 1:length(samplenames), row.names = samplenames)
phenoData <- new("AnnotatedDataFrame", data = pdata,
varMetadata = data.frame(labelDescription = "arbitrary numbering",
More information about the Qpcr-commits
mailing list