Nothing
# last modified: 3 June 2008
# This reproduces Table 2 of Wichert et al. (2004)
# Note that for time series with even sample size the numbers
# below differ from the ones published in Wichert et al (2004)
# because in the original publication the g-statistic was erroneously
# computed with the intensity at frequency \pi included.
# This is fixed in version 1.0.4 of GeneCycle.
# Note: in the following multiple testing is performed in four ways:
# - as in Wichert et al. using the original Benjamini-Hochberg (1995)
# algorithm (this yields the "conservative" estimate)
# - using the FDR approach of Storey and Tibshirani (2003)
# (this yields the "less conservative" estimate)
# - using semiparametric Fdr and local fdr, as estimated by the
# algorithm in fdrtool
#
# Other FDR controlling approaches may lead to different sets of
# significant genes.
library("GeneCycle")
library("qvalue")
# load data sets
data(caulobacter)
# download these files from http://www.strimmerlab.org/data.html
load("spellmann-yeast.rda.gz")
load("fibroblasts.rda.gz")
load("humanhela.rda.gz")
# calculate p-values and determine number of significant genes
find.periodic.genes <- function(dataset)
{
cat("Computing p-values ...\n")
pval = fisher.g.test(dataset)
n1 = sum(qvalue(pval, lambda=0)$qvalues < 0.05)
n2 = sum(qvalue(pval)$qvalues < 0.05)
fdr.out <- fdrtool(pval, statistic="pvalue", plot=FALSE)
n3 <- sum(fdr.out$qval < 0.05)
n4 <- sum(fdr.out$lfdr < 0.2)
cat("Conservative estimate (Wichert et al.) =", n1, "\n")
cat("Less conservative estimate =", n2, "\n")
cat("Semiparametric Fdr < 0.05 (fdrtool) =", n3, "\n")
cat("Semiparametric fdr < 0.2 (fdrtool) =", n4, "\n")
}
### data analysis ###
dim(cdc15) # 24 x 4289
find.periodic.genes(cdc15)
# Conservative estimate (Wichert et al.) = 221
# Less conservative estimate = 294
# Semiparametric Fdr < 0.05 (fdrtool) = 275
# Semiparametric fdr < 0.2 (fdrtool) = 319
dim(cdc28) # 17 x 1365
find.periodic.genes(cdc28)
# Conservative estimate (Wichert et al.) = 27 (!! note: typographic error in ms. !!)
# Less conservative estimate = 56
# Semiparametric Fdr < 0.05 (fdrtool) = 63
# Semiparametric fdr < 0.2 (fdrtool) = 179
dim(alpha) # 18 x 4415
find.periodic.genes(alpha)
# Conservative estimate (Wichert et al.) = 170
# Less conservative estimate = 203
# Semiparametric Fdr < 0.05 (fdrtool) = 213
# Semiparametric fdr < 0.2 (fdrtool) = 284
dim(elution) # 14 x 5695
find.periodic.genes(elution)
# Conservative estimate (Wichert et al.) = 72
# Less conservative estimate = 101
# Semiparametric Fdr < 0.05 (fdrtool) = 125
# Semiparametric fdr < 0.2 (fdrtool) = 490
dim(caulobacter) # 11 x 1444
find.periodic.genes(caulobacter)
# Conservative estimate (Wichert et al.) = 45
# Less conservative estimate = 45
# Semiparametric Fdr < 0.05 (fdrtool) = 52
# Semiparametric fdr < 0.2 (fdrtool) = 94
dim(humanN2) # 13 x 4574
find.periodic.genes(humanN2)
# Conservative estimate (Wichert et al.) = 0
# Less conservative estimate = 0
# Semiparametric Fdr < 0.05 (fdrtool) = 0
# Semiparametric fdr < 0.2 (fdrtool) = 0
dim(humanN3) # 12 x 5079
find.periodic.genes(humanN3)
# Conservative estimate (Wichert et al.) = 0
# Less conservative estimate = 0
# Semiparametric Fdr < 0.05 (fdrtool) = 0
# Semiparametric fdr < 0.2 (fdrtool) = 0
dim(score1) # 12 x 14728
find.periodic.genes(score1)
# Conservative estimate (Wichert et al.) = 2
# Less conservative estimate = 2
# Semiparametric Fdr < 0.05 (fdrtool) = 2
# Semiparametric fdr < 0.2 (fdrtool) = 1
dim(score2) # 26 x 15472
find.periodic.genes(score2)
# Conservative estimate (Wichert et al.) = 72
# Less conservative estimate = 211
# Semiparametric Fdr < 0.05 (fdrtool) = 185
# Semiparametric fdr < 0.2 (fdrtool) = 430
dim(score3) # 48 x 39724
find.periodic.genes(score3)
# Conservative estimate (Wichert et al.) = 4250
# Less conservative estimate = 4505
# Semiparametric Fdr < 0.05 (fdrtool) = 4549
# Semiparametric fdr < 0.2 (fdrtool) = 4447
dim(score4) # 19 x 39192
find.periodic.genes(score4)
# Conservative estimate (Wichert et al.) = 57
# Less conservative estimate = 60
# Semiparametric Fdr < 0.05 (fdrtool) = 60
# Semiparametric fdr < 0.2 (fdrtool) = 122
dim(score5) # 9 x 34890
find.periodic.genes(score5)
# Conservative estimate (Wichert et al.) = 0
# Less conservative estimate = 0
# Semiparametric Fdr < 0.05 (fdrtool) = 0
# Semiparametric fdr < 0.2 (fdrtool) = 0
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.