A csaw/edgeR pipeline – Figure 2
We use edgeR based testing instead of DESeq2 for the testing. For this we extract the raw counts from the data and then run the tesing pipeline just like csaw does it, using the quasi–likelihood based testing pipeline implemented in edgeR.
Not that it is not possible to define a log–fold–change cutoff for testing in edgeR, so we use a post–hoc cutoff to determine significance.
We then produce a plot similar to Figure 1.
# extract counts
ccc <- counts(DESeq2Data(importedData))
condition <- colData(DESeq2Data(importedData))$condition
y <- DGEList(ccc,
group = condition)
design <- model.matrix(~ 0 + condition)
colnames(design) <- levels(condition)
y <- estimateDisp(y, design = design)
plotBCV(y)
fit <- glmQLFit(y, design, robust=TRUE)
plotQLDisp(fit)
contrast <- makeContrasts(WT-mutant, levels=design)
res <- glmQLFTest(fit, contrast=contrast)
resEdgeR <- topTags(sort.by = "none", res, n=Inf)$table
FDRedgeR <- fdrtool(resEdgeR$PValue, statistic = "pvalue", plot = FALSE)
Step 1... determine cutoff point
Step 2... estimate parameters of null distribution and eta0
Step 3... compute p-values and estimate empirical PDF/CDF
Step 4... compute q-values and local fdr
resEdgeR$ldfr <- FDRedgeR$lfdr
table(resEdgeR$ldfr < 0.2)
FALSE TRUE
1157 1344
lfdrEdgeR <- FDRedgeR$lfdr[1000:2501]
sum(lfdrEdgeR < 0.2)
[1] 1127
countsLog2 <- as.data.frame(log2(counts(DESeq2Data(importedData), norm = TRUE)))
sampleTable <- colData(DESeq2Data(importedData))
pos <- seq(from =-unique(sampleTable$upstream)
, to = unique(sampleTable$downstream), by = 1)
suppressMessages(c.ggplot2 <- data.frame(pos = pos,
melt(countsLog2,
variable.name= "sample",
value.name= "PosSignal")))
### get mean per group for ratio plots
m.f <- function(x) { smhuber(x)$mu }
m.per.group <- t(aggregate.data.frame(t(countsLog2),
by = list(sampleTable$condition), FUN = m.f )[,-1])
colnames(m.per.group) <- levels(sampleTable$condition)
# introduce post hoc l2fc cutoff
idxSig <- resEdgeR$ldfr < 0.2 & (abs(resEdgeR$logFC / log(2)) > 0.05)
m.per.group <- as.data.frame(m.per.group)
m.per.group$pos <- pos
m.per.group$significant <- idxSig
dataGG <- gather(m.per.group,
key = "experimental_Group_and_significance",
value = "mean_log2_counts", 1:2)
dataGG$experimental_Group_and_significance <- ifelse(dataGG$significant,
"significant",
as.character(dataGG$experimental_Group_and_significance))
dataGG$experimental_Group_and_significance <- mapvalues(dataGG$experimental_Group_and_significance,
from = c("mutant", "WT", "significant"), to = c("Mutant", "Wild type", "Significant Difference"))
dataGG$experimental_Group_and_significance <- factor(dataGG$experimental_Group_and_significance)
pl <- (ggplot(data = dataGG,
aes_string(x = "pos",
y = "mean_log2_counts",
color = "experimental_Group_and_significance")) +
geom_point() +
labs(x ="Distance from TSS (bp)",
y="Normalized counts (log2)") +
scale_color_manual(values = c( "#d95f02", "#7570b3", "black"), name = "") +
theme(panel.background = element_blank(), panel.grid.minor=element_blank(),
axis.line = element_line(colour = "black", size = 0.5)))
pl
#svg(filename = 'K4me2_significance_edgeR.svg',height = 6, width = 15)
(pl + theme( axis.ticks.length = unit(10, "points"),
axis.line = element_line(size = 1),
axis.title = element_text(size = 16, face = "bold"),
axis.text = element_text(size = 12, face = "bold"),
plot.title = element_text(size = 16, face = "bold"),
legend.position = c(0.4, .7),
legend.key = element_rect(fill = "white"),
legend.text = element_text(size = 16, face = "bold")
)
+ ggtitle("Call for significantly enriched regions based on csaw/edgeR"))
#dev.off()
