1 | limma_voom_DE(infile, a, b, method)
|
infile |
Input file should be in order of condition(first set will be one condition and next condition will be next set) |
a |
Number of replicates for first set of condition |
b |
Number of replicate for second set of condition |
method |
FDR |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (infile, a, b, method)
{
file = basename(infile)
data = read.table(infile, header = T, row.names = 1)
group = as.factor(c(rep("x", a), rep("y", b)))
head(data)
nf <- calcNormFactors(data, method = "TMM")
design <- model.matrix(~group)
design
y <- voom(data, design, plot = TRUE, lib.size = colSums(data) *
nf)
fit <- lmFit(y, design)
fit <- eBayes(fit)
toptable = topTable(fit, coef = ncol(design), number = 20000,
adjust.method = "BH")
p.adjusted <- p.adjust(fit$p.value[, 2], method = method)
toptable_out = paste0(Sys.Date(), file, "_toptable_rsults.txt")
write.table(toptable, toptable_out, quote = FALSE, row.names = FALSE)
results_limma <- cbind(fit$coeff, fit$p.value[, 2], p.adjusted)
colnames(results_limma) <- c("av_expr", "2LogFC", "pvalue",
"adjusted_pvalue")
results_limma <- results_limma[order(p.adjusted), ]
results_limma_genes = cbind(rownames(results_limma), combat)
dim(results_limma_genes)
colnames(results_limma_genes) <- c("Genes", colnames(results_limma))
output = paste0(Sys.Date(), file, "_", "DE_results.txt")
write.table(results_limma, output, quote = FALSE, row.names = FALSE)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.