Expression data analysis use limma

source(“http://bioconductor.org/biocLite.R”)
options(BioC_mirror=”http://mirrors.ustc.edu.cn/bioc/”)
biocLite(c(“Biobase”, “affyPLM”, “limma”))

library(“Biobase”)
library(“affyPLM”)
library(limma)

exprsFile <- "counts.txt" raw_exprs <- as.matrix(read.table(exprsFile, header=TRUE, sep="\t", row.names=1, as.is=TRUE)) class(raw_exprs) dim(raw_exprs) colnames(raw_exprs) head(raw_exprs[,1:5]) boxplot(raw_exprs) log2exp <- log2(raw_exprs) boxplot(log2exp) #raw_eset <- new("ExpressionSet", exprs=raw_exprs) log2_eset <- new("ExpressionSet", exprs=log2exp) # qt_exps <- normalize.ExpressionSet.quantiles(raw_eset) # qt_exps <- normalize.ExpressionSet.scaling(raw_eset, trim=0.02, baseline=-1) qt_exps <- normalize.ExpressionSet.quantiles(log2_eset) # qt_exps <- normalize.ExpressionSet.scaling(log2_eset, trim=0.02, baseline=-1) boxplot(qt_exps) exp <- exprs(qt_exps) write.csv(exp, file = "exp_quantile_scaling2.csv") # write.table(exp, file = "exp_quantile_scaling2.txt", sep="\t") # anti_log2exp <- 2^exp # write.csv(anti_log2exp, file = "anti_log2exp.csv") # run limma # D2 vs D0 eset <- qt_exps[,c(1,2,3,4,5,6,7,8)] colnames(eset) strain <- c("D0","D0","D0","D0","D2","D2","D2","D2") design <- model.matrix(~factor(strain)) colnames(design) <- c("D0","D2VD0") design fit <- lmFit(eset, design) fit2 <- eBayes(fit) options(digits=2) result <- topTable(fit2, coef=2, n=38535, adjust="BH") write.csv(result, file = "limma.stat.fit-d2vd0.csv") # write.table(result, file = "limma.stat.fit-d2vd0.txt", sep="\t") # D7 vs D0 eset <- qt_exps[,c(1,2,3,4,9,10,11,12)] colnames(eset) strain <- c("D0","D0","D0","D0","D7","D7","D7","D7") design <- model.matrix(~factor(strain)) colnames(design) <- c("D0","D7VD0") design fit <- lmFit(eset, design) fit2 <- eBayes(fit) options(digits=2) result <- topTable(fit2, coef=2, n=38535, adjust="BH") write.csv(result, file = "limma.stat.fit-d7vd0.csv") # write.table(result, file = "limma.stat.fit-d7vd0.txt", sep="\t") # D7 vs D2 eset <- qt_exps[,c(5,6,7,8,9,10,11,12)] colnames(eset) strain <- c("D2","D2","D2","D2","D7","D7","D7","D7") design <- model.matrix(~factor(strain)) colnames(design) <- c("D2","D7VD2") design fit <- lmFit(eset, design) fit2 <- eBayes(fit) options(digits=2) result <- topTable(fit2, coef=2, n=38535, adjust="BH") write.csv(result, file = "limma.stat.fit-d7vd2.csv") # write.table(result, file = "limma.stat.fit-d7vd2.txt", sep="\t")

Leave a Comment