[Qpcr-commits] r101 - in pkg/NormqPCR: . R inst/doc
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu May 13 23:31:12 CEST 2010
Author: jperkins
Date: 2010-05-13 23:31:12 +0200 (Thu, 13 May 2010)
New Revision: 101
Added:
pkg/NormqPCR/R/gM_ddAvgCt.R
Modified:
pkg/NormqPCR/NAMESPACE
pkg/NormqPCR/R/deltaDeltaCt.R
pkg/NormqPCR/inst/doc/NormqPCR.Rnw
Log:
changed delta delta cd - now calculates s.d. properly, as per the shmittigan methods paper (2001) also now this can do same and separate well calculations
Modified: pkg/NormqPCR/NAMESPACE
===================================================================
--- pkg/NormqPCR/NAMESPACE 2010-04-15 00:09:25 UTC (rev 100)
+++ pkg/NormqPCR/NAMESPACE 2010-05-13 21:31:12 UTC (rev 101)
@@ -1,2 +1,2 @@
importClasses(qPCRSet)
-export(geomMean,stabMeasureM,stabMeasureRho,selectHKs,deltaCt,replaceNAs,deltaDeltaCt,deltaDeltaAvgCt,makeAllNAs,combineTechReps,replaceAboveCutOff, makeAllNAs, replaceNAs, gM_ddCt)
+export(geomMean,stabMeasureM,stabMeasureRho,selectHKs,deltaCt,replaceNAs,deltaDeltaCt,deltaDeltaAvgCt,makeAllNAs,combineTechReps,replaceAboveCutOff, makeAllNAs, replaceNAs, gM_ddCt, gM_ddAvgCt)
Modified: pkg/NormqPCR/R/deltaDeltaCt.R
===================================================================
--- pkg/NormqPCR/R/deltaDeltaCt.R 2010-04-15 00:09:25 UTC (rev 100)
+++ pkg/NormqPCR/R/deltaDeltaCt.R 2010-05-13 21:31:12 UTC (rev 101)
@@ -1,9 +1,9 @@
setGeneric("deltaDeltaCt",
- function(qPCRBatch, maxNACase=0, maxNAControl=0, hkg, contrastM, case, control)
+ function(qPCRBatch, maxNACase=0, maxNAControl=0, hkg, contrastM, case, control, paired=TRUE)
standardGeneric("deltaDeltaCt")
)
setMethod("deltaDeltaCt", signature = "qPCRBatch", definition =
- function(qPCRBatch, maxNACase, maxNAControl, hkg, contrastM, case, control) {
+ 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")
@@ -18,6 +18,8 @@
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)
+
ddCts <- vector(length=length(featureNames(qPCRBatch)))
minddCts <- vector(length=length(ddCts))
maxddCts <- vector(length=length(ddCts))
@@ -37,17 +39,19 @@
dCtControl <- NA
} else {
dCtCase <- geomMean(VCase, na.rm=TRUE) - geomMean(hkgVCase, na.rm=TRUE)
- sdCase <- sd(VCase - hkgVCase, na.rm=TRUE)
+ if (paired == TRUE) {
+ sdCase <- sd(VCase - hkgVCase, na.rm=TRUE)
+ } else {
+ sdCase <- sqrt(sd(VCase, na.rm=TRUE)^2 + sdHkgCase^2)
+ }
}
if(length(VControl) == 1) {
warning("Only one Detector for Control")
dCtControl <- VControl
- sdControl <- NA
} else if(! FALSE %in% is.na(VControl)) {
warning("No Detector for Control")
dCtControl <- rep(NA, length = VControl)
- sdControl <- NA
} else {
dCtControl <- geomMean(VControl, na.rm=TRUE) - geomMean(hkgVControl, na.rm=TRUE)
}
@@ -73,9 +77,9 @@
maxddCts[i] <- NA
}
else {
- sdCt <- sqrt(sdCase)
- minddCts[i] <- 2 ^ -(ddCt + sdCt)
- maxddCts[i] <- 2 ^ -(ddCt - sdCt)
+# sdCt <- sqrt(sdCase)
+ minddCts[i] <- 2 ^ -(ddCt + sdCase)
+ maxddCts[i] <- 2 ^ -(ddCt - sdCase)
ddCts[i] <- 2^-ddCt
}
}
Added: pkg/NormqPCR/R/gM_ddAvgCt.R
===================================================================
--- pkg/NormqPCR/R/gM_ddAvgCt.R (rev 0)
+++ pkg/NormqPCR/R/gM_ddAvgCt.R 2010-05-13 21:31:12 UTC (rev 101)
@@ -0,0 +1,104 @@
+setGeneric("gM_ddAvgCt",
+ function(qPCRBatch, maxNACase=0, maxNAControl=0, hkgs, contrastM, case, control)
+ standardGeneric("gM_ddAvgCt")
+)
+setMethod("gM_ddAvgCt", signature = "qPCRBatch", definition =
+ function(qPCRBatch, maxNACase, maxNAControl, hkgs, contrastM, case, control) {
+ if(length(hkgs)<=1) stop("More than 1 houskeeping gene required")
+ hkgs <- make.names(hkgs)
+ if(FALSE %in% (hkgs %in% featureNames(qPCRBatch))) stop("invalid housekeeping gene")
+
+
+
+ case <- row.names(contrastM)[contrastM[,case] == 1]
+ control <- row.names(contrastM)[contrastM[,control] == 1]
+ expM <- exprs(qPCRBatch)
+ caseM <- expM[,case]
+ controlM <- expM[,control]
+ hkgM <- expM[hkgs,]
+#return(hkgM)
+ if (TRUE %in% apply(apply(hkgM, 1, is.na),2,sum)>0) {
+ warning("NAs present in housekeeping genes. NAs will be excluded when combining housekeepers to make pseudogenes")
+ if (0 %in% apply(! apply(hkgM, 1, is.na),2,sum)) stop("Need at least 1 non NA for each housekeeper")
+ }
+ hkgMCase <- caseM[hkgs, ]
+ hkgMControl <- controlM[hkgs, ]
+
+hkgVCase <- apply(hkgMCase, 2, geomMean, na.rm=TRUE)
+hkgVControl <- apply(hkgMControl, 2, geomMean, na.rm=TRUE)
+ ##############################################
+ # Here we use the geoMean function to find the geometric mean of the chosen housekeepers
+
+
+ if(! FALSE %in% is.na(hkgVCase) || ! FALSE %in% is.na(hkgVControl)) stop("Need at least 1 non NA for the housekeeper")
+ gMhkgControl <- geomMean(hkgVControl, na.rm=TRUE)
+ gMhkgCase <- geomMean(hkgVCase, na.rm=TRUE)
+
+ ddCts <- vector(length=length(featureNames(qPCRBatch)))
+ minddCts <- vector(length=length(ddCts))
+ maxddCts <- vector(length=length(ddCts))
+
+ i <- 1
+ for (detector in featureNames(qPCRBatch)) {
+ VCase <- caseM[detector,]
+ VControl <- controlM[detector,]
+ if(! FALSE %in% is.na(VCase)) {
+ warning("No Detector for Case")
+ dCtCase <- rep(NA, length(VCase))
+ sdCase <- NA
+ }
+# if(length(VCase) == 1) {
+# warning("Only one Detector for Case")
+# dCtCase <- VCase
+# sdCase <- NA
+# }
+ else {
+ dCtCase <- geomMean(VCase, na.rm=T) - gMhkgCase
+ sdCase <- sd(VCase, na.rm=TRUE)
+ }
+# if(length(VControl) == 1) {
+# warning("Only one Detector for Control")
+# dCtControl <- VControl
+# }
+ if(! FALSE %in% is.na(VControl)) {
+ warning("No Detector for Control")
+ dCtControl <- rep(NA, length(VControl))
+ sdControl <- NA
+ }
+ else {
+ dCtControl <- geomMean(VControl, na.rm=TRUE) - gMhkgControl
+ sdControl <- sd(VControl, na.rm=TRUE)
+ }
+ if(sum(is.na(VCase)) > maxNACase) {
+ dCtCase <- NA
+ }
+ if(sum(is.na(VControl)) > maxNAControl) {
+ dCtControl <- NA
+ }
+ ddCt <- (dCtCase - dCtControl)
+
+ if(is.na(ddCt)) {
+ if(is.na(dCtCase) && ! is.na(dCtControl)) ddCt <- "-"
+ else if(is.na(dCtControl) && ! is.na(dCtCase)) ddCt <- "+"
+ else if(is.na(dCtControl) && is.na(dCtCase)) ddCt <- NA
+ minddCts[i] <- NA
+ maxddCts[i] <- NA
+ ddCts[i] <- ddCt
+ }
+ else {
+ if(is.na(sdCase) | is.na(sdControl)) {
+ minddCts[i] <- NA
+ maxddCts[i] <- NA
+ }
+ else {
+ sdCt <- sqrt((sdCase^2) + (sdControl^2))
+ minddCts[i] <- 2 ^ -(ddCt + sdCt)
+ maxddCts[i] <- 2 ^ -(ddCt - sdCt)
+ ddCts[i] <- 2^-ddCt
+ }
+ }
+ i <- i+1
+ }
+ return(cbind(featureNames(qPCRBatch),ddCts,minddCts,maxddCts))
+ }
+)
Modified: pkg/NormqPCR/inst/doc/NormqPCR.Rnw
===================================================================
--- pkg/NormqPCR/inst/doc/NormqPCR.Rnw 2010-04-15 00:09:25 UTC (rev 100)
+++ pkg/NormqPCR/inst/doc/NormqPCR.Rnw 2010-05-13 21:31:12 UTC (rev 101)
@@ -378,7 +378,7 @@
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))
+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
@@ -392,8 +392,7 @@
<< ddCt >>=
hkg <- "Actb-Rn00667869_m1"
-ddCt.taqman <- deltaDeltaCt(qPCRBatch.taqman,1,1,hkg, contM,
-"interestingPhenotype","wildTypePhenotype")
+ddCt.taqman <- deltaDeltaCt(qPCRBatch = qPCRBatch.taqman, maxNACase=1, maxNAControl=1, hkg=hkg, contrastM=contM, case="interestingPhenotype", control="wildTypePhenotype")
head(ddCt.taqman)
@
@@ -408,8 +407,7 @@
<< ddCt Avg >>=
hkg <- "Actb-Rn00667869_m1"
-ddCtAvg.taqman <- deltaDeltaAvgCt(qPCRBatch.taqman,1,1,hkg, contM,
-"interestingPhenotype","wildTypePhenotype")
+ddCtAvg.taqman <- deltaDeltaAvgCt(qPCRBatch = qPCRBatch.taqman, maxNACase=1, maxNAControl=1, hkg=hkg, contrastM=contM, case="interestingPhenotype", control="wildTypePhenotype")
head(ddCtAvg.taqman)
@
%-------------------------------------------------------------------------------
@@ -419,12 +417,36 @@
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.
+geometric mean of these values to form a "pseudo-housekeeper" from which the other
+values are subtracted. For the dataset above:
+<<taqman gM>>=
+path <- system.file("exData", package = "ReadqPCR")
+taqman.example <- paste(path, "/example.txt", sep="")
+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")
+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
+to the \code{deltaDeltaAvgCt} method shown above, intended for use when the samples are being compared from different experiments.
+<<taqman gM Avg>>=
+path <- system.file("exData", package = "ReadqPCR")
+taqman.example <- paste(path, "/example.txt", sep="")
+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")
+ddAvgCt.gM.taqman <- gM_ddAvgCt(qPCRBatch = qPCRBatch.taqman, maxNACase=1, maxNAControl=1, hkgs=hkgs, contrastM=contM, case="interestingPhenotype", control="wildTypePhenotype")
+head(ddAvgCt.gM.taqman)
+@
+
+
%-------------------------------------------------------------------------------
\begin{thebibliography}{1}
More information about the Qpcr-commits
mailing list