Comparing results from two different DESeq analyes

Question

[luzglongoria]

Hi there,

I am analysing some RNA expression from an experiment with a set up like:

sample     condition       moment
    R1         A             1     
    R2         A             1    
    R3         A             1     
    R4         A             2     
    R5         A             2     
    R6         A             2  
    C1         B             1     
    C2         B             1    
    C3         B             1     
    C4         B             2     
    C5         B             2     
    C6         B             2

I have compared RNA expression between moments. It is mean : (R1,R2,R3) vs (R4,R5,R6)
and (C1,C2,C3) vs (C4,C5,C6) by doing:

library(DESeq2)
library(tidyverse)

#### Load data
library(readxl)

setwd("~/Documents/path/to/txt/file/")
data= read.table("Expression_level.txt", header = T)
View(data)
                          R1   R2     R3       R4     R5       R6
gene-CpipJ_CPIJ008101 484021 412077 445173  154707  148776  169263
gene-CpipJ_CPIJ001132 334997 391789 435968  445623  504466  445865
gene-CpipJ_CPIJ006209 326414 260289 301946  169859  149214  141446
gene-CpipJ_CPIJ002271 320207 282722 326901  203648  170398  134834
gene-CpipJ_CPIJ005941 316818 252593 273103   55266   43730   26304
gene-CpipJ_CPIJ009303 269236 357244 386633  426546  531801  483546
gene-CpipJ_CPIJ010326 233568 226659 254108  362953  278742  325969
gene-CpipJ_CPIJ008915 230936 276916 277624  355937  357974  239651
gene-CpipJ_CPIJ009571 223388 187980 207711  128457  139515   87437

annotation.info <- read.table("~/Documents//path/to/txt/file/",header = T)

  Viewannotation.info)
  sample     condition       moment
    R1         A             1     
    R2         A             1    
    R3         A             1     
    R4         A             2     
    R5         A             2     
    R6         A             2  

## Create Data Set
dds <- DESeqDataSetFromMatrix(countData = data,
                              colData = annotation.info,
                              design = ~ moment)

#do the analyses
dds <- DESeq(dds)
res <- results(dds)
res

And then I did the same for Condition B ---> res.1

Obviously, the results are different for each analysis:

summary(res) 
out of 18958 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 2108, 11%
LFC < 0 (down)     : 1873, 9.9%
    outliers [1]       : 2, 0.011%
    low counts [2]     : 2911, 15%
    (mean count < 2)

summary(res.1)
    out of 17783 with nonzero total read count
    adjusted p-value < 0.1
    LFC > 0 (up)       : 4636, 26%
    LFC < 0 (down)     : 5751, 32%
    outliers [1]       : 6, 0.034%
    low counts [2]     : 2039, 11%
    (mean count < 1)

My question is the following: Is there any way of statistically compare these two results? I would like to know which genes are differentially expressed between condition. I would also like to know whether these diferencial expressed genes are up or down regulated. Is it possible?

Thank you so much in advance.

Answer 1

[i.sudbery]

The correct approach here is not to compare the two results, but to use a new design that specifically asks the question you are interested in.

You are interested in whether the different between the time dependent change is different between the different conditions. There are two ways that people conceptualize this.

Combing factors into groups

The limma/edgeR manual has long advocated using an approach where you combine your two columns to create a single "group" column with four different values - A1, A2, B1 and B2.

Your first experiment in your post how calculated the log Fold Change and its significance A1-A2. In the second you do B1-B2. The fold change you are looking for is (B1-B2) - (A1-A2). You would retrieve this by using constrast = list(c("B1", "A2"), c("B2", "A1")) in your call to results.

Adding an interaction term.

A statistician would say that your experimental design is an interaction design, that that the formula should be specified ~condition + moment + condition:moment. In this design a logFoldChange is fitted for the difference between 1 and 2 in A, and then a second logFoldChange is fitted for the difference between logFoldChanges for A and B conditions. See the "interactions" part of the DESeq2 manual: http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#interactions