Difference between revisions of "GGB-correlation-snippet.R"

Latest revision as of 00:37, 6 June 2007

Code snipits and programs written in R, S or S-PLUS

m <- p <- 1000 nreps <- 5

n <- nreps * 2 pi0 <- 0.95

m1 <- round(m * (1-pi0)) m0 <- p - m1

params <- c(2.74886, 1.36546, 4.12844) # "IPTG-a parameters" k <- 1000

a.shape <- params[1] a0.shape <- params[2] scale <- params[3]

kDE <- rgamma(2*m1, shape=k, rate=k) # p49 2001 NEWTON paper DEscales <- rgamma(2*m1, shape=a0.shape, rate=scale) DEscales <- DEscales * kDE

kEE <- rgamma(m0, shape=k, rate=k) EEscales <- rgamma(m0 , shape=a0.shape, rate=scale) EEscales <- EEscales / kEE

scales <- c(rep(DEscales, each=nreps), rep(EEscales, each=2*nreps)) X <- rgamma(n* p, a.shape, rate=scales) dim(X) <- c(n,p)

X <- t(X)

hist(kDE, prob=TRUE) lines(density(kDE))

hist(kEE, prob=TRUE) lines(density(kEE))

Xbar <- rowMeans(log2(X)) hist(Xbar, breaks=30, prob=TRUE) lines(density(Xbar)) pairs(log2(X))

@@ Line 1: / Line 1: @@
-p <- 1000
+# {{R}}
-n <- 10
+# p (genes) by n (slides) matrix
-a.shape  <- parameters[1]
+# m=p keeping fdr notation
-a0.shape <- parameters[2]
-scale    <-   parameters[3]
-k <- 20
-kcorrelation <- rgamma(2*m1 shape=k, rate=k) # p49 2001 NEWTON paper
+m <- p  <- 1000
-DEscales <- rep(rgamma(m1, shape=a0.shape, rate=scale), each=2)
+nreps <- 5
-DEscales <- DEscales * kcorrelation
-kcorrelation <- rgamma(m0, shape=k, rate=k)
+n   <- nreps * 2
+pi0 <- 0.95
+m1 <- round(m * (1-pi0))
+m0 <- p - m1
+params <- c(2.74886, 1.36546, 4.12844) # "IPTG-a parameters"
+k <- 1000
+a.shape  <- params[1]
+a0.shape <- params[2]
+scale    <- params[3]
+kDE <- rgamma(2*m1, shape=k, rate=k) # p49 2001 NEWTON paper
+DEscales <- rgamma(2*m1, shape=a0.shape, rate=scale)
+DEscales <- DEscales * kDE
+kEE <- rgamma(m0, shape=k, rate=k)
 EEscales <- rgamma(m0 , shape=a0.shape, rate=scale)
-EEscales <- EEscales / kcorrelation
+EEscales <- EEscales / kEE
-scales   <- c(rep(DEscales, each=slides), rep(EEscales, each=2*slides))
+scales   <- c(rep(DEscales, each=nreps), rep(EEscales, each=2*nreps))
+X <- rgamma(n* p, a.shape, rate=scales)
+dim(X) <- c(n,p)
+X <- t(X)
+# Graphical check
+hist(kDE, prob=TRUE)
+lines(density(kDE))
+hist(kEE, prob=TRUE)
+lines(density(kEE))
-X <- rgamma(ngenes*p, a.shape, rate=scales)
+Xbar <- rowMeans(log2(X))
-dim(X) <- c(p, n)
+hist(Xbar, breaks=30, prob=TRUE)
-X <- t(X)  # Data already on raw scale
+lines(density(Xbar))
+pairs(log2(X))