Example of ORA and GSEA

The airway dataset¶

The methods of the previous chapter are easiest to appreciate on data where we already know what the answer should look like. We will therefore use the airway dataset of Himes et al. (2014), deposited at GEO as GSE52778.

Himes and co-workers cultured primary human airway smooth muscle (ASM) cells from four donors and treated them with dexamethasone, a synthetic glucocorticoid used to control asthma, or left them untreated. Each of the four cell lines contributes one treated and one untreated sample, so the experiment is a paired two-condition design with $4+4=8$ RNA-seq libraries. This is about as controlled as a human transcriptomics experiment gets: the only thing that differs within a pair is the drug.

It is also an experiment with a built-in positive control. Dexamethasone acts through the glucocorticoid receptor (NR3C1), a transcription factor with a well characterized set of direct target genes – FKBP5, PER1, TSC22D3 and DUSP1 among them. If our pathway analysis is working, that signal should be impossible to miss.

GEO distributes the study as a gene-level FPKM matrix with one column per library, next to annotation and per-condition summary columns. We download it once and cache it under dsbook/data/.

The full series also contains cells treated with the $\beta_2$-agonist albuterol, alone and in combination with dexamethasone. We ignore those and keep the eight libraries of the untreated/dexamethasone comparison. The library names carry the treatment, and the cell lines were sequenced in blocks of four, so LL01–LL04 come from donor N61311, LL05–LL08 from N052611, and so on. That lets us line up the samples pairwise.

Two remarks on that filter. First, the matrix is quantified as FPKM rather than as raw counts; FPKM is a poor basis for a count-based differential-expression model, but the $t$-test we use below and the signal-to-noise metric that GSEA ranks the genes by are both computed from the log-transformed expression values themselves and do not assume count data, so it serves our purpose. Second, the gene-level FPKM values contain a handful of all-or-nothing artefacts, where a gene is reported at several thousand FPKM in one condition and at exactly zero in the other. Requiring at least 1 FPKM in all eight libraries removes them and leaves us with around eleven thousand robustly expressed genes.

Differential expression, dexamethasone versus untreated¶

Because the samples are paired by cell line, we use a paired $t$-test rather than the two-sample test used elsewhere in this book: each gene is tested on the four within-donor differences. This removes the donor-to-donor variation, which in primary cells is substantial.

A volcano plot shows a clear response: some four hundred genes move by more than a factor of two, with about twice as many going up as down.

Before going further, we check that the four glucocorticoid-receptor target genes are up-regulated by dexamethasone.

All four are up by roughly a factor of 8–16 with $q$ values around 1%. The experiment worked.

Note also how modest the $q$ values are. With four pairs the $t$ statistic rests on three degrees of freedom, so even a large change that is fairly consistent across the donors cannot be declared significant with much confidence: the smallest $q$ value in this experiment is about $2 \cdot 10^{-3}$, and only some two hundred genes reach a $q$ value below 10^-2. That is a property of the sample size rather than of the effect sizes, and it is why the thresholds we use below are far less stringent than they would be for a cohort of hundreds of patients.

Over-representation analysis¶

We use the gseapy module to run an overrepresentation analysis. gseapy can run the test locally if we give it the gene sets themselves, but when we name one of the gene-set libraries of the Enrichr web server, as we do here, it sends our gene list and background to Enrichr, which performs the test.

In the analysis here we use the KEGG database’s definition of biological pathways. This choice can easily be changed to other databases such as GO.

As input we use the genes with $q$ values below 0.01. The background is the set of genes we actually measured, not all human genes – a pathway can only be over-represented relative to what could have been detected.

We keep only the $p$ values and term names, correct for multiple hypothesis testing, and list the terms with an FDR less than 20%. Enrichr only reports pathways that share at least one gene with our list, so the pathways with $p = 1$ are missing from the input to the correction; this biases the estimate of $\pi_0$ downwards and makes the $q$ values somewhat optimistic.

The analysis seems to find overrepresentation of relatively few pathways, given how convincing the differential expression was. The strongest term, Mineral absorption, is driven by the metallothioneins MT2A, MT1X and MT1E, which are indeed classically induced by glucocorticoids, together with the metalloreductases STEAP1 and STEAP2 – but it is a somewhat oblique way of describing a drug response, and nothing that we would call a glucocorticoid pathway comes out at all.

This is the characteristic failure mode of ORA. The 200-odd genes that pass the threshold are a thin slice of a response in which thousands of genes shift a little, and by cutting at $q<0.01$ we threw away all the information about how much each gene moved and in which direction.

Gene set enrichment analysis¶

We now use gseapy to run a gene set enrichment analysis (GSEA) on the whole ranked list instead of a thresholded one. This time we compare to the hallmark collection of MSigDB, fifty carefully curated sets that each summarize a coherent biological state.

We put the dexamethasone samples first so that a positive enrichment score means “up in dexamethasone”.

Why the permutation null matters here¶

Look closely at the nominal $p$ values: several of them cluster tightly just above 0.016, and none falls below it. That is not a coincidence. Permuting the phenotype labels of eight samples into two groups of four can produce at most $\binom{8}{4}=70$ distinct assignments, of which only 35 are distinct up to swapping the group names. The permutation $p$ value can therefore not be resolved beyond roughly $1/70$, no matter how many permutations we ask for, and the FDR estimates inherit that coarseness. Strictly, our design is paired, and permutations that respect the pairing can only swap the labels within each donor, which gives just $2^4 = 16$ assignments and an even coarser resolution.

With very few samples per group (the usual rule of thumb is fewer than seven), phenotype permutation runs out of resolution, and the recommended fallback is to permute the gene sets instead: the ranked list is kept fixed and random gene sets of the same size are drawn to build the null. This is a weaker null: unlike phenotype permutation, it treats genes as independent and so tends to be anti-conservative for sets of co-regulated genes, but it is the only one with enough resolution here.

We list the pathways with an FDR below 0.25, the threshold conventionally used for GSEA.

Where ORA found almost nothing, GSEA reports about twenty hallmark sets, all of them up in dexamethasone. Two of the strongest are exactly the positive controls we hoped for, although they are not labeled the way a pharmacologist would label them:

• Adipogenesis is a glucocorticoid signature in disguise. Dexamethasone is a standard ingredient of the cocktail used to drive adipocyte differentiation in culture, so the hallmark set is populated with genes that the glucocorticoid receptor induces directly.

• Androgen Response is a steroid-receptor set, and the nuclear receptors share a large part of their target repertoire; FKBP5 sits in its leading edge.

Reading the leading edge, and reading the label¶

The second hit, TNF-alpha Signaling via NF-kB, deserves a closer look, because at face value it says the opposite of what one expects from an anti-inflammatory drug. Inspecting its leading edge explains the discrepancy.

These are not markers of an activated inflammatory program. PER1, DUSP1 and KLF9 are well-known direct targets of the glucocorticoid receptor, ATF3 and NR4A1 are immediate-early genes that respond rapidly to many stimuli, and further down the list TNFAIP3 (A20) and NFKBIA ($I\kappa B\alpha$) are inhibitors of NF-$\kappa$B whose induction is one of the mechanisms by which glucocorticoids shut the pathway down. The gene set is named after the stimulus used to define it, not after the direction of regulation, and a set can be enriched because its brake was applied just as easily as because its accelerator was.

The cells here were never challenged with a cytokine, so there was no inflammatory tone to suppress; consistently, the Inflammatory Response hallmark itself is flat.

The enrichment score tells you that a coherent group of genes moved together; it takes the leading edge to tell you which genes did it, and hence what the hit means biologically.

Visualizing the enrichment¶

gseapy has built-in routines for plotting the enrichment of a pathway.

We can also compare the enrichment of several pathways in one plot.

A heatmap shows the normalized expression of the genes in a pathway, here the first 25 leading-edge genes of Adipogenesis, in each sample. The four columns on the left are the dexamethasone-treated cultures and the four on the right the untreated ones; the clean split of the leading-edge genes into two blocks is what the enrichment score is measuring.