Hello!

Me and some colleagues are trying to work out if we need to account for multiple testing when using DESEQ2 or if this is something the package covers already. And also, how to handle global testing rather than pairwise.

Our data structure: RNA-seq gene count data with the variables: tissue type (6 levels), treatment (2 levels control/treated), and replicate number (3 levels).

The types of questions we want to ask: (1) which genes are differentially expressed in each tissue between control and treatment? (2) which genes are differentially expressed from one tissue to the next? (3) which genes are differentially expressed in one tissue compared to the others? (global)

On a previous forum, someone suggested creating another variable "group" which combines tissue type and treatment to allow us to extract certain groups needed for the first question. Therefore "group" variable has 12 levels.

Contrast examples:

TissueA:control versus TissueA:treated

Tissue B:control versus TissueB:treated

etc

Tissue A versus Tissue B

Tissue A versus Tissue C

etc

Tissue A versus all other tissues.

Tissue B versus all other tissues.

etc

I am aware that DESEQ2 using BH adjustment to generate p-adjusted values. According to Michael Love "The adjusted p-value column only correct across genes." but not across different contrasts/comparisons. I have also looked in the vignette and have not found an answer.

Could anyone please advise on whether we need to do any extra steps for DESEQ2 to account for multiple testing problems, and if we do, what?

The default is to adjust for multiplicity within a contrast, using BH, so you expect at most a 5% false positive rate for each contrast. If you have a 5% maximum false positive in N contrasts, you also have a 5% max false positive rate over the combination of all contrasts as well.

This elides one thing though - in order to achieve any significant results, you need to achieve Bonferroni significance for at least one test. If you have say 15,000 genes, and for sake of argument, 15 contrasts, that means you need at least one gene with a p<3e-6 in a given contrast to have a chance for any BH significant genes. But if you do the global test by concatenating the p-values from all 15 contrasts, then you need at least one gene to have a p<2e-7. If you don't have that, you are going to say that none of the genes in any contrast are significant.

Anyway, the convention is to just adjust by contrast. I think that's fine. By accepting any false positives, you are by definition saying that the set of significant genes is preliminary in some sense (since you don't know which ones are the false positives after all) and will need some follow-up for validation.