This is a test version of Biostars. For the public version, visit https://www.biostars.org.
[DESeq2 multiple treatments vs. 1-by-1] Inconsistant p-adj

So we have RNAseq data for 3 treatments (A, B and C) that we compared to control. Each treatment should be compared to control (as in A vs. CT, B vs. CT, C vs. CT) and from then we should know the number of differentially expressed genes (padj < 0.1) and which genes are upregulated and downregulated.

When I run the DESeq2 pipeline with a table containing Control (CT) and the 3 treatments, we found about 1,000 differentially expressed genes (padj < 0.1) but the result CSV spreadsheet only had one log2fold column and padj column, where I was expecting to have 3 column log2fold and padj for each comparison (A vs. CT, B vs. CT and C vs. CT). Since I couldn't extract the log2fold and padj for each comparison, I re-run the DESeq2 pipeline on each treatment vs control separately, and now the number of differentially expressed genes is completely different, I obtained:

  • A vs CT : 1000 genes
  • B vs CT: 2000 genes
  • C vs CT : 2500 genes

In total that's far beyond the original 1,000 genes I had when running DESeq2 with the 3 treatments vs Control in one spreadsheet.

  1. What causes this difference?
  2. And now I'm wondering which pipeline is the right one?
  3. Should I run it independently or all treatments together on one spreadsheet?
  4. If so, how can I extract the padj and log2 fold changes for each treatment if they are run together??
rna-seq deseq2
• 4,967 views • link
updated by Ram Ram • written by madkitty madkitty

You'll need to show the code you used in the first vs. one of the second cases for us to help. My guess is that in the first instance you ended up comparing the full model against something like ~1, which isn't what you want. It's completely possible to extract fold-changes and adjusted p-values while keeping everything in. In fact, you'll get more reliable results that way, due to better variance estimation.

• 0 views • link
updated by Ram Ram • written by Devon Ryan Devon Ryan

2 answers

You want to use the "contrast" argument to the results() function, in order to build a results table for your three desired contrasts. You can run DESeq() on the dataset containing all the samples, which can improve variance estimation as Devon mentioned.

First read over the section on contrasts in the DESeq2 vignette and the help for the "contrast" argument in ?results.

Also, as Devon suggested, it helps to provide your code to get more precise answers to your questions.

• 0 views • link
updated by Ram Ram • written by Michael Love Michael Love

I too have the same of similar doubt ,so im comparing multiple cell types each having its own control sample based on cell hierarchy, for normalisation purpose i use all the sample for pca ,clustering ,correlation etc , but when it comes to doing a differential expression i have to do each of them separately ,like Stem cell vs progenitor , Common myeloid progenitor vs Granulocyte monocyte progenitor[GMP] , then GMP vs Monocyte .

Its a conceptual doubt since my control is not always the same as in one case its stem cell and other one is progenitor cell .So how can i make multiple in this case as each comparison will have different foldchange and the calculated p value ,unless i have my control same for every test.

I have used the contrast when i was making comparison stem cell with everything downstream ,but how to do when controls are different

• 1 views • link
written by 1769mkc 1769mkc

I have a data set of six sinalge samples with no replica. I know, that DESeq can work with these, as I have done it already. I would like to know, what to do in this case, when working with multiple comparisons.

I have ran the analysis in a pair-wise mode. DESeq than gives this output:

estimating size factors # workflow of the DESeq2 function pair-wise modus
estimating dispersions
same number of samples and coefficients to fit, estimating dispersion by treating samples as replicates
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

But when doing the DESeq() command with all the samples in one go, I sort of trick DESeq to think there are multiple replica in the data set and I get this output:

estimating size factors #multiple comparisons in one count table
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

is it better in this case to run DESeq in a pair-wise modus?

Don't I sort of cheat the results, when letting DESeq thinks there are multiple replica in the data set?

thanks

Assa

• 0 views • link
updated by Ram Ram • written by Assa Yeroslaviz Assa Yeroslaviz

Please post this as a new question rather than as an answer.

• 1 views • link
written by Devon Ryan Devon Ryan

Log in to answer this question.

More posts like this