What is dexamethasone doing to airway smooth muscle?
A complete Inflexa session on GSE52778: paired differential expression from a matrix that was never meant for it, a cross-check against the submitters' own pipeline, and a blind test that asks the signature to name the drug that produced it.
GSE52778 (Himes et al. 2014, PMID 24926665). Bulk RNA-seq of primary human airway smooth muscle cells from four donors, each treated with dexamethasone, albuterol, both, or left untreated for 18 hours. 16 samples, 4 supplementary files, 3.3 MB total. This tutorial uses the 8 dexamethasone and untreated samples, which form 4 donor pairs. GEO accession
This is the smallest useful mechanism-of-action analysis we know of. A glucocorticoid with a well-characterised receptor, applied to primary cells from four donors, with a published answer you can grade your own run against. It runs in under an hour on a laptop and the download is measured in megabytes.
It is also, quietly, an awkward dataset, which is the other reason to start here. What GEO hosts is an FPKM matrix rather than counts, so the tools most people reach for first do not apply. What happens next is the part worth watching: whether the analysis notices, says so, and picks a method that actually fits, or whether it runs DESeq2 on normalised values and reports a confident number.
Everything below is one session. The prompts are quoted exactly as they were typed, including the parts that were vague, and the screenshots are from the run that produced the figures on this page.
- The drug effect is larger than the donor effect. PC1 carries 47.3% of the variance and separates treated from untreated almost perfectly. For primary cells from four different people, that is unusual, and it is the first sign the experiment worked.
- Every canonical glucocorticoid target came back up and significant, none of them named in the prompt: FKBP5, KLF15, TSC22D3, DUSP1, PER1, SPARCL1, and CRISPLD2, the gene the original paper was built around.
- The cross-check agrees on direction and disagrees on how much is real. Fold changes correlate at Spearman 0.969, but the submitters' pipeline calls 295 genes significant where this one calls 4,126. That gap is left open rather than argued away.
- One claim was withdrawn on inspection. A tempting story about coordinated NF-κB repression did not survive a check against the actual fold changes, and the report says so instead of keeping it.
- The signature can partly name its own drug. Scored blind against 718,055 drug-perturbation experiments, dexamethasone's median instance sits at the 80th percentile. Two independent lines of evidence converge on the glucocorticoid receptor, one corroborates, and one turned out to be unanswerable with the reference data available.
The dataset, and what is awkward about it
Himes and colleagues treated airway smooth muscle cells from four donors with dexamethasone, a glucocorticoid, and albuterol, a beta-2 agonist. Those are the two mainstays of asthma therapy. The paper used the experiment to find glucocorticoid-responsive genes in these cells, and its headline result was CRISPLD2 as an asthma and glucocorticoid-response gene. That published answer is why this dataset is a good place to learn: you know what a correct run looks like before you start.
Four things about it are genuinely inconvenient, and all four matter:
- There are no counts. The GEO supplement is an FPKM matrix, 23,273 genes by 16 samples, produced by the Tuxedo pipeline against hg19. DESeq2 and edgeR model counts and cannot be pointed at this file honestly. The canonical raw counts exist in the Bioconductor
airwaypackage, and the reads are in SRA under SRP033351, but re-aligning turns a fifteen-minute exercise into an afternoon. - The design is paired and small. Four donors, one of each treatment per donor. Donor is the only axis of replication, so anything that ignores it is throwing away the structure that makes n=4 workable.
- Only one contrast is precomputed. The submitters included Cuffdiff results for dexamethasone versus untreated. The albuterol arms have no precomputed comparison, which is one reason to leave them out.
- ERCC spike-ins were added to only two of the four donors. Not fatal here, but it rules out spike-in normalisation across the whole set.
Before you start
This tutorial assumes Inflexa is installed, the sandbox is provisioned, and a model is configured. If any of that is not true yet, the getting started guide takes about fifteen minutes, most of which is the sandbox image downloading.
Requirements beyond that are small: Docker or Podman running, roughly 2 GB of disk for the dataset and outputs, and no special hardware. The second half of this tutorial downloads a 2.1 GB reference library, so leave room for that if you plan to follow it all the way through.
Make a directory, and start Inflexa in it:
mkdir GSE52778 && cd GSE52778 && inflexaEverything after this happens inside the terminal interface. No file was downloaded by hand and no script was written by hand.
1. Get the data in
You do not need to find the FTP path. Ask for the accession by name and let the agent resolve it.
download GSE52778 and add the files to the analysis inputs
Two things happen. The agent runs the CLI's own geo download subcommand, which pulls the four supplementary files, then registers all four as analysis inputs. Registering is what makes them visible to every later step and to the provenance record, and it triggers profiling, so by the time the download finishes the agent already knows the shape of what arrived.
Note the approval prompts. The first time a subcommand or a tool runs you are asked, and you can approve it once or for the rest of the analysis. Nothing runs unattended that you have not permitted.
2. Ask what it is before you ask what it means
This step costs one question and saves the analysis. Before deciding what to run, find out what you actually have.
What is this dataset about?
The answer comes from the profile of the files themselves plus the series metadata, not from the model's memory of the paper. It reports the four-by-four factorial design, the 18-hour treatment, the sequencing setup, and then a section headed "Caveats worth knowing before analysis" that lists exactly the four inconveniences described above: no raw counts, no precomputed albuterol contrasts, n=4 with donor as the only replication axis, and ERCC spike-ins in only two donors.
That last block is the reason this step exists. Nobody prompted for caveats. The awkwardness of the dataset was surfaced before a single analytical choice was made, which is the difference between choosing a method and defending one afterwards.
3. The prompt, and the plan it produced
Here is the whole scientific request. It states a question and two constraints, and deliberately does not state a method, a design, or a gene to look for.
I'm interested in what dex is actually doing in these cells. Differential expression first, then a pathway-level view of what responds. Ignore the albuterol arms. Start from FPKM, but compare with the provided Dex-vs-Untreated results to increase confidence.
What comes back is a plan, not a result. Four steps: QC and subset to the eight dexamethasone and untreated samples, paired donor-blocked differential expression, the Cuffdiff cross-check, then enrichment.
The paired design was proposed, not requested. Nothing in the prompt mentioned donors or blocking. It follows from the profile: with four donors and one sample per condition each, treating the samples as eight independent observations would spend the entire experiment estimating differences between people.
Then the part that matters most on this page. Before asking for approval, the plan flags two things about itself:
- limma on log2(FPKM+1), not DESeq2, with the reason stated plainly: no counts exist here. This is the fork in the road for this dataset, and it is resolved in the open, before anything runs.
- The cross-check is not independent replication. Cuffdiff processed the same eight samples through a different pipeline. Agreement demonstrates methodological consistency, and calling it validation would overstate it. The plan says so, and says it will degrade gracefully if that file's coverage turns out to be thin.
A tool that wanted to look good would have run DESeq2 on FPKM and reported four thousand significant genes. The interesting output of this step is not the plan. It is that the plan argues against its own strongest claim before you approve it.
4. Approve, and let it run
You can edit the plan, ask for a different design, or drop a step. Here it was right, so it was approved as written.
looks good, execute
The run detaches and reports back. Each step runs in the sandbox with no network access, writing figures, tables, and a provenance record as it goes. The sidebar tracks progress; you can keep working in the chat pane while it does.
This run took a few minutes. The whole session, including the perturbation analysis further down the page, took about an hour.
5. The gene-level answer
Start with the QC, because it decides whether anything after it is worth reading. PC1 carries 47.3% of the variance and separates dexamethasone from untreated almost perfectly: eta-squared for treatment on that axis is 0.993, for donor 0.005. Donor gets its own axis, PC2, at 16.3% and eta-squared 0.986. The two effects are close to orthogonal, which is the cleanest possible case for blocking.
Treatment beating donor as the top axis of variation is not what you normally see in primary cells from four unrelated people. It says the drug effect is strong and consistent. The per-gene picture is the complement of that: averaged across genes, donor accounts for more variance than treatment does (0.492 against 0.330), meaning donor is diffuse background across the whole transcriptome while dexamethasone is concentrated in a specific responsive set. Both facts together are what justify the paired design rather than merely permitting it.
The model is limma-trend with robust empirical-Bayes moderation, fitted as treatment plus donor on genes with FPKM above 1 in at least four of the eight samples. That leaves 12,094 genes, used both as the modelling matrix and, importantly, as the background universe for every enrichment test later.
It also leaves three residual degrees of freedom. Five design columns, eight samples. Empirical-Bayes moderation is the only reason inference is usable at that df, and it is worth holding onto that number, because nothing later in the analysis improves it. A reference library of 718,055 perturbation experiments does not add replicates to your experiment.
The result: 4,174 genes at FDR below 0.05, 2,369 up and 1,805 down, with 408 clearing an absolute log2 fold change of 1. And the part you can grade:
| Gene | logFC | FDR | What it is |
|---|---|---|---|
| SPARCL1 | +3.88 | 7.8e-05 | Matricellular protein, strongly glucocorticoid-induced |
| FKBP5 | +3.70 | 8.5e-04 | Co-chaperone in the GR complex, the textbook readout |
| KLF15 | +3.39 | 2.6e-03 | Transcription factor behind glucocorticoid metabolic effects |
| TSC22D3 (GILZ) | +3.01 | 3.0e-03 | Glucocorticoid-induced leucine zipper, direct GR target |
| DUSP1 | +2.90 | 1.9e-05 | MAPK phosphatase, the anti-inflammatory arm |
| PER1 | +2.60 | 4.3e-05 | Circadian gene, links GR to clock control |
| CRISPLD2 | +2.51 | 1.1e-03 | The gene the original paper was built around |
| CCL2 | -0.51 | 0.19 | Expected to fall, and it does not reach significance |
Seven canonical glucocorticoid targets, all up, all significant, none of them requested. CRISPLD2 at +2.51 is the strongest single piece of evidence that the pipeline behaves, because it is the result the original paper was built around and it was recovered without being asked for.
The eighth row is worth as much as the other seven. CCL2, a chemokine you would expect a steroid to suppress, trends down at -0.51 and misses significance at FDR 0.19. It is reported as a trend and not as a hit. The threshold was not loosened to make the anti-inflammatory story tidier, and the summary says so out loud.
6. The cross-check, and the gap it opened
The submitters' Cuffdiff table, restricted to rows flagged OK with finite fold changes, covers 13,370 testable genes, of which 11,877 join to the limma table. Two traps had to be handled before a single number could be compared, and both are the kind of thing that quietly corrupts an analysis:
- Cuffdiff reports the contrast the other way round, as log2(untreated / dex). Comparing without negating would have produced a perfect anticorrelation and a very confident wrong conclusion.
- 143 rows carry Cuffdiff's DBL_MAX sentinel for an infinite fold change. That value passes a naive finite-number check, so it survives the obvious filter and lands in the correlation as a real observation. They were removed by magnitude.
With those fixed, the two pipelines agree closely on effect: Spearman 0.969 on aligned fold changes, and 94.9% of joined genes moving in the same direction.
On how many genes are significant, they do not agree at all. In that shared universe limma calls 4,126 and Cuffdiff calls 295, roughly a fourteen-fold gap from the same eight samples.
The direction of the disagreement is informative rather than alarming. 289 of Cuffdiff's 295 sit inside limma's set, so the two are nested rather than contradictory: Cuffdiff is simply the more conservative caller at n=4. Those 289 genes, sign-concordant between pipelines by construction, became the high-confidence set that everything downstream was run on. 200 up, 89 down, led by VCAM1 (-2.83), MAOA (+3.11), TIMP4 (+2.25), ACSS1 (+1.64), and PRSS35 (-2.21).
The fourteen-fold gap could have been chased. We chose not to. Adjudicating limma's extra 3,800 calls needs either raw counts through a count-based model, which means the Bioconductor airway package or reprocessing SRA, or an independent replication cohort. Neither is a re-analysis of this input set, and neither fits in an hour.
So the report calls the gap unarbitrated rather than implying the extra calls are confirmed, and the biology downstream was run on the 289 genes both pipelines agree about instead of the larger and more flattering list. You might make that call differently. What matters is that the uncertainty is stated where a reader will see it.
7. Pathways, and how thin the certification is
Enrichment ran two ways against four collections, with the 12,094-gene expressed set as background throughout: over-representation on the high-confidence genes, and rank-based GSEA on the full moderated-t vector so no cutoff is needed.
GO output was collapsed with rrvgo before anything was reported. 74 significant GO biological process terms reduce to 20 semantic clusters, which is the difference between one finding and twenty near-duplicates of it dressed as twenty findings.
The top clusters on the up-regulated half are cellular response to hormone stimulus and endocrine process. That is the positive control landing: no glucocorticoid-specific term was sought, and hormone response is what came out on top. Below those sit regulation of glucose transmembrane transport, regulation of actin filament-based process, and regulation of cell adhesion, which is metabolic reprogramming and cytoskeletal remodelling, both expected of a glucocorticoid in a contractile mesenchymal cell.
GSEA is where the sample size shows. Across all four collections combined, four gene sets clear collection-wide FDR 0.05. Four. The one worth naming is Rho protein signal transduction (NES 2.151, adjusted p 0.0246, leading edge COL1A2, RAC1, RHOB, ITGB1), which in a contractile cell treated with a drug given for airway constriction is a better result than its p-value suggests. Three WikiPathways sets make up the rest.
Several biologically apposite terms sit just outside: glucose import at adjusted p 0.098, iron ion transport at 0.098, muscle hypertrophy at 0.112. With three residual degrees of freedom, correction across thousands of sets is punishing. The ranked signal is real and thinly certified, which is a statement about power rather than about biology.
Over-representation on the full 4,174-gene list instead of the 289-gene intersection produces a far longer result whose top collapsed cluster is "circulatory system process" at a score of 9.19. Bigger numbers, more terms, and a drift toward diffuse enrichment that says less about dexamethasone than the shorter list does.
It is the same data and the same code. The only change is which gene list you feed it, and that choice is worth more than any parameter in the enrichment call.
8. A claim that did not survive checking
This is the most useful thing on the page, and it started as a footnote.
In its own written summary, the enrichment step noted that some of its gene-level statements about NF-κB had been computed from the upstream differential expression table rather than from an output file belonging to that step. Nothing was obviously wrong with them. They simply were not backed by an artifact, and the step flagged that about itself rather than waiting to be caught. You can see the flag at the end of the run summary.
Asked to either produce the artifact behind those statements or strike them, the verification pass did neither. It went to the limma table, pulled every NF-κB core component, and withdrew the claim:
NFKBIA is up (+0.78). So are RELA (+0.21), CHUK (+0.33), IKBKB (+0.22), NFKB2, and NFKBIB. Only NFKB1 (-0.21) and RELB (-0.31) are down, and marginally. Every one of these effects is under 0.8 in absolute log2 fold change. NFKBIA rising is consistent with glucocorticoid transrepression on its own, since inducing the IκBα inhibitor is a classical mechanism, but a rising inhibitor is not a repression signature and the rest of the components do not line up behind it.
The report's conclusion, quoted exactly: "There is no coordinated NF-κB suppression signature in these data, and NF-κB core component directions must not be presented as supporting the anti-inflammatory finding."
What the anti-inflammatory story actually rests on, after that, is narrower and sturdier. VCAM1 is the strongest-FDR gene in the entire 289-gene set: down 2.83 in limma at FDR 1.4e-06, down 3.58 in Cuffdiff at q 2.0e-05. It is an adhesion molecule that recruits leukocytes into the airway wall and a canonical target of glucocorticoid transrepression. It then turns up again, unprompted, as the leading-edge gene of the repressed-by-glucocorticoid set in the blind test below. Two observations of the same repressed target is what the evidence supports, and a pathway-wide claim is not.
A verification pass that only ever confirms is not a verification pass. This one cost a paragraph of the story and made the rest of it worth more.
9. Closing the loop: can the signature name its own drug?
Everything so far confirms biology we already knew was there. The harder question, and the one that actually tests a mechanism-of-action workflow, is the reverse: hand the signature to something that has never been told what produced it, and see whether it points back at glucocorticoids.
Take the 289-gene set toward drug/target context. What I'd like to see is the loop closing: treat the signature as if the perturbagen were unknown and tell me whether it points back to glucocorticoids / GR connectivity against a signature library, GR regulon or ChIP-derived target overlap, whatever you can defend with what you have locally.
The response begins by refusing to plan. Whether it can defend an answer depends on which reference libraries are actually staged locally, so it checks first. Six queries against the reference store later, it reports what is installed, what is catalogued but not downloaded, and what each option would buy.
This is the behaviour that separates an analysis tool from a chat window. The answer to "can you do this" was not assumed and it was not guessed. It was looked up, and the honest position, in its own words, was that the single best tool for the question was catalogued but not installed.
What followed was a recommendation with trade-offs attached. Three catalogued datasets go at the question directly: the LINCS L1000 chemical perturbation library at 2.1 GB, MSigDB's published perturbation signatures at 2.9 MB, and GTRD ChIP-derived transcription-factor target sets at 4.1 MB. The advice was to take all three, because three unrelated priors that can each fail independently make a negative result interpretable instead of ambiguous.
One caveat was raised before anything was downloaded, which is the right time to raise it. LINCS is dominated by cancer cell lines rather than primary airway smooth muscle, and L1000 measures around 978 genes directly with the rest inferred. If dexamethasone ranks top anyway, that is strong. If it ranks mid-pack, cell-context mismatch is a live explanation rather than an automatic negative.
Remember that sentence. It decides how to read the main result before the main result exists.
let's download all 3
10. The blind test, and the rules it set for itself
The plan that came back is six steps: prepare the signature once, run four arms in parallel, then synthesise. The rules attached to it are what make it a test rather than a demonstration:
- No arm may filter its library to glucocorticoid terms before scoring. Every arm scores its full library and reports where the right answer landed. Rank out of library is the deliverable, not a p-value on a term chosen in advance. In the connectivity arm, compound and cell-line names were not even parsed until after the ranking was fixed, so parsing choices could not bias the scores.
- The two curated-knowledge arms do not count as two confirmations. Transcription-factor regulons and ChIP-derived target sets both encode existing literature about the glucocorticoid receptor. The synthesis was required to say so, and to flag that connectivity is the only arm reading observed drug responses.
- A missing reference is reported as unavailable, not as a null result. If GTRD has no NR3C1 set, that arm answers nothing, and the plan refuses to let it read as evidence of absence.
It also set an expectation in advance: a mid-pack dexamethasone rank would not be a clean negative, for the cell-line reason above. Committing to how a weak result will be read, before seeing it, is the difference between an analysis and a story.
11. What the four arms found
| Arm | Reference | Library | Result | Weight |
|---|---|---|---|---|
| Connectivity | LINCS L1000 | 718,055 drug-perturbation experiments | Dexamethasone's median instance ranks 143,716, the 80th percentile. Better than every one of 30 permuted nulls. | Independent, converges |
| Published perturbations | MSigDB C2:CGP | 3,540 signatures scored of 3,555 | The repressed-by-glucocorticoid set lands at the far negative tail, which is the concordant end. | Independent, converges |
| TF activity | CollecTRI, DoRothEA | 601 and 257 transcription factors | NR3C1 positive in both: rank 24 and rank 30, upper tier in each and top of neither. | Dependent, corroborates |
| ChIP targets | GTRD | 506 target sets | NR3C1 is not in the collection at all. No rank exists to report. | Unavailable |
Connectivity: real and moderate
This is the arm that carries weight, because it is the only one whose reference is an observed drug-response transcriptome rather than curated knowledge about the receptor.
The headline is the median, and the report is insistent about that. Across all 467 dexamethasone instances in LINCS, the median rank is 143,716 of 718,055, the 80th percentile. The best single instance ranks 5th, and the worst ranks 715,930. Leading with the best would inflate the apparent recovery by a wide margin, so it is reported second.
Against a permutation null matched on signature size, the observed mean rank of 215,516 falls below all 30 permuted means, giving an empirical p below 0.033. That is resolution-limited by the number of permutations rather than a precise p-value, and it is described that way.
So: dexamethasone sits meaningfully above chance in a library of nearly three-quarters of a million perturbation experiments, without being the top hit.
Then the detail that explains the median. Every top-ranked dexamethasone instance, ranks 5, 6, 9, 12, 14, 16, 19, 20, 22, 23, and 24, comes from the same fibroblast-lineage plate series, all at 6 hours and 0.37 micromolar. Airway smooth muscle is mesenchymal. The drug's signature matches best in the one lineage that shares the query's cell context.
This was identified from metadata after the ranking was fixed, which is what keeps it a finding rather than a rescue. It also implies the median understates what a context-matched screen would return, since most of the 467 instances come from cancer lines unrelated to the query tissue. The prediction made before the download is the effect that appeared.
Two structural limits stay attached to this arm regardless. No LINCS cell line is airway smooth muscle, so the fibroblast panel is the nearest lineage and not the same tissue. And 74 of the 289 signature genes were never observed anywhere in the L1000 gene space during scoring.
Published perturbations: the result that looks like a failure
All 3,555 published perturbation signatures in MSigDB's C2:CGP collection were scored by preranked GSEA, 3,540 of them computably. Seven of them name-match glucocorticoid, dexamethasone, or corticosteroid keywords, and the non-cherry-picked headline is the median of those seven: rank 2,169 of 3,555, the 39th percentile. That is the weakest number in the whole analysis and it is reported first.
The reason it is so weak is a lesson about ranked tables. The library is sorted by descending enrichment score, so rank 1 is the strongest positive concordance and rank 3,555 the strongest negative one, with indifference in the middle. Rank is a position on a signed axis, not a strength ordering, and a median over seven sets scattered across both ends of that axis means very little.
Which brings us to the set that matters: DEBOSSCHER_NFKB_TARGETS_REPRESSED_BY_GLUCOCORTICOIDS, ranked 3,522 of 3,555. Read carelessly that is thirty-third from the bottom, and it looks like the most on-the-nose glucocorticoid set in the collection failing outright.
It is the opposite. The set is named for genes glucocorticoids repress. Our query is a glucocorticoid response. For the set to agree, its members must go down, and a rank at the negative extreme is exactly that. Its enrichment score of -2.00 is the largest magnitude of any of the seven, and its nominal p of 0.0023 is the only one below 0.01. Its leading edge is VCAM1 and IL12A, and VCAM1 is the strongest-FDR gene in the entire high-confidence set.
The expected result, reached from the far end of the table, and one that anyone treating rank as a strength ordering would have written up as a miss.
One more thing this arm exposes about its own instrument. Seven sets is a lower bound, because name-matching only finds glucocorticoid biology that somebody put a drug name on. The clearest case is CHEN_LVAD_SUPPORT_OF_FAILING_HEART_UP, which ranks 2nd of 3,555 overall and whose leading edge is FKBP5, TSC22D3, KLF9, ZBTB16, PER1, and CEBPD, all canonical direct receptor targets. Nothing in its title suggests glucocorticoids. Searching set contents rather than set titles would be the fairer instrument, and would change this arm's denominator.
Transcription-factor activity: right direction, confounded significance
NR3C1, the glucocorticoid receptor itself, scores positive in both regulon resources: rank 24 of 601 in CollecTRI (score +3.81, adjusted p 0.0025) and rank 30 of 257 in DoRothEA (score +1.03, adjusted p 0.85). Upper tier in both, top of neither.
Both significance figures are confounded, in opposite directions. NR3C1's CollecTRI regulon holds 147 in-universe targets, in the 96th percentile of regulon sizes, and the standard errors of this method shrink as regulons grow, so part of that low p-value is a size advantage. In DoRothEA the same receptor has only five usable targets, sitting exactly at the inclusion floor, which is why it is directionally right and statistically powerless.
The two resources are also not fully independent of each other: some CollecTRI edges carry DoRothEA provenance. So the durable finding here is not either p-value. It is that two curated resources put the same receptor in the same direction near the top of a thousand candidates, and that this is corroboration of prior knowledge rather than a new line of evidence.
ChIP targets: no answer, which is not the same as a negative
NR3C1 is absent from all 506 sets in the installed GTRD collection. Exact-name and alias searches both returned zero matches, and no proxy transcription factor was substituted for it. There is no rank to report. All 506 sets were scored against the signature anyway, and none reaches FDR 0.05.
The report's guardrail on this is the sharpest sentence in it: the arm must never be reported as "GR showed no ChIP-derived target overlap", because that sentence asserts a measurement that was never made. A missing reference removes a confirmation slot. It does not supply a negative, and a scoring scheme that lets a missing file count against the hypothesis is broken.
The verdict
On the biology, the run matches the published result and adds to it. Every canonical glucocorticoid target was recovered without being named, CRISPLD2 included. The paired design the data profile implied turned out to be the right call, and the analysis showed its work on why. The programmes that came out of the high-confidence set, transrepression at VCAM1, matrix and cytoskeletal remodelling through Rho signalling, and metabolic reprogramming, are what a pharmacologist would expect a glucocorticoid to do in a contractile airway cell.
On the blind test, the report's own summary is the right one to quote: "Net: two independent arms converge, one dependent arm corroborates, one arm is unavailable. PARTIAL SUCCESS. No arm's non-cherry-picked headline metric places GR content at the literal top of its library."
That is a more useful sentence than a clean win would have been. Connectivity and published perturbations converge from genuinely independent starting points. Transcription-factor activity agrees but is reading back curated knowledge about the same receptor, so it does not get counted twice. The ChIP arm could not run. Calling that a four-way confirmation would be wrong, and the run does not.
What remains uncertain is stated rather than smoothed. The fourteen-fold gap between the two differential-expression methods is unarbitrated. The enrichment result rests on four gene sets clearing correction. The connectivity median is mid-pack, for a reason predicted in advance that is still a limitation. And the number that bounds all of it does not move: four donor pairs, three residual degrees of freedom, no matter how large the reference libraries downstream are.
If you run this yourself, those are the numbers to compare against. A run that recovers CRISPLD2 and FKBP5, and puts dexamethasone above chance in a connectivity library, is working. A run that reports six hundred significant pathways and no caveats is telling you something about the tool rather than about the biology.
What is still open
The report closes with three questions it cannot settle from this input, each with a decision attached about whether it is worth pursuing. They are the honest to-do list for anyone taking this further:
- Adjudicating the limma and Cuffdiff gap. This needs raw counts through a count-based model, from the Bioconductor
airwaypackage or by reprocessing SRA, or an external replication cohort. It is not a re-analysis of these files. - Whether the mesenchymal cell-context effect is general. Do the fibroblast-lineage plates dominate the top ranks across all 28 corticosteroids, or only for dexamethasone? The per-instance metadata to test it is already in hand, and the answer distinguishes a lineage effect from a compound-specific one.
- Replacing name-matching with content-matching. Scoring the perturbation collection by canonical receptor-target membership rather than by set title would change that arm's denominator, and possibly its headline. CHEN_LVAD at rank 2 is the proof that the current instrument undercounts.
None of them were run. The session was scoped to about an hour and stopped where the evidence stopped, which is also why the write-up can be specific about what it does not know.
The full report, and the data package
This page is the walkthrough. The full report and the data package carry the rest, both in a single Zenodo deposit, and they exist so that you do not have to take the write-up on trust.
The generated report
"Dexamethasone in Primary Human Airway Smooth Muscle (GSE52778): The Transcriptional Response, and a Bounded Blind Recovery of Glucocorticoid Biology". Every figure, every table with its data downloadable, the methods as executed, the evidence accounting, and the caveats in full rather than in summary. It is the document this page was written from.
The report ships as a self-contained HTML file inside the Zenodo package below. Open it in a browser after downloading.
The data package
Inputs, the scripts as they were generated and run, and every output file, archived on Zenodo with a DOI. Enough to reproduce the run without this page, or to check any number on it.
10.5281/zenodo.22081987Your turn
The interesting version of this is not repeating it. It is aiming the same question at a compound you actually care about. The prompt barely changes:
I'm interested in what [your compound] is actually doing in these cells. Differential expression first, then a pathway-level view of what responds. [Any arms to ignore.] [Where to start from, and anything to cross-check against.]
Three habits from this run transfer to any dataset.
Ask what the data is before asking what it means, because the caveats that come back constrain every method choice after them. When a result confirms what you expected, make something try to disconfirm it: the blind arms here were more informative than the differential expression precisely because they could have failed. And when a claim cannot be traced to an output file, treat that as a reason to check it rather than a formality, because on this run that check removed a paragraph that would otherwise have been published.
If you want to start with this exact dataset, it is four files and 3.3 MB, and the first prompt is on this page.
Run this on your own data.
Inflexa is free and open source under Apache 2.0. It runs on your machine, with the model you choose, and every result carries a record of how it was made.