QC and Filtering

A count matrix fresh out of a pipeline is mostly genes with counts at or near zero in every sample. Those rows carry no information about differences between conditions, and they add noise to the variance estimate. The first step is to drop them. The second is to look at how the samples relate before any model is fit, when a swapped label or an outlier is still easy to catch.

Filter low-count genes

The common rule keeps a gene when it has at least 10 reads summed across all samples. The threshold applies to the raw counts, before normalization, and it works as a floor rather than a parameter to tune. A gene under that floor has too few observations for a differential expression test to have any power.

library(DESeq2)
library(airway)

data(airway)
dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds$dex <- relevel(dds$dex, ref = "untrt")

# Keep genes with at least 10 reads total.
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]

# Genes before and after filtering.
nrow(airway)
nrow(dds)

[1] 63677
[1] 22369

The design formula ~ cell + dex asks DESeq2 to estimate the treatment effect while accounting for the donor cell line. Calling relevel sets untreated as the reference level, so the reported coefficient compares treatment against untreated.

Variance-stabilizing transformation

Raw counts make a poor input for PCA. Samples differ in sequencing depth, so the counts are not directly comparable, and the variance of a count grows with its mean, which lets the most highly expressed genes dominate every distance. The variance-stabilizing transformation in vst removes that mean-variance dependence, so a distance between two samples reflects the samples themselves rather than the size of their libraries.

# Variance-stabilizing transform, blind to the design for QC.
vsd <- vst(dds, blind = TRUE)

# Transformed dimensions match the filtered gene set.
dim(assay(vsd))

[1] 22369     8

With blind = TRUE, the transformation is estimated without reference to the design, which is the right choice for an unsupervised look at the samples.

PCA

When the samples separate by condition on a PCA of the transformed counts, the treatment has a detectable effect. A sample sitting far from its own group points to a swapped label or a bad library.

# PCA data: one row per sample. plotPCA returns every column of the sample
# metadata, so select the two axes and the covariates that the design uses.
pca_data <- plotPCA(vsd, intgroup = c("dex", "cell"), returnData = TRUE)
head(pca_data[, c("PC1", "PC2", "group", "dex", "cell")])

# The fraction of variance on each axis.
attr(pca_data, "percentVar")

                  PC1        PC2         group   dex    cell
SRR1039508 -20.261251 -4.3138021  untrt:N61311 untrt  N61311
SRR1039509   9.600245 -0.9134282    trt:N61311   trt  N61311
SRR1039512 -11.226224 -4.4215017 untrt:N052611 untrt N052611
SRR1039513  20.313084 -2.5411925   trt:N052611   trt N052611
SRR1039516 -17.032150 17.8512194 untrt:N080611 untrt N080611
SRR1039517  11.716744 23.1974284   trt:N080611   trt N080611
[1] 0.4094782 0.2613781

The percentVar attribute holds the fraction of the total variance that each principal axis carries. A first axis that separates treated from untreated samples, and that carries a large fraction of the variance, is the signal that the treatment moved expression across the transcriptome. The visualization page draws this PCA as a figure; the numbers here are the evidence behind it.