Dr Cheng ZhangDigital Liver Group · translational systems medicine

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 stat and pvalue columns?

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?