Dr Cheng ZhangDigital Liver Group · translational systems medicine

Bioinformatics, step by step · Workshop

Module 6 · Multiple testing and the false discovery rate

Video: T8 20,000 tests at once. Script: scripts/06_multiple_testing.R (run module 5 first)

Goal. See false positives appear when nothing has changed. Then reproduce Benjamini–Hochberg (BH) by hand, and understand the NAs in padj.

A. A null experiment

We split the four SHAM samples into two fake groups, each with one day-1 and one day-3 sample.

Q6.1 How many genes have a raw p < 0.05? How many have padj < 0.05? What does the p-value histogram look like?

Q6.2 Now put both day-1 samples in group A and both day-3 samples in group B. What changes, and why? What does this tell you about confounding?

B. The real comparison

Q6.3 How many tests were run at day 1? If no gene had changed, how many raw p < 0.05 would you expect? How many did we get?

Q6.4 Describe the shape of the p-value histogram. The day 3 histogram also has a spike near p = 1. What kind of genes might produce that? (Hint: look at baseMean for genes with p > 0.95.)

C. BH by hand

Q6.5 Run bh_by_hand() on the ten toy p-values from the video. Does it match p.adjust(p, "BH")? How many pass 0.05 with raw p, with Bonferroni and with BH?

D. Independent filtering

Q6.6 p.adjust(r1$pvalue, "BH") does not give the same number of significant genes as DESeq2's padj. Why? Show that the two do agree once the filtered genes are removed.

Q6.7 How many genes have NA in pvalue? How many have a p-value but NA in padj? What is different about them?