Difference between revisions of "GGBsnippet.R"
From Organic Design wiki
m (# {{R}}) |
|||
(6 intermediate revisions by the same user not shown) | |||
Line 1: | Line 1: | ||
− | p <- 1000 | + | # {{R}} |
− | n <- | + | # p (genes) by n (slides) matrix |
− | a.shape <- | + | # m=p keeping fdr notation |
− | a0.shape <- | + | |
− | scale <- | + | 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" | ||
+ | |||
+ | a.shape <- params[1] | ||
+ | a0.shape <- params[2] | ||
+ | scale <- params[3] | ||
DEscales <- rgamma(2*m1, shape=a0.shape, rate=scale) | DEscales <- rgamma(2*m1, shape=a0.shape, rate=scale) | ||
EEscales <- rgamma(m0 , shape=a0.shape, rate=scale) | EEscales <- rgamma(m0 , shape=a0.shape, rate=scale) | ||
− | scales <- c(rep(DEscales, each= | + | scales <- c(rep(DEscales, each=nreps), rep(EEscales, each=2*nreps)) |
− | |||
X <- rgamma(n* p, a.shape, rate=scales) | X <- rgamma(n* p, a.shape, rate=scales) | ||
− | dim(X) <- c(p | + | dim(X) <- c(n,p) |
+ | |||
X <- t(X) | X <- t(X) | ||
+ | |||
+ | # Graphical check | ||
+ | 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)) |
Latest revision as of 00:38, 6 June 2007
Code snipits and programs written in R, S or S-PLUS
- p (genes) by n (slides) matrix
- m=p keeping fdr notation
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"
a.shape <- params[1] a0.shape <- params[2] scale <- params[3]
DEscales <- rgamma(2*m1, shape=a0.shape, rate=scale) EEscales <- rgamma(m0 , shape=a0.shape, rate=scale)
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))
Xbar <- rowMeans(log2(X)) hist(Xbar, breaks=30, prob=TRUE) lines(density(Xbar)) pairs(log2(X))