DeLimma <- function(mtrx, grps, paired=FALSE, logged=TRUE) {
require(DEGandMore);
require(limma);
prepared <- PrepareDe(mtrx, grps, paired);
mtrx <- prepared[[1]];
grps <- prepared[[2]];
paired <- prepared[[3]];
n <- sapply(grps, length);
if (paired & n[1]==n[2]) {
pair <- factor(rep(1:n[1], 2));
comp <- factor(rep(names(grps), sapply(grps, length)));
design <- model.matrix(~pair+comp);
fit <- lmFit(mtrx, design);
fit <- eBayes(fit);
res <- topTable(fit, coef=colnames(design)[ncol(design)], number=nrow(mtrx));
res <- res[rownames(mtrx), , drop=FALSE];
} else {
design <- cbind(Ctrl=rep(1, ncol(mtrx)), Comp=rep(1:2, sapply(grps, length)));
fit <- lmFit(mtrx, design);
fit <- eBayes(fit);
res <- topTable(fit, coef=2, number=nrow(mtrx));
res <- res[rownames(mtrx), , drop=FALSE];
}
m1 <- rowMeans(mtrx[, grps[[1]], drop=FALSE]);
m2 <- rowMeans(mtrx[, grps[[2]], drop=FALSE]);
pv <- res[, 4];
qv <- p.adjust(pv, method='BH');
if (logged) l2 <- m2-m1 else l2 <- log2(m2/pmax(m1, min(m1[m1>0])/2));
m1[is.na(m1)] <- 0;
m2[is.na(m2)] <- 0;
l2[is.na(l2)] <- 0;
pv[is.na(pv)] <- 1;
qv[is.na(qv)] <- 1;
s <- cbind(m1, m2, m2-m1, l2, pv, qv);
colnames(s) <- c(paste('Mean', names(grps), sep='_'), 'Mean_Change', 'LogFC', 'Pvalue', 'FDR');
rownames(s) <- rownames(mtrx);
list(stat=s, group=grps, limma=fit);
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.