Hello,

I am new to DeSeq2 and I am trying to use a multi-factor design to answer this question:

Is captured/input in condition 1 different from captured/input in condition 2 ?

I used a combination of the DeSeq2 manual and DESeq2 testing ratio of ratios (RIP-Seq, CLIP-Seq, ribosomal profiling) to write the following script:

>print(MF_coldata)

                          MF_condition           MF_assay
1_captured_rep1           1                      captured
1_captured_rep2           1                      captured
1_captured_rep3           1                      captured
1_input_rep1              1                      input
1_input_rep2              1                      input
1_input_rep3              1                      input
2_captured_rep1           2                      captured
2_captured_rep2           2                      captured
2_captured_rep3           2                      captured
2_input_rep1              2                      input
2_input_rep2              2                      input
2_input_rep3              2                      input
3_captured_rep1           3                      captured
3_captured_rep2           3                      captured
3_captured_rep3           3                      captured
3_input_rep1              3                      input
3_input_rep2              3                      input
3_input_rep3              3                      input
4_captured_rep1           4                      captured
4_captured_rep2           4                      captured
4_captured_rep3           4                      captured
4_input_rep1              4                      input
4_input_rep2              4                      input
4_input_rep3              4                      input
5_captured_rep1           5                      captured
5_captured_rep2           5                      captured
5_captured_rep3           5                      captured
5_input_rep1              5                      input
5_input_rep2              5                      input
5_input_rep3              5                      input


>ddsMF <- DESeqDataSetFromMatrix(countData=countdata, colData=MF_coldata, design=~ MF_assay + MF_condition + MF_assay:MF_condition)
>ddsMF <- ddsMF[ rowSums(counts(ddsMF)) > 1, ]

>ddsMF <- DESeq(ddsMF, test="LRT", reduced= ~ MF_assay + MF_condition)
>print(results(ddsMF)) # OUTPUT1

>resMF = results(ddsMF, contrast=c("MF_condition", '1','2'))
>print(resMF) # OUTPUT2

My problem is when comparing OUTPUT1 and OUTPUT2, I obtain different values of log2FoldChange but p-values are the same. Here are the headers for the two OUTPUTs:

OUTPUT1:

log2 fold change (MLE): MF assayinput.MF condition2
LRT p-value: '~ MF_assay + MF_condition + MF_assay:MF_condition' vs '~ MF_assay + MF_condition'

OUTPUT2:

log2 fold change (MLE): MF_condition 1 vs 2 
LRT p-value: '~ MF_assay + MF_condition + MF_assay:MF_condition' vs '~ MF_assay + MF_condition'

Am I doing something wrong with the contrast feature?

Thanks in advance!

> sessionInfo( )
R version 4.1.0 (2021-05-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /usr/local/intel/compilers_and_libraries_2020.2.254/linux/mkl/lib/intel64_lin/libmkl_rt.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
[1] compiler_4.1.0

My problem is when comparing OUTPUT1 and OUTPUT2, I obtain different values of log2FoldChange but values are the same. Here are the headers for the two OUTPUTs:

Well, yeah, you are asking for different comparisons. I also don't think the printed L2FC values means much in a LRT test. Are you absolutely sure that LRT is what you want?

What exact comparison are you trying to do?

I strongly recommend you out your normalized counts into something, get averages for your different groups, and figure out what those different designs and ways of calling results are giving you.

See ?results. You need to change test argument to generate new pvalues.