# LMGene
# http://www.bioconductor.org/packages/release/bioc/html/LMGene.html
# LMGene {LMGene}
DeLMGene <- function(mtrx, grps, paired=FALSE) {
require(DEGandMore);
require(LMGene);
require(Biobase);
require(tools);
prepared <- PrepareDe(mtrx, grps, paired);
mtrx <- prepared[[1]];
grps <- prepared[[2]];
paired <- prepared[[3]];
n <- sapply(grps, length);
if (paired & n[1]==n[2]) {
prs <- rep(paste('pair', 1:n[1], sep='_'), 2);
cnd <- rep(names(grps), n);
df <- data.frame(condition=cnd, pair=prs);
mt <- data.frame(labelDescription=c('condition', 'pair'));
} else {
cnd <- rep(names(grps), n);
df <- data.frame(condition=cnd);
mt <- data.frame(labelDescription=c('condition'));
}
lvls <- lapply(names(df), function(nm) df[[nm]]);
rownames(df) <- colnames(mtrx);
names(lvls) <- names(df);
anno <- AnnotatedDataFrame(data=df, varMetadata = mt);
eset <- ExpressionSet(mtrx, phenoData = anno);
esti <- tranest(eset, mult = TRUE, method=2);
tran <- transeS(neweS(mtrx, lvls), lambda = esti[[1]], alpha = esti[[2]]);
res <- genediff(tran, method='MOMlog', verbose = FALSE); # use method "MLE" sometimes throw converge error
ex <- exprs(tran);
pv <- res$Gene.Specific[, 1];
m1 <- rowMeans(ex[, grps[[1]], drop=FALSE]);
m2 <- rowMeans(ex[, grps[[2]], drop=FALSE]);
m1 <- exp(m1)/2;
m2 <- exp(m2)/2;
l2 <- m2-m1;
qv <- p.adjust(pv, method='BH');
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, LMGene=res);
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.