Dear all,

I am trying to use the DESeq2 to conduct the differential expression analysis in time course experiments. I want to find changes of gene expression at different time points after compound treatment. Time point is the only phenotype to be included. I followed an example at http://www.bioconductor.org/help/workflows/rnaseqGene/#time-course-experiments.

My group design contains 0, 0.5, 1, 2, 4, 8, 12, and 24 h, and N=3 for each group.

dds <- DESeqDataSetFromMatrix(countData = expr, colData = coldata, design= ~ Time)

dds_diff <- DESeq(dds, test="LRT", reduced=~1)

time_levels <- levels(dds_diff$Time)

> time_levels
#[1] "0h"   "0.5h" "1h"   "2h"   "4h"   "8h"   "12h"  "24h" 

# Loop through each pair of time levels and perform pairwise comparisons

for (i in 1:(length(time_levels) - 1)) {
  for (j in (i + 1):length(time_levels)) {
    contrast <- c("Time", time_levels[j], time_levels[i])
    res <- results(dds_diff, test = "Wald", contrast = contrast)
    # extract differential results
    res <- res[order(res$padj),]
    pairwise_results[[paste(time_levels[j], "vs", time_levels[i])]] <- res
  }
}

head(coldata)
#GS_ID  Num MaterialType    RQN Time
#s_016  1   total RNA   7.5 1h
#s_017  2   total RNA   7.6 8h
#s_018  3   total RNA   7.4 24h
#s_019  4   total RNA   7.3 12h
#s_020  5   total RNA   7.6 2h
#s_021  6   total RNA   7.6 0.5h
#s_022  7   total RNA   7.6 0.5h
#s_023  8   total RNA   7.4 8h
#s_024  9   total RNA   7.5 0h
#s_025  10  total RNA   7.7 4h
#s_026  11  total RNA   6.4 1h
#s_027  12  total RNA   7.8 12h

contrast <- c("Time", time_levels[2], time_levels[1])
> contrast
#[1] "Time" "0.5h" "0h"

Actullay, I have also used the limma-trend pipline to complete the analysis, and I want to compare whether I can get some similar significantly different genes by the use of these two piplines.

The comparison levels setting in limma-trend is here. input was generated by using DGEList()

treatment = factor(coldata$Time, level=c("0h","0.5h","1h","2h","4h","8h","12h","24h"))
design <- model.matrix(~ 0 + treatment)

logCPM <- cpm(input, log = TRUE, prior.count = 3)
fit <- lmFit(logCPM, design)
fit2 = contrasts.fit(fit,my.contrasts)
fit.eb = eBayes(fit2, trend = TRUE, robust= FALSE)

my.contrasts = makeContrasts(
  # versus baseline: 1:7
  t05h_vs_t0h = t05h - t0h, 
  t1h_vs_t0h  = t1h - t0h, 
  t2h_vs_t0h  = t2h - t0h, 
  t4h_vs_t0h  = t4h - t0h, 
  t8h_vs_t0h  = t8h - t0h, 
  t12h_vs_t0h = t12h - t0h, 
  t24h_vs_t0h = t24h - t0h, 
  # vs 0.5h: 8:13
  t1h_vs_t05h = (t1h - t0h) - (t05h - t0h),
  t2h_vs_t05h = (t2h - t0h) - (t05h - t0h),
  t4h_vs_t05h = (t4h - t0h) - (t05h - t0h),
  t8h_vs_t05h = (t8h - t0h) - (t05h - t0h),
  t12h_vs_t05h = (t12h - t0h) - (t05h - t0h),
  t24h_vs_t05h = (t24h - t0h) - (t05h - t0h),
  # vs 1h: 14:18
  t2h_vs_t1h = (t2h - t0h) - (t1h - t0h),
  t4h_vs_t1h = (t4h - t0h) - (t1h - t0h),
  t8h_vs_t1h = (t8h - t0h) - (t1h - t0h),
  t12h_vs_t1h = (t12h - t0h) - (t1h - t0h),
  t24h_vs_t1h = (t24h - t0h) - (t1h - t0h),
  # vs 2h: 19:22
  t4h_vs_t2h = (t4h - t0h) - (t2h - t0h),
  t8h_vs_t2h = (t8h - t0h) - (t2h - t0h),
  t12h_vs_t2h = (t12h - t0h) - (t2h - t0h),
  t24h_vs_t2h = (t24h - t0h) - (t2h - t0h),
  # vs 4h: 23:25
  t8h_vs_t4h = (t8h - t0h) - (t4h - t0h),
  t12h_vs_t4h = (t12h - t0h) - (t4h - t0h),
  t24h_vs_t4h = (t24h - t0h) - (t4h - t0h),
  # vs 8h: 26:27
  t12h_vs_t8h = (t12h - t0h) - (t8h - t0h),
  t24h_vs_t8h = (t24h - t0h) - (t8h - t0h),
  # vs 12: 28
  t24h_vs_t12h = (t24h - t0h) - (t12h - t0h),
  levels = design
)
head(my.contrasts)[1:6,3:10]
# Contrasts
#Levels t2h_vs_t0h t4h_vs_t0h t8h_vs_t0h t12h_vs_t0h t24h_vs_t0h t1h_vs_t05h t2h_vs_t05h
#t0h          -1         -1         -1          -1          -1           0           0
#t05h          0          0          0           0           0          -1          -1
#t1h           0          0          0           0           0           1           0
#t2h           1          0          0           0           0           0           1
#t4h           0          1          0           0           0           0           0
#t8h           0          0          1           0           0           0           0

Here are my questions:

1. Did I use the correct design and test methods in both DESeq2 and limma-trend?
2. I want to know whether the comparison levels I set in DESeq2 were consistent with those in the limma-trend pipeline?
3. Should I also use makeContrasts() in DESeq2 to set the same comparison levels as in the limma-trend?

Thank you for your helping!

Jiahao Tian

If your question is "what genes are changed in any of the time points?" LRT is better than sifting through every possible time comparison.

Are you sure you want every time point to every time point? I would think you might want every time point to time 0, or every time point to the time before it. Or maybe cluster the line graphs, to see what time patterns are present.

Your limma-trend analysis seems correct in principle, although your contrasts are repeating t0h unnecessarily. For example, defining

t1h_vs_t05h = (t1h - t0h) - (t05h - t0h)

is exactly the same as:

t1h_vs_t05h = t1h - t05h