Hi,

I'm analyzing 36 RNAseq samples in order to obtain differentially expressed genes. The study design is a single time course with 12 timepoints ("t1" to "t12") and 3 replicates per timepoint (A, B and C)

First of all I want to confirm the strategy I'm following it's the correct one. I'm using edgeR-glm for doing the comparisons 't2-t1', 't3-t2', 't4-t3', 't5-t4' and so on. What I pretend here is to obtain the DE genes between each time. I know about the existence of packages prepared to do time course analyses, such as 'maSigPro'. However I had a problem running this analyses and I think it's due to the large amount of timepoints. Specifically, the command that fails is:

design <- make.design.matrix(eDesign, degree = 2)
min.obs <- (2+1)*(ncol(eDesign)-2)+1
#command that fails:
fit <- p.vector(data = counts_norm, design = design, Q = 0.05, counts = TRUE, min.obs = min.obs) 

Here a little snippet of the 'eDesign' matrix:

I have tried up to degree = 5, and still failing. I think the problem is the timepoints number (12) because if I group them in only 5 the command runs smoothly.

That's the major cause I'm doing this pairwise comparisons between times ('t2-t1', 't3-t2', 't4-t3'...), cause I can't directly use maSigPro with the 12 timepoints. Moreover, with the edger-glm pairwise approach I can know exactly what genes are DE between each paire of times. Haven't seen this is possible with maSigPro.

For this 5 timepoint grouping I have been based on the plotMDS results and the DE analyses between times. And here comes my next question:

- From what data should MDS plots be done: raw counts or normalized counts? I would say normalized data cause they take into account lib. sizes. The point is I don't obtain the same results when running next commands:

# First one
DGE <- calcNormFactors(DGE, method = "TMM")
plotMDS(DGE)

# Second one
DGE_norm <- tmm(DGE$counts, long = 1000, lc = 0, k = 0)
plotMDS(DGE_norm)

I know calcNormFactors()  doesn't actually normalize the counts, just calculates the normalization factor. However in the edgeR manual (https://www.bioconductor.org/packages/devel/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf), sections 4.1.4 and 4.1.5, the pass to limma::plotMDS() this normalized DGElist, such as the first comand above. Doing that I obtain the same MDS plot that with the raw counts (the raw DGElist without doing calcNormFactors()). The second command does generate a different plot. I could attach the plots if someone ask for them.

I'm interested in knowing how my samples group because I need to reduce the timepoints from 12 to 4-6 in order to use maSigPro.

Many thanks,

Dani

I can't speak for whatever MaSigPro is doing, but edgeR itself can handle 12 time points without any problems. There is no need to reduce the number of time points prior to a DE analysis. In your case, it seems you want to do comparisons between specific time points, so parametrizing each time point as a group would be sufficient:

design <- model.matrix(~0 + Time, data=eDesign)

... or, if the samples with the same Replicate value are related to each other:

design <- model.matrix(~0 + Time + Replicate, data=eDesign)

You should then be able to do comparisons between specific groups (time points) as you please, using edgeR's GLM machinery to set up and test contrasts of interest. In fact, you seem to acknowledge that this is possible in your post, so I don't see why you need to consolidate the time points to get what you want.

As for the MDS plots; I don't know what the tmm() function does, so I can't speak for that. All I can say is that calling plotMDS() on a DGEList will first call cpm() and then create an MDS plot on the log-CPMs generated from the count matrix. This will use any normalization factors that are present, so the result will be different from that generated with the "raw" DGElist. If you're not observing any differences, I would suggest that it is because your normalization factors are so close to 1 that they're not having much effect.

Hi Aaron, many thanks for your quick response. I must apologize for one thing: MDS plots do actually vary due to the TMM normalization. The samples distribution is quite similar in both plots, and this is why I got mistaken, just didn't look enough to them. 

Here both plots. First number means the timepoint and letter is the replicates (A, B or C). I don't know if the replicates pertain to different batches. Last number (i.e. 1A1, 4B2) has not importance.

Here a summary of lib. sizes and norm. factors:

> summary(DGE$samples$lib.size)
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
 8232274  9206213 11469978 11822144 13279158 18278786
> summary(DGE$samples$norm.factors)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
 0.5085  0.8216  1.0399  1.0521  1.3012  1.5670

And lastly, plots showing counts distribution before (first one) and after normalization (second one).

Based on these initial quality control plots I would say data quality seems to be quite good. However, I'm thinking in removing the '12C1' as it clearly doesn't group with the other 't12' samples.

There is something you said that I have not clear. When you say '... or, if the samples with the same Replicate value are related to each other...', what do you mean exactly? As far as I understand these replicates are supposed to be identical (in theory at least), so I'm not taking into account the Replicates information. Should I include this information in my analysis?

Once normalization and dispersion stimation have been done, that's how I'm doing the contrasts (notice I build design in a different manner):

# times are 't1', 't2', 't3'...
group <- as.factor( paste("t", (rep(c(1:12), each=3)), sep=""))
# design creation
design <- model.matrix(~0+group)
colnames(design) <- gsub("group", "", colnames(design))
# Contrasts
contr.matrix <- makeContrasts(
  t2vst1 = t2-t1,
  t3vst2 = t3-t2,
  t4vst3 = t4-t3,
  t5vst4 = t5-t4,
  t6vst5 = t6-t5,
  t7vst6 = t7-t6,
  t8vst7 = t8-t7,
  t9vst8 = t9-t8,
  t10vst9 = t10-t9,
  t11vst10 = t11-t10,
  t12vst11 = t12-t11,
  levels = colnames(design)
)
# DE analysis
fit <- glmQLFit(DGE, design)
qlf <-glmQLFTest(fit, contrast = contr.matrix[,"t2vst1"])

I run topTags(qlf, n = nrow(DGE$counts))$table after every constrast, and that's how I obtain DE genes for every contrast. This is the way I have done my time course analysis, but I'm not sure if is the best, or even if it's correct.

Is there any better way than doing these pairwise comparisons?