Dexamethasone response in A549 — bulk RNA-seq differential expression and pathway analysis

Study: GSE229084 — Transcriptomic profiling of A549 cells in response to Dexamethasone, in the context of GR inhibition based on antagonists as well as a novel GR-PROTAC (PMID 38071198). Human A549 lung adenocarcinoma cells (ATCC CCL-185), 51 RNA-seq libraries, GRCh38. Contrast analysed: dexamethasone 100 nM, 18 h vs untreated — 4 vs 4 replicates. Software: PyDESeq2 0.5.4 (the DESeq2 negative-binomial GLM in Python), gseapy 1.3.0, MSigDB 2024.1.Hs.

1. Headline

Dexamethasone produces a large, well-replicated transcriptional response in A549: 5,296 genes at FDR < 0.05, of which 898 pass a |log2FC| > 1 threshold (400 up, 498 down). The signature is the canonical glucocorticoid-receptor response — FKBP5, TSC22D3 (GILZ), ZBTB16, ANGPTL4, CIDEC, MT2A, PER1, DUSP1 — together with coordinate repression of proliferation genes.

At the pathway level the ranked list splits into two coherent blocks: down, the cell-cycle / DNA-replication machinery; up, mitochondrial oxidative phosphorylation and a lipogenic–adipogenic metabolic programme.

The strength of the evidence is very different for those two layers, and that difference is the most important thing in this report. The gene-level result is solid. The individual pathway statistics are not — the set-level FDRs from preranked GSEA do not survive a permutation control. Section 6 gives the numbers.

2. Data and design

2.1 What the series contains

GSE229084 has 51 libraries in 13 treatment arms (~4 replicates each), of which the ones relevant here are:

The remaining ten arms are GR antagonists (mifepristone, CORT113176) and the GR degrader KH-103 with or without dexamethasone; they are outside this analysis.

Two design features matter:

• Replicate and batch coincide. Batches set1–set4 correspond to replicate 1–4 (from !Sample_characteristics_ch1), so the batch effect is the replicate effect. Replicate is therefore entered as a block, giving the paired design ~ replicate + group. There is no separate batch factor to correct for, and none is estimable.
• The arms are balanced across blocks. Each of replicates 1–4 contributes one sample per arm (18 h DMSO is missing its replicate-2 sample), so the block × treatment structure is close to complete and the block term is estimable.

2.2 Counts — one caveat to record

The supplementary file GSE229084_geneCounts.tab.gz is described by the submitters as "a tab-delimited matrix of estimated read counts per gene and sample". The values are not integers: ~21% of entries carry a fractional part.

I checked whether that fraction is a normalisation before using the file:

• maximum deviation from the nearest integer is 0.500 — i.e. every value is a whole count plus sub-unit noise;
• the median of value / round(value) over values > 50 is exactly 1.000000 in all 51 columns, so there is no per-sample scaling.

The jitter is additive and sub-unit, not multiplicative, so rounding recovers the integer counts essentially exactly. All DESeq2 input was round()ed. This is a deviation from a strict "raw integer counts" reading of the file and is flagged here rather than buried.

2.3 QC

• Library sizes 17.5M–24.9M reads; 14,760 genes have ≥ 10 counts in at least 3 samples.
• PCA on log2(CPM+1): PC1 (59.2%) separates DEX from both control arms; PC2 (11.1%) tracks the replicate/batch axis. Untreated and vehicle overlap — see results/qc_pca.png.
• Sex-marker check is uninformative here (single male cell line, XIST absent, RPS4Y1 detected), so no sample-swap signal is available from it.

3. Method

1. Model — PyDESeq2 0.5.4, design ~ replicate + group over the 11 samples of the three arms above, median-of-ratios size factors, parametric dispersion fit, Wald test, Benjamini–Hochberg across all tested genes.
2. Shrinkage — empirical-Bayes LFC shrinkage (PyDESeq2 lfc_shrink) is used for ranking and for the volcano; unshrunk LFC and its standard error are what the ranked table and the marker-gene numbers below report, since shrinkage biases low-count genes toward zero.
3. Three contrasts from the one model — dexamethasone vs untreated (primary), dexamethasone vs vehicle (corroboration), vehicle vs untreated (empirical null). Each is its own BH family.
4. Enrichment — preranked GSEA (gseapy, 1,000 gene-set permutations, seed 0) on the Wald statistic over the 16,069 genes that passed independent filtering, against MSigDB 2024.1.Hs Hallmark (50 sets) and Reactome (1,008 sets after the 15–500 size filter); plus over-representation on the 898-gene thresholded list with the 16,069 tested genes as background.
5. Redundancy collapse — significant sets clustered by Jaccard > 0.5 on leading-edge genes; 172 significant sets collapse to 93 clusters.
6. Negative control — the entire pipeline (DE → ranking → GSEA) re-run on treatment labels permuted within each replicate block, which is the correct randomisation for this blocked design.

4. Gene-level results

Of 16,069 genes tested:

Top of the ranked list (unshrunk LFC ± SE):

FKBP5, TSC22D3, ZBTB16, PER1, DUSP1, PDK4, ANGPTL4, CIDEC and SGK1 are all established direct GR targets; their joint induction at these effect sizes is the expected positive control for the assay, and it is reproduced here.

4.1 Two independent control arms agree

If the DEX effect were an artefact of the design, the two control arms would behave differently. They do not:

• Vehicle vs untreated: 1 gene of 22,988 at FDR < 0.05. 0.1% DMSO is transcriptionally inert here, which is the empirical null for this dataset.
• DEX vs vehicle: 4,726 genes at FDR < 0.05 — the dexamethasone signature reproduces against a second, independent control.
• Of the 848 genes significant (FDR < 0.05, |log2FC| > 1) in both treatment contrasts, 100% move in the same direction, and the two contrasts' log2FC correlate at r = 0.99997.

5. Pathway results — which pathways move

GSEA over Hallmark + Reactome was run on the Wald statistic of all 16,069 tested genes. 172 of 1,058 sets reach FDR < 0.05 (20 up, 152 down). The direction split is dominated by the fact that Reactome devotes many redundant sub-terms to cell-cycle and DNA-replication processes.

Hallmark, the readable summary — 10 of 50 sets at FDR < 0.05. The last three columns are the permutation-control verdict (§6.2): the null mean and maximum of |NES| for that term over 19 label-permuted runs, and how many of those runs matched or exceeded the observed value.

Only four of these ten — ADIPOGENESIS, SPERMATOGENESIS, HYPOXIA, ANDROGEN_RESPONSE — exceed the observed |NES| of every clean null run, and with 19 permutations the finest tail probability available is p = 0.05. The top-ranked term, E2F_TARGETS, is matched or beaten by a permuted ranking six times out of nineteen.

5.1 Reading the leading edge, not the label

The set names are misleading in two places, and the leading-edge genes say so:

• HALLMARK_OXIDATIVE_PHOSPHORYLATION (127 leading-edge genes) is the real metabolic result: PDK4, NDUFA8, NDUFA5, ATP5F1B, ATP5MC3, ATP5PB, CYC1, UQCRH, SDHD, ACO2, LDHB. PDK4 is a direct GR target that shuts off pyruvate dehydrogenase and pushes metabolism toward fatty-acid oxidation; the accompanying induction of respiratory-chain subunits and of the lipogenic genes in the adipogenesis set (ANGPTL4, CAVIN2) is a coherent GR-driven metabolic remodelling.
• HALLMARK_HYPOXIA (53 leading-edge genes) is not hypoxia. Its leading edge is ANGPTL4, MT2A, CAV1, TGM2, DUSP1, DDIT4, SERPINE1, ZFP36 — GR target genes that happen to be members of the hypoxia set. The same warnings apply to the ORA result below; enrichment of a set is not evidence that the process is active.

5.2 Over-representation on the 898-gene list

With the 16,069 tested genes as background: 20 sets up, 46 down at FDR < 0.05.

• Up, most specific: REACTOME_METALLOTHIONEINS_BIND_METALS 5/5 and REACTOME_RESPONSE_TO_METAL_IONS 5/8 — the complete metallothionein family (MT1A, MT1B, MT1F, MT1X, MT2A). A five-for-five, family-complete hit is the cleanest single pathway result in this analysis. HALLMARK_HYPOXIA 18/182 and HALLMARK_TNFA_SIGNALING_VIA_NFKB 19/181 follow; the NF-κB set's leading edge is NFKBIA, SGK1, PER1, DUSP1, KLF9, ZFP36, TNFAIP3, G0S2 — again mostly canonical GR targets, not evidence of NF-κB activation.
• Down: HALLMARK_G2M_CHECKPOINT 42/199, HALLMARK_E2F_TARGETS 38/199, REACTOME_CELL_CYCLE 67/652, REACTOME_CELL_CYCLE_MITOTIC 57/527. The leading edge is TOP2A, CDK1, BIRC5, AURKB, PLK1, MKI67, PCNA, MCM5, MCM7 — the glucocorticoid anti-proliferative response.

5.3 Two of the significant sets are proliferation programmes in disguise

HALLMARK_SPERMATOGENESIS (leading edge: CDK1, AURKA, CCNB2, NEK2, BUB1, EZH2, TTK) and HALLMARK_ANDROGEN_RESPONSE are not germ-cell or androgen biology in a lung adenocarcinoma line; they are well-known proliferation-correlated gene sets. They should be read as further cell-cycle signal, not as their labels.

6. How confident should we be

6.1 Gene level — high confidence

Permuting treatment labels within each replicate block and re-running the whole DE pipeline, 24 times:

Every draw that behaves as a null returns 0–5 genes; BH at FDR < 0.05 returns essentially nothing when there is nothing, which is exactly what it should do. The single outlier (5,118) is not a calibration failure — it is a signal-preserving relabelling. Because each block contains exactly one sample per arm, some within-block permutations simply swap the group names and reproduce the treatment contrast (that run put all four DEX samples into one group; its log2FC profile correlates with the real one at −0.60, i.e. the same effect with the sign flipped). It is excluded from the null in §6.2 for that reason, and it is worth stating plainly: this permutation scheme has a non-trivial heavy tail, so a small number of permutations is not enough to characterise it.

The gene-level result is also supported independently by the two control arms agreeing on 848/848 genes in direction (§4.1) and by the recovery of the canonical GR target set.

6.2 Pathway level — limited confidence, and the FDRs should not be quoted as they stand

The same permutation exercise applied to preranked GSEA is a warning:

The clean null runs produce more "significant" Hallmark sets than the real analysis does. Preranked GSEA's gene-set permutation does not account for gene–gene correlation, so its FDR is badly anti-conservative here: it cannot be used as the confidence measure, and the FDR column in §5 is a ranking device only.

What I can defend instead is the per-term empirical tail probability against the 19 clean null runs (the runs ≥ observed column of §5). It is only partly encouraging:

• Above every null run — ADIPOGENESIS, SPERMATOGENESIS, HYPOXIA, ANDROGEN_RESPONSE. With 19 permutations the floor is p = 0.05; no term can be shown to be better than that.
• Not distinguishable from the null — E2F_TARGETS (6 of 19), DNA_REPAIR (4 of 19), OXPHOS and MITOTIC_SPINDLE and HEDGEHOG (3 of 19 each), G2M_CHECKPOINT (2 of 19).

The reason is visible in the null means. E2F_TARGETS has a null |NES| averaging 2.43 against an observed 3.12: it is a large, highly co-expressed module whose rank statistic is intrinsically noisy, and a permuted ranking reaches an equally extreme value often. A pooled null — pooling all 50 terms × 19 runs, which treats terms as exchangeable — would instead give E2F p = 0.006 and OXPHOS p = 0.021. I do not report that as the headline result, because the per-term nulls plainly are not exchangeable (null mean 2.43 for E2F vs 1.24 for ANDROGEN_RESPONSE) and pooling would import exactly the assumption the per-term test exposes as false.

So the two most prominent pathway hits by GSEA FDR — E2F and G2M — are the two that the control removes. The cell-cycle biology is not in doubt, because it rests on hundreds of individually significant genes (TOP2A, CDK1, MKI67, PLK1, PCNA, MCM5/7 all at FDR < 1e-38); what is in doubt is the claim that the E2F gene set is non-redundantly enriched beyond what a permuted ranking produces.

Bottom line on pathways: treat §5 as a ranked description of what the dexamethasone-regulated genes are, not as a set of established enrichments. The gene lists are trustworthy; the set-level p-values are not.

6.3 Caveats to carry forward

1. Non-integer counts — rounded before modelling (§2.2). Sub-unit jitter, so the loss is negligible, but it is a deviation from strict raw counts.
2. n = 4 per arm. Dispersion is estimated from 5 residual degrees of freedom; the model is adequate but not generous, and the effect sizes that matter (FKBP5, TSC22D3, TOP2A) are large enough that this is not the limiting factor.
3. One cell line, one time point, one dose. Nothing here generalises to other lines, other glucocorticoids or other exposure times without replication.
4. Two replicate-numbering conventions collide. The 18h_DMSO columns are numbered 1/3/4 in the supplementary counts file but 1/2/3 in the series matrix. The block labels for those three samples are taken from the counts file. Because the vehicle arm is used only for the corroborating contrast, this cannot affect the primary result.
5. No RNA quality metrics were available for the 11 samples used (RIN, mapping rate), so sample-level technical quality beyond library size and PCA was not assessed.
6. The permutation control itself needed care. Permuting labels within replicate blocks is the correct randomisation for this design, but because every block holds exactly one sample per arm, a minority of permutations relabel the groups rather than randomise them — 1 of 24 here. Those runs must be excluded before the null is read (they are not nulls), and the exclusion rule is stated rather than hidden: a run is kept as a null if it returns at most 25 genes at FDR < 0.05. A different threshold would not change the conclusion, since the two populations are separated by three orders of magnitude (0–5 vs 5,118).

7. Files

Reproduction, in order: refit_unshrunk.py (model, Wald LFC/SE for all three contrasts) → make_figures.py (figures + ranked table) → perm_control.py and perm_sets.py/perm_sets2.py (24 within-block permutations) → diag22.py (diagnosis of the signal-preserving permutation).