Question: Deseq2 Multivariate Time Series: DEGs for each time point comparison?
0
gravatar for jamesfischer047
6 weeks ago by
jamesfischer0470 wrote:

Hi all,

I have a set of RNA-seq data. It's a time series with 4 stages of leaf development, 2 treatments, and 3 bio-reps per stage x treat combo. I've already found DEGs for the treatments, and DEGs for the stage:treatment interaction, but what I would like to do is find DEGs at each stage comparison (e.g. genes that are specifically differentially expressed at stage 2, etc.). I've tried doing clustering with maSigPro, but it did not really provide the resolution that I need.

Any ideas on this? My code is below.

Thanks in advance!

  reads <- as.matrix(read.csv("genereads_ONLY3.txt", sep = '\t', row.names = 1, header = TRUE)) 
  meta <- read.table("metatest2.txt", header = TRUE)




  mm <- model.matrix(~ stage + treat + stage:treat, data = meta)
    mm2 <- model.matrix(~ stage + treat, data = meta)
      mm3 <- model.matrix(~ stage, data = meta)

    ddsd <- DESeqDataSetFromMatrix(countData = reads,
                                  colData = meta,
                                  design = mm)

    dds <- ddsd[ rowSums(counts(ddsd)) > 10, ]

  #running DESeq

   dds.interaction <- DESeq(dds, test = "LRT", full = mm, reduced = mm2) 
                #Retrieves genes that are DEGs due to interaction of stage and treatment
    dds.interact.results <- results(dds.interaction)
      sum(dds.interact.results$padj < 0.05, na.rm = TRUE)
        sub.interact <- subset(dds.interact.results, padj<0.05)

           dds.interact.DEGs.genes <-  row.names(sub.interact) 
              dds.interact.DEGs.genes.df <-  as.data.frame(dds.interact.DEGs.genes) 
                colnames(dds.interact.DEGs.genes.df) <- "GENES"

    dds.treatonly <-  DESeq(dds, test = "LRT", full = mm, reduced = mm3)
       dds.treat.DEGS <- results(dds.treatonly)
           sum(dds.treat.DEGS$padj < 0.05, na.rm = TRUE)
              sub.treat <- subset(dds.treat.DEGS, padj<0.05)

          dds.treat.DEGs.genes <-  row.names(sub.treat) 
            dds.treat.DEGs.genes.df <- as.data.frame(dds.treat.DEGs.genes)
                colnames(dds.treat.DEGs.genes.df) <- "GENES"

           DIFFERENCE <- setdiff(dds.treat.DEGs.genes.df, dds.interact.DEGs.genes.df)
deseq2 rna-seq • 112 views
ADD COMMENTlink modified 6 weeks ago by Michael Love24k • written 6 weeks ago by jamesfischer0470
Answer: Deseq2 Multivariate Time Series: DEGs for each time point comparison?
1
gravatar for Michael Love
6 weeks ago by
Michael Love24k
United States
Michael Love24k wrote:

Have you seen the time series example in the workflow? You can pull out time/stage specific differences usually using results() with name or contrast=list(..., ...).

ADD COMMENTlink written 6 weeks ago by Michael Love24k

Thank you for the answer! Whenever I use resultsNames(), I get the same thing for my different tests (the different tests have different designs, as can be seen above).

I get:

[1] "Intercept"       "stage"           "treatX800"       "stage.treatX800"

In which "X800" represents my treatment.

I know I've seen other outcomes of "results()" online that seem much lengthier. How would I go about designing my 'contrast = " with that output to get a comparison of a single stage between my control and treat?

Here is my full model matrix in case that helps:

  (Intercept) stage treatX800 stage:treatX800
1            1     1         0               0
2            1     1         0               0
3            1     1         0               0
4            1     2         0               0
5            1     2         0               0
6            1     2         0               0
7            1     3         0               0
8            1     3         0               0
9            1     3         0               0
10           1     4         0               0
11           1     4         0               0
12           1     4         0               0
13           1     1         1               1
14           1     1         1               1
15           1     1         1               1
16           1     2         1               2
17           1     2         1               2
18           1     2         1               2
19           1     3         1               3
20           1     3         1               3
21           1     3         1               3
22           1     4         1               4
23           1     4         1               4
24           1     4         1               4
ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by jamesfischer0470
1

When you run DESeqDataSetFrom...() it should be printing a message to the console that you have encoded the variables as numeric, and that you may have intended them to be factors instead. You should heed that message.

(Edit: the message comes upon dataset construction, not DESeq())

ADD REPLYlink modified 6 weeks ago • written 6 weeks ago by Michael Love24k

Thank you! For some reason, that message never popped up. But that was in fact the issue!

ADD REPLYlink written 6 weeks ago by jamesfischer0470

Here is what happens when you create a dds with an integer-valued numeric in the design:

> DESeqDataSetFromMatrix(counts(dds), colData(dds), ~condition)
  the design formula contains a numeric variable with integer values,
  specifying a model with increasing fold change for higher values.
  did you mean for this to be a factor? if so, first convert
  this variable to a factor using the factor() function

What version of DESeq2 are you using?

ADD REPLYlink written 6 weeks ago by Michael Love24k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 270 users visited in the last hour