Graphing timecourse data from DESeq2
0
0
Entering edit mode
15 months ago

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.

DESeq • 302 views
ADD COMMENT

Login before adding your answer.

Traffic: 2591 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6