Bioinformatics, step by step · Workshop
Module 5 · Differential expression with DESeq2
Video: T7 Testing one gene: the GLM and the Wald test. Script: scripts/05_differential_expression.R
Goal. Run DESeq2 for MI vs SHAM at day 1 and at day 3. Then take one gene apart and reproduce its Wald test by hand.
dds <- DESeqDataSetFromMatrix(countData = count, colData = coldata, design = ~conds)
dds <- DESeq(dds)
res1 <- results(dds, contrast = c("conds", "MI_1D", "SHAM_1D")) # numerator first, reference second
Q5.1 Why do we give DESeq2 counts and not TPM? Why do we
round()the counts first?Q5.2 At padj < 0.05, how many genes go up and how many go down at day 1? At day 3? Is the effect of MI larger at day 1 or at day 3?
Q5.3 What are the most significant genes at day 1 and at day 3? Look one of them up. Does it make sense in a heart after an infarction?
One gene by hand
For Serpine1:
Q5.4 Compute z = log2FoldChange ÷ lfcSE, and p = 2 × pnorm(−|z|). Do they match the
statandpvaluecolumns?Q5.5 Compare the log2 fold change with log2(mean of MI_1D ÷ mean of SHAM_1D), using the normalised counts.
Q5.6 Repeat for Nppa and Postn. Which are significant? Use
plotCounts()to explain why.
Plots
Q5.7 In the MA plot, where are the genes with the largest fold changes? Why do low-count genes show such extreme fold changes, and what does
lfcShrink()do about it?Q5.8 In the volcano plot, find a gene with a big fold change that is not significant. Why isn't it?