There are two steps- first you create a matrix that is 4x14- the 4 replicates by the 14 (in your case) genes. Also create an experiment design: a data frame with one column, c(1, 1, 0, 0)
, indicating the first two columns of the matrix are test and the second two are control).
library(limma)
m = as.matrix(cbind(test[, 2:3], control[, 2:3]))
rownames(m) = test[, 1]
d = data.frame(condition=c(1, 1, 0, 0))
Then you use lmFit
to fit the linear model, and eBayes
to adjust it and calculate p-values:
fit = lmFit(m, d)
e = eBayes(fit)
The e
object, which is of the class MArrayLM
, contains the p-values and estimated coefficients, along with some other information:
print(e$coefficients)
condition
# Gene_1 3643.03848
# Gene_2 772.31186
# ...
print(e$p.value)
# condition
# Gene_1 6.299507e-07
# Gene_2 1.333970e-04
names(e)
# [1] "coefficients" "rank" "assign" "qr"
# [5] "df.residual" "sigma" "cov.coefficients" "stdev.unscaled"
# [9] "pivot" "genes" "Amean" "method"
# [13] "design" "df.prior" "s2.prior" "var.prior"
# [17] "proportion" "s2.post" "t" "df.total"
# [21] "p.value" "lods" "F" "F.p.value"
Incidentally, while limma
is better for multiple reasons (it calculates other parameters, provides more options, and performs useful corrections), if you did want to perform a t.test
to get p-values for each gene in a matrix, you could do this:
pvals = apply(m, 1, function(r) t.test(r ~ d$condition)$p.value)