[Qpcr-commits] r138 - in pkg/NormqPCR: . R inst/doc man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Jun 30 00:26:15 CEST 2011
Author: jperkins
Date: 2011-06-30 00:26:15 +0200 (Thu, 30 Jun 2011)
New Revision: 138
Modified:
pkg/NormqPCR/DESCRIPTION
pkg/NormqPCR/R/allGenerics.R
pkg/NormqPCR/R/deltaCt.R
pkg/NormqPCR/R/deltaDeltaCt.R
pkg/NormqPCR/R/plotDCt.R
pkg/NormqPCR/R/selectHKs.R
pkg/NormqPCR/inst/doc/NormqPCR.Rnw
pkg/NormqPCR/man/selectHKs.Rd
Log:
made selectHKs S4 & changed examples so it dispatched on qPCR.Batch (in fact it now only dispatches on qPCR.Batch, this is valid since ReadqPCR can read in the data to a qPCRBatch easily)
Modified: pkg/NormqPCR/DESCRIPTION
===================================================================
--- pkg/NormqPCR/DESCRIPTION 2011-06-23 19:59:04 UTC (rev 137)
+++ pkg/NormqPCR/DESCRIPTION 2011-06-29 22:26:15 UTC (rev 138)
@@ -1,5 +1,5 @@
Package: NormqPCR
-Version: 0.99.1
+Version: 0.99.0
Date: 2011-06-23
Title: Functions for normalisation of RT-qPCR data
Description: Functions for normalisation of real-time quantitative PCR data
Modified: pkg/NormqPCR/R/allGenerics.R
===================================================================
--- pkg/NormqPCR/R/allGenerics.R 2011-06-23 19:59:04 UTC (rev 137)
+++ pkg/NormqPCR/R/allGenerics.R 2011-06-29 22:26:15 UTC (rev 138)
@@ -27,3 +27,8 @@
function(qPCRBatch, ...)
standardGeneric("makeAllNAs")
)
+
+setGeneric("selectHKs",
+ function(qPCRBatch, ...)
+ standardGeneric("selectHKs")
+)
Modified: pkg/NormqPCR/R/deltaCt.R
===================================================================
--- pkg/NormqPCR/R/deltaCt.R 2011-06-23 19:59:04 UTC (rev 137)
+++ pkg/NormqPCR/R/deltaCt.R 2011-06-29 22:26:15 UTC (rev 138)
@@ -1,18 +1,15 @@
setMethod("deltaCt", signature = "qPCRBatch", definition =
function(qPCRBatch, hkgs, combineHkgs=FALSE, calc="arith") {
hkgs <- make.names(hkgs)
-# cat("hkgs",hkgs)
if(FALSE %in% (hkgs %in% featureNames(qPCRBatch))) stop ("given housekeeping gene, ", hkgs," not found in file. Ensure entered housekeeping genes appear in the file")
expM <- exprs(qPCRBatch)
hkgM <- expM[hkgs, ]
-# print(hkgM)
if(length(hkgs) > 1) {
if (TRUE %in% apply(hkgM, 1, is.na)) {
warning("NAs present in housekeeping genes readings")
if (0 %in% apply(! apply(hkgM, 1, is.na),2,sum)) stop("Need at least 1 non NA for each housekeeper")
}
-# hkgV <- apply(hkgM,2,geomMean,na.rm=TRUE) CANT CURRENTLY DO THIS BECAUS EOF PROBLEMS WITH DOUBLE NAS
hkgV <- vector(length = dim(hkgM)[2])
if(calc == "arith") {
for(i in 1:dim(hkgM)[2]) {
@@ -26,15 +23,13 @@
hkgV[i] <- geomMean(hkgM[,i], na.rm=TRUE)
}
}
-# cat("hkgV",hkgV,"\n")
}
else {
if(TRUE %in% is.na(hkgM)) {
- warning("NAs present in housekeeping gene readings")
- if(! FALSE %in% is.na(hkgM)) stop("Need at least 1 non NA for the housekeeper")
+ warning("NAs present in housekeeping gene readings")
+ if(! FALSE %in% is.na(hkgM)) stop("Need at least 1 non NA for the housekeeper")
}
- hkgV <- hkgM # Because it's a vector really anyway
-# cat(hkgV)
+ hkgV <- hkgM # Because it's a vector really anyway: we only have one NA
}
exprs(qPCRBatch) <- t(t(exprs(qPCRBatch)) - hkgV)
return(qPCRBatch)
Modified: pkg/NormqPCR/R/deltaDeltaCt.R
===================================================================
--- pkg/NormqPCR/R/deltaDeltaCt.R 2011-06-23 19:59:04 UTC (rev 137)
+++ pkg/NormqPCR/R/deltaDeltaCt.R 2011-06-29 22:26:15 UTC (rev 138)
@@ -1,10 +1,8 @@
setMethod("deltaDeltaCt", signature = "qPCRBatch", definition =
function(qPCRBatch, maxNACase=0, maxNAControl=0, hkgs, contrastM, case, control, paired=TRUE, hkgCalc="arith", statCalc="arith") {
hkgs <- make.names(hkgs)
+ if(length(hkgs) < 1) stop("Not enough hkgs given")
-# if(combineHkgs == TRUE) {
- if(length(hkgs) < 1) stop("Not enough hkgs given")
-# }
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))
@@ -18,8 +16,6 @@
if(length(hkgs) > 1) {
hkgMCase <- caseM[hkgs, ]
hkgMControl <- controlM[hkgs, ]
- #hkgVCase <- apply(hkgMCase, 2, geomMean, na.rm=TRUE)
- #hkgVControl <- apply(hkgMControl, 2, geomMean, na.rm=TRUE)
if(hkgCalc == "arith") {
hkgVCase <- apply(hkgMCase, 2, mean, na.rm=TRUE)
hkgVControl <- apply(hkgMControl, 2, mean, na.rm=TRUE)
@@ -49,7 +45,6 @@
for (detector in featureNames(qPCRBatch)) {
VCase <- caseM[detector,]
VControl <- controlM[detector,]
-#stop("length of case is:",VCase,"_",length(VCase))
if(length(VCase) == 1) {
warning("Only one Detector for Control")
dCtCase <- VCase
@@ -57,7 +52,6 @@
} else if(! FALSE %in% is.na(VCase)) {
warning("No Detector for Case")
dCtCase <- rep(NA, length = VCase)
-# dCtControl <- NA
sdCase <- NA
} else {
if(statCalc == "geom") {
@@ -81,7 +75,6 @@
} else if(! FALSE %in% is.na(VControl)) {
warning("No Detector for Control")
dCtControl <- rep(NA, length = VControl)
-# dCtCase <- NA
sdControl <- NA
} else {
if(statCalc == "geom") {
@@ -137,7 +130,6 @@
}
}
if(statCalc == "geom") {
-# ddCts[i] <- (dCtCase/dCtControl)
dCtCases[i] <- dCtCase
sdCtCases[i] <- sdCase
dCtControls[i] <- dCtControl
@@ -156,6 +148,5 @@
ddCtTable <- as.data.frame(cbind(featureNames(qPCRBatch),format(dCtCases, digits=4),format(sdCtCases, digits=4),format(dCtControls, digits=4),format(sdCtControls, digits=4),format(ddCts,digits=4),format(minddCts, digits=4),format(maxddCts, digits=4)))
names(ddCtTable) <- c("ID", paste("2^-dCt",case,sep="."), paste(case,"sd",sep="."), paste("2^-dCt",control,sep="."), paste(control,"sd",sep="."),"2^-ddCt","2^-ddCt.min", "2^-ddCt.max")
return(ddCtTable)
-# stop(ddCtTable)
}
)
Modified: pkg/NormqPCR/R/plotDCt.R
===================================================================
--- pkg/NormqPCR/R/plotDCt.R 2011-06-23 19:59:04 UTC (rev 137)
+++ pkg/NormqPCR/R/plotDCt.R 2011-06-29 22:26:15 UTC (rev 138)
@@ -17,10 +17,5 @@
plotU <- plotCts + plotSds
plotL <- plotCts - plotSds
}
-# print(plotCts)
-
-# cat(plotU,"\t")
-# cat(plotL,"\n")
-
barplot2(..., height=t(plotCts), plot.ci=TRUE, ci.u=t(plotU), ci.l=t(plotL), beside=T)
}
Modified: pkg/NormqPCR/R/selectHKs.R
===================================================================
--- pkg/NormqPCR/R/selectHKs.R 2011-06-23 19:59:04 UTC (rev 137)
+++ pkg/NormqPCR/R/selectHKs.R 2011-06-29 22:26:15 UTC (rev 138)
@@ -8,21 +8,11 @@
## Symbols: character, symbols for variables/columns
## trace: locical, print information
## na.rm: remove NA values
-selectHKs <- function(x, group, method = "geNorm", minNrHKs = 2,
- log = TRUE, Symbols, trace = TRUE, na.rm = TRUE){
- if(method == "geNorm") {
- if(class(x) == "qPCRBatch") x <- t(exprs(x))
- }
- if(method == "NormFinder") {
- if(class(x) == "qPCRBatch"){
- if (missing(group)) group <- pData(x)[,"Group"]
- x <- t(exprs(x))
- }
- else stop("'x' must be of class qPCRBatch")
- }
- if(!is.matrix(x) & !is.data.frame(x))
- stop("'x' needs to be of class matrix or data.frame")
- if(is.data.frame(x)) x <- data.matrix(x)
+setMethod("selectHKs", signature = "qPCRBatch", definition =
+ function(qPCRBatch, group, method = "geNorm", minNrHKs = 2,
+ log = TRUE, Symbols, trace = TRUE, na.rm = TRUE){
+ x <- t(exprs(qPCRBatch))
+ x <- data.matrix(x)
n <- ncol(x)
if(n < 3)
stop("you need data from at least 3 variables/columns")
@@ -67,14 +57,14 @@
}
if(trace){
- cat("###############################################################\n")
- cat("Step ", n-i+1, ":\n")
- cat("stability values M:\n")
+ message("###############################################################\n")
+ message("Step ", n-i+1, ":\n")
+ message("stability values M:\n")
print(sort(M))
- cat("average stability M:\t", meanM[n-i+1], "\n")
+ message("average stability M:\t", meanM[n-i+1], "\n")
if(i > 2){
- cat("variable with lowest stability (largest M value):\t", Symbols[ind], "\n")
- cat("Pairwise variation, (", i-1, "/", i, "):\t", V[n-i+1], "\n")
+ message("variable with lowest stability (largest M value):\t", Symbols[ind], "\n")
+ message("Pairwise variation, (", i-1, "/", i, "):\t", V[n-i+1], "\n")
}
}
x <- x[,-ind]
@@ -95,11 +85,11 @@
rho.min[1] <- rho[b[1]]
if(trace){
- cat("###############################################################\n")
- cat("Step ", 1, ":\n")
- cat("stability values rho:\n")
+ message("###############################################################\n")
+ message("Step ", 1, ":\n")
+ message("stability values rho:\n")
print(sort(rho))
- cat("variable with highest stability (smallest rho value):\t", Symbols[b[1]], "\n")
+ message("variable with highest stability (smallest rho value):\t", Symbols[b[1]], "\n")
}
for(i in 2:minNrHKs){
rho[b[i-1]] <- NA
@@ -116,11 +106,11 @@
R[i] <- Symbols[b[i]]
rho.min[i] <- rho[b[i]]
if(trace){
- cat("###############################################################\n")
- cat("Step ", i, ":\n")
- cat("stability values rho:\n")
+ message("###############################################################\n")
+ message("Step ", i, ":\n")
+ message("stability values rho:\n")
print(sort(rho[!is.na(rho)]))
- cat("variable with highest stability (smallest rho value):\t", Symbols[b[i]], "\n")
+ message("variable with highest stability (smallest rho value):\t", Symbols[b[i]], "\n")
}
}
names(R) <- 1:minNrHKs
@@ -129,4 +119,5 @@
}else{
stop("specified method not yet implemented")
}
-}
+ }
+)
Modified: pkg/NormqPCR/inst/doc/NormqPCR.Rnw
===================================================================
--- pkg/NormqPCR/inst/doc/NormqPCR.Rnw 2011-06-23 19:59:04 UTC (rev 137)
+++ pkg/NormqPCR/inst/doc/NormqPCR.Rnw 2011-06-29 22:26:15 UTC (rev 138)
@@ -216,23 +216,21 @@
\code{minNrHK} genes remain.
<<geNorm>>=
-geNorm <- t(exprs(geNorm.qPCRBatch))
-
tissue <- as.factor(c(rep("BM", 9), rep("FIB", 20), rep("LEU", 13),
rep("NB", 34), rep("POOL", 9)))
-res.BM <- selectHKs(geNorm[tissue == "BM",], method = "geNorm",
- Symbols = colnames(geNorm), minNrHK = 2, log = FALSE)
-res.POOL <- selectHKs(geNorm[tissue == "POOL",], method = "geNorm",
- Symbols = colnames(geNorm), minNrHK = 2,
+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[tissue == "FIB",], method = "geNorm",
- Symbols = colnames(geNorm), minNrHK = 2,
+res.FIB <- selectHKs(geNorm.qPCRBatch[,tissue == "FIB"], method = "geNorm",
+ Symbols = featureNames(geNorm.qPCRBatch), minNrHK = 2,
trace = FALSE, log = FALSE)
-res.LEU <- selectHKs(geNorm[tissue == "LEU",], method = "geNorm",
- Symbols = colnames(geNorm), minNrHK = 2,
+res.LEU <- selectHKs(geNorm.qPCRBatch[,tissue == "LEU"], method = "geNorm",
+ Symbols = featureNames(geNorm.qPCRBatch), minNrHK = 2,
trace = FALSE, log = FALSE)
-res.NB <- selectHKs(geNorm[tissue == "NB",], method = "geNorm",
- Symbols = colnames(geNorm), minNrHK = 2,
+res.NB <- selectHKs(geNorm.qPCRBatch[,tissue == "NB"], method = "geNorm",
+ Symbols = featureNames(geNorm.qPCRBatch), minNrHK = 2,
trace = FALSE, log = FALSE)
@
We obtain the following ranking of genes (cf. Table~3 in Vand02)
@@ -266,7 +264,7 @@
legend("topright", legend = c("BM", "POOL", "FIB", "LEU", "NB"),
fill = mypalette)
@
-\myinfig{1}{NormqPCR-fig2.pdf}QCqPCR/inst/doc/QCqPCR.Rnw
+\myinfig{1}{NormqPCR-fig2.pdf}
\par
Second, we plot the pairwise variation for each cell type (cf. Figure~3~(a) in
Vand02)
@@ -333,11 +331,13 @@
procedure
is available via function \code{selectHKs}.
<<NormFinder2>>=
+group <- pData(Colon.qPCRBatch)[,"Group"]
selectHKs(Colon.qPCRBatch,
- log = FALSE, trace = TRUE,
+ log = FALSE, trace = TRUE, group = group,
Symbols = featureNames(Colon.qPCRBatch), minNrHKs = 12,
method = "NormFinder")$ranking
-selectHKs(Bladder.qPCRBatch,
+group <- pData(Bladder.qPCRBatch)[,"Group"]
+selectHKs(Bladder.qPCRBatch, group = group,
log = FALSE, trace = FALSE,
Symbols = featureNames(Bladder.qPCRBatch), minNrHKs = 13,
method = "NormFinder")$ranking
Modified: pkg/NormqPCR/man/selectHKs.Rd
===================================================================
--- pkg/NormqPCR/man/selectHKs.Rd 2011-06-23 19:59:04 UTC (rev 137)
+++ pkg/NormqPCR/man/selectHKs.Rd 2011-06-29 22:26:15 UTC (rev 138)
@@ -1,16 +1,21 @@
\name{selectHKs}
\alias{selectHKs}
+\alias{selectHKs-methods}
+\alias{selectHKs,qPCRBatch-method}
\title{ Selection of reference/housekeeping genes }
\description{
This function can be used to determine a set of reference/housekeeping (HK)
genes for gene expression experiments.
}
\usage{
-selectHKs(x, group, method = "geNorm", minNrHKs = 2, log = TRUE, Symbols,
+selectHKs(qPCRBatch, \dots)
+
+\S4method{selectHKs}{qPCRBatch}(qPCRBatch, group, method = "geNorm", minNrHKs = 2, log = TRUE, Symbols,
trace = TRUE, na.rm = TRUE)
}
\arguments{
- \item{x}{ matrix or data.frame containing the data }
+ \item{qPCRBatch}{ qPCRBatch, containing the data }
+ \item{\dots}{ Extra arguments, detailed below }
\item{group}{ optional factor not used by all methods, hence may be missing }
\item{method}{ method to compute most stable genes }
\item{minNrHKs}{ minimum number of HK genes that should be considered }
@@ -67,8 +72,12 @@
%}
%\seealso{}
\examples{
-data(geNorm)
-res.BM <- selectHKs(geNorm[1:9,], method = "geNorm", Symbols = names(geNorm),
- minNrHK = 2, log = FALSE, trace = TRUE, na.rm = TRUE)
+path <- system.file("exData", package = "NormqPCR")
+geNorm.example <- file.path(path, "qPCR.geNorm.txt")
+geNorm.qPCRBatch <- read.qPCR(geNorm.example)
+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)
}
\keyword{data}
More information about the Qpcr-commits
mailing list