yifanzhang — Apr 13, 2014, 9:18 PM
library(edgeR)
Loading required package: limma
library(limma)
library(splines)
library(DESeq)
Loading required package: BiocGenerics
Loading required package: parallel
Attaching package: 'BiocGenerics'
The following objects are masked from 'package:parallel':
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, parApply, parCapply, parLapply,
parLapplyLB, parRapply, parSapply, parSapplyLB
The following object is masked from 'package:limma':
plotMA
The following object is masked from 'package:stats':
xtabs
The following objects are masked from 'package:base':
anyDuplicated, append, as.data.frame, as.vector, cbind,
colnames, duplicated, eval, evalq, Filter, Find, get,
intersect, is.unsorted, lapply, Map, mapply, match, mget,
order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
rbind, Reduce, rep.int, rownames, sapply, setdiff, sort,
table, tapply, union, unique, unlist
Loading required package: Biobase
Welcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
Loading required package: locfit
locfit 1.5-9.1 2013-03-22
Welcome to 'DESeq'. For improved performance, usability and
functionality, please consider migrating to 'DESeq2'.
library(VennDiagram)
Loading required package: grid
dat <- read.table("bottomly_count_table.tsv", header = TRUE,
row.names = 1)
des <- read.table("bottomly_phenodata.tsv", header = TRUE,
row.names = 1)
# count equal tot zero across all samples
dat1<- dat[rowSums(dat != 0) > 0,]
dat2<-dat
# remove probes count equal to zero in at least one sample in each genotype group
for (i in 1:ncol(dat2)){
dat2<-subset(dat2, dat2[,i]!=0)
}
group <- factor(c(rep("1", 10), rep("2", 11)))
dge.glm1 <- DGEList(counts = dat1, group = group)
design <- model.matrix(~group)
dge.glm.com.disp1 <- estimateGLMCommonDisp(dge.glm1, design, verbose = TRUE)
Disp = 0.03893 , BCV = 0.1973
dge.glm.trend.disp1 <- estimateGLMTrendedDisp(dge.glm.com.disp1, design)
dge.glm.tag.disp1 <- estimateGLMTagwiseDisp(dge.glm.trend.disp1, design)
plotBCV(dge.glm.tag.disp1)
fit1 <- glmFit(dge.glm.tag.disp1, design)
lrt1 <- glmLRT(fit1, coef = 2)
tt.glm1 <- topTags(lrt1, n = Inf)
de.glm1 <- decideTestsDGE(lrt1, p = 0.05, adjust = "BH")
tags.glm1 <- rownames(dge.glm.tag.disp1)[as.logical(de.glm1)]
plotSmear(lrt1, de.tags = tags.glm1)
abline(h = c(-2, 2), col = "blue")
dge.glm2 <- DGEList(counts = dat2, group = group)
dge.glm.com.disp2 <- estimateGLMCommonDisp(dge.glm2, design, verbose = TRUE)
Disp = 0.03704 , BCV = 0.1925
dge.glm.trend.disp2 <- estimateGLMTrendedDisp(dge.glm.com.disp2, design)
dge.glm.tag.disp2 <- estimateGLMTagwiseDisp(dge.glm.trend.disp2, design)
plotBCV(dge.glm.tag.disp2)
fit2 <- glmFit(dge.glm.tag.disp2, design)
lrt2 <- glmLRT(fit2, coef = 2)
tt.glm2 <- topTags(lrt2, n = Inf)
de.glm2 <- decideTestsDGE(lrt2, p = 0.05, adjust = "BH")
tags.glm2 <- rownames(dge.glm.tag.disp2)[as.logical(de.glm2)]
plotSmear(lrt2, de.tags = tags.glm2)
abline(h = c(-2, 2), col = "blue")
deSeqDat <- newCountDataSet(dat, group)
deSeqDat <- estimateSizeFactors(deSeqDat)
deSeqDat <- estimateDispersions(deSeqDat)
results <- nbinomTest(deSeqDat, levels(group)[1], levels(group)[2])
cutoff <- 0.05
DESeq.genes <- rownames(subset(results, padj < cutoff))
norm.factor <- calcNormFactors(dat)
dat.voomed <- voom(dat, design, plot = TRUE, lib.size = colSums(dat) * norm.factor)
fit <- lmFit(dat.voomed, design)
fit <- eBayes(fit)
voom.limma.genes <- rownames(topTable(fit, p.value=cutoff))
dge.glm <- DGEList(counts = dat, group = group)
dge.glm.com.disp <- estimateGLMCommonDisp(dge.glm, design, verbose = TRUE)
Disp = 0.03893 , BCV = 0.1973
dge.glm.trend.disp <- estimateGLMTrendedDisp(dge.glm.com.disp, design)
dge.glm.tag.disp <- estimateGLMTagwiseDisp(dge.glm.trend.disp, design)
fit <- glmFit(dge.glm.tag.disp, design)
lrt <- glmLRT(fit, coef = 2)
tt.glm <- topTags(lrt, n = Inf)
edgeR.genes <- rownames(tt.glm$table[tt.glm$table$FDR < cutoff,])
de.genes <- list(edgeR = edgeR.genes, DESeq = DESeq.genes, voom.limma = voom.limma.genes)
plot.new()
venn.plot <- venn.diagram(de.genes, filename = NULL, fill = c("red", "blue", "green"))
grid.draw(venn.plot)
