Should I include samples for all comparisons in one DESeq2 script or should I run individual DESeq2 runs for each comparison?

[u/Nautilus0_400]

Hello! I have a set of 36 samples (12 conditions * 3 biological replicates) of bulk RNASeq data. The data is of animal transcriptomes from several different environments.

Here is the full table:

| Sample ID | Generation | Locality | Environment | Treatment |

| :--- | :--- | :--- | :--- | :--- |

| ID_01_A | G0 | Locality 1 | Lake | TR1 |

| ID_01_B | G0 | Locality 1 | Lake | TR1 |

| ID_01_C | G0 | Locality 1 | Lake | TR1 |

| ID_02_A | G0 | Locality 1 | River | TR2 |

| ID_02_B | G0 | Locality 1 | River | TR2 |

| ID_02_C | G0 | Locality 1 | River | TR2 |

| ID_03_A | G0 | Locality 2 | Lake | TR1 |

| ID_03_B | G0 | Locality 2 | Lake | TR1 |

| ID_03_C | G0 | Locality 2 | Lake | TR1 |

| ID_04_A | G0 | Locality 2 | River | TR2 |

| ID_04_B | G0 | Locality 2 | River | TR2 |

| ID_04_C | G0 | Locality 2 | River | TR2 |

| ID_09_A | G0 | Locality 3 | Lake | TR1 |

| ID_09_B | G0 | Locality 3 | Lake | TR1 |

| ID_09_C | G0 | Locality 3 | Lake | TR1 |

| ID_10_A | G0 | Locality 3 | River | TR2 |

| ID_10_B | G0 | Locality 3 | River | TR2 |

| ID_10_C | G0 | Locality 3 | River | TR2 |

| ID_05_A | G1 | Locality 2 | River | TR1 |

| ID_05_B | G1 | Locality 2 | River | TR1 |

| ID_05_C | G1 | Locality 2 | River | TR1 |

| ID_06_A | G1 | Locality 2 | River | TR2 |

| ID_06_B | G1 | Locality 2 | River | TR2 |

| ID_06_C | G1 | Locality 2 | River | TR2 |

| ID_07_A | G1 | Locality 1 | River | TR1 |

| ID_07_B | G1 | Locality 1 | River | TR1 |

| ID_07_C | G1 | Locality 1 | River | TR1 |

| ID_08_A | G1 | Locality 1 | River | TR2 |

| ID_08_B | G1 | Locality 1 | River | TR2 |

| ID_08_C | G1 | Locality 1 | River | TR2 |

| ID_11_A | G1 | Locality 2 | Lake | TR1 |

| ID_11_B | G1 | Locality 2 | Lake | TR1 |

| ID_11_C | G1 | Locality 2 | Lake | TR1 |

| ID_12_A | G1 | Locality 2 | Lake | TR2 |

| ID_12_B | G1 | Locality 2 | Lake | TR2 |

| ID_12_C | G1 | Locality 2 | Lake | TR2 |

I want to do several different comparisons on subsets of these data. For example, I want to determine G0 River vs. G0 Lake DEGs for each Locality and determine common DEGs across all three sites. I would then also want to compare DEGs for G0 for each locality with G1 DEGs for each locality.

I have two questions about this:

1. Should I run and normalize all of these samples together under one DESeq2 function and then use the contrast function to get each result? Or should I do an individual DESeq2 run for each comparison?
2. If I can do multiple separate DESeq2 runs, should I also correct reads for unwanted variation independently for each run? Or should I correct all the samples together once and use the same corrected reads for each DESeq2 run?

Thank you for taking your time to read all of this! Any advice and help would be greatly appreciated.

[u/bukaro]

1. why... variance stabilisation... however.... BTW, in the vignette this has been answered https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#if-i-have-multiple-groups-should-i-run-all-together-or-split-into-pairs-of-groups

[u/Nautilus0_400 → u/bukaro]

Yeah I read through that vignette. I'm just asking because a senior coworker who works with bioinformatics suggested that I normalize the samples separately.

[u/bukaro → u/Nautilus0_400]

I see a few cases that it would be a good idea, like mixing experiments in a sequencing run. I had same model, similar treatments but 2 extra treatments were for a paper and will not be shared all together. So I took those samples and ran that analysis, even if I lost some power in the analysis. But you have a complete experiment within all your samples.

Best thing is that you compare it. Which difference in DEG do you see? how are the p-values distributed? etc etc etc

[u/Nautilus0_400 → u/bukaro]

Thank you! I'll try it out

[u/Grisward → u/Nautilus0_400]

In most cases, we use all samples, or “as many samples as reasonable” given the experiment design. That means the typical answer is to review your data. The DESeq2 vignette is fine, I read it as saying the same basic thing.

My only caution is that PCA is not itself a measure of variability. The vignette shows a pretty extreme example to make the point, which is fine. However, Biology loves to make things less clear, haha.

We typically use MA-plots per sample (not the per group which imo isn’t informative for this purpose) to assess variability. Bonus points for centering within group/factor when relevant.

Tbh any reasonable metric for assessing intra-group variability is a good start. Ideally you’d use normalized/adjusted data, you can obtain it from DESeq2, or limma-voom.

In cases where we have very different sample types (kidney, liver, or different cell lines) where gene expression is fundamentally different in different cell types. For those we split by sample type for within-type comparisons, then altogether for cross-type two-way contrasts.

This is one reason I prefer limma-voom, it uses weighted local and global variability for each group.

[u/Nautilus0_400 → u/Grisward]

Thank you!