Drawing DESeq2 fitted spline curves on top of scaled counts

deseq2_curves.R

library(splines)
library(DESeq2)

# make some demo data
dds <- makeExampleDESeqDataSet(n=100, m=40)
dds$condition <- sort(runif(40))
# make one gene where expression has a curve shape (just for demo)
s_shape <- round(500 * sin(dds$condition*2*pi) + 1000 + rnorm(40,0,50))
mode(s_shape) <- "integer"
counts(dds)[1,] <- s_shape

# just for visualization
# next line just for demo, use standard size factors for real data
sizeFactors(dds) <- rep(1,ncol(dds))
library(ggplot2)
plotCounts(dds, 1, returnData=TRUE) %>%
  ggplot(aes(condition, count)) +
  geom_point()

# how to fit curves in DESeq2
design(dds) <- ~ns(condition, df=5)
dds <- DESeq(dds, test="LRT", reduced=~1)

# expected results?
res <- results(dds)
hist(res$pvalue)
plot(res$pvalue, log="y")

# build the needed matrices for fitted spline
coef_mat <- coef(dds)
design_mat <- model.matrix(design(dds), colData(dds))

# build the fitted splines and log counts
library(ggplot2)
library(dplyr)
dat <- plotCounts(dds, 1, returnData=TRUE) %>%
  mutate(logmu = design_mat %*% coef_mat[1,],
         logcount = log2(count))

# plot it all together
dat %>%
  ggplot(aes(condition, logcount)) +
  geom_point() +
  geom_line(aes(condition, logmu), col="dodgerblue")