I've got a set of samples, treated and control across multiple time points. I did DESeq with Day as design to get the genes which change nicely across the time span. I wanted to make a graph of some of those genes. The Log2 of the normalized counts, versus the time, colored by treatment, with the regression lines for each data set added.
Here's my code. I know it's not doing exactly what DESeq is doing, but it's close enough for a picture
treated_expression <- log(counts(dds_treated, normalized = TRUE)[gene_of_interest,],2)
control_expression <- log(counts(dds_Control, normalized = TRUE)[gene_of_interest,],2)
lm_treated <- lm(treated_expression ~ dds_treated@colData$Day)
lm_control <- lm(control_expression ~ dds_Control@colData$Day)
d<- plotCounts(dds, gene=gene_of_interest, intgroup=c("Treatment", "Day"), returnData = TRUE)
p <- ggplot(d, aes(x=Day, y=log(count,2), color = Condition)) +
geom_point(position=position_jitter(w=0.1,h=0)) +
labs(title = "gene_of_interest") +
theme(plot.title = element_text(hjust = 0.5))
p + geom_abline(intercept = lm_treated$coefficients[1], slope = lm_treated$coefficients[2], color = "turquoise") +
geom_abline(intercept = lm_control$coefficients[1], slope = lm_control$coefficients[2], color = "red")
And it's working more or less as I want. Unless a sample has a read count of 0. I know I could add a pseudocount before logging, or just omit the zero counts, but I was wondering if there was a a way of fixing this closer to what DESeq does. DESeq results will return me the right slope for the regression line, but I don't see how to find the intercept, which is why I recalculate the regression line.
