What is IFN-β doing to PBMCs, cell type by cell type?

Two and a half minutes later, nine files and 74.4 MB are on disk, itemised by what each one is: the 72.7 MB tar of per-sample matrices, the gene tables, and the per-batch cell metadata carrying the t-SNE coordinates and the author cell-type labels.

Note what it does not do. The files are downloaded but not registered, and it says so and asks. Registering is what makes an input visible to every later step and to the provenance record, so it is a decision, not a formality.

Eleven seconds. All nine paths registered, and profiling starts in the background, so by the time you have a question about the design the answer already exists.

The last sentence is the one that does the work. It sets a standard and leaves the method open. What comes back is a five-step plan on two tracks, generated from the data profile rather than from a template: build the batch-2 singlet object and the paired design table, run paired pseudobulk DE per lineage, then split into a driver-gene track and an interaction track that feeds the pathway step.

Approved as written, the run detaches under identifier 7a82f521-10f9-4bd8-8086-d7ac6e4867a4 and completes five steps in 2 hours 6 minutes. Each step runs in the sandbox with no network access, writing figures, tables, logs, scripts and a provenance record as it goes.

What comes back when it finishes is the whole result in one screen. Every number in it is unpacked in the sections below, and every one of them is traceable to a file in the run directory. The summary opens by reading one of those files back.

4,386 droplets are dropped as doublets or ambiguous, 6 more carry no cell-type label, and 111 fail QC. The two arms end 40 cells apart, so the paired contrast is not being driven by differential recovery.

Two things in this step are worth more than the table. The first is the barcode join. The batch-2 metadata TSV has seven header names and eight fields per row, and 313 stimulated barcodes carry an R make.unique() mangling that would have made them fail to match the matrix. After correction the join is exact: 14,619 of 14,619 and 14,446 of 14,446, both arms at 100%, with a log line confirming no cell was lost and therefore no silent merge occurred.

The second is the mitochondrial filter, and it is the best example on this page of a check that earns its keep by failing.

From there the design is enumerated rather than assumed. The unit is donor × cell type × condition, giving a complete grid of 128. A unit is testable at 10 cells or more; a cell type qualifies only when at least 6 of the 8 donors have both arms testable. Seven of the eight types qualify. Megakaryocytes do not: 208 cells, only 4 donors paired. The summary names them and says which threshold they missed, so they are excluded rather than quietly absent.

Dendritic cells qualify at exactly the boundary: 6 of 8 donors, because donor 1039 has 9 control DCs and donor 107 has 8, both under the ten-cell floor. That number follows them through every remaining step on this page, and every step says so.

All eight donors (101, 107, 1015, 1016, 1039, 1244, 1256, 1488) appear in both arms, with none present in only one. The stimulation contrast is fully paired within donor, which is the property the whole design rests on. Nineteen verification checks confirm the object: shape, singlets only, batch-2 GSMs only, zero nulls in donor, condition and cell type, integer non-negative counts identical to the matrix, and the design table's cell counts and UMI totals re-derived from the object and matched.

That ordering matters for everything that follows. Lineage identity, not the stimulus, is the dominant axis of variation in this experiment. The paired design is what makes the smaller axis, the one you actually care about, estimable at all.

The ranking statistic is a fixed 13-gene ISG panel (ISG15, IFI6, IFIT1, IFIT2, IFIT3, MX1, ISG20, OAS1, STAT1, IRF7, LY6E, RSAD2, CXCL10), scored per donor as the mean stimulated-minus-control difference in TMM-normalised log2 CPM, then tested with a one-sample t on the paired differences. Because it is a within-donor CPM ratio it is independent of library depth, and therefore of cell number.

Read the fourth and fifth columns against the second. Dendritic cells sit second on effect size from 384 cells, while FCGR3A+ monocytes sit fifth despite having the second-highest DE-gene count. Ranking on counts would have swapped them. Effect size and DE count agree at only Spearman 0.714, and DE count is not even monotone in cell number, at Spearman 0.214, because both the number of cells and the intrinsic response feed it. That is what makes it unusable as a magnitude scale, and it is the reason the plan ruled it out before the numbers existed.

Then the check that decides how much of the ranking survives. Every analysed unit was downsampled to 15 cells per donor per condition, the size of the smallest analysed unit, and the entire pipeline re-run on three independent seeds.

The top of the table holds and the middle does not. Top-1 and top-2 are identical; top-3 is not; maximum rank displacement is three positions, with CD4 T cells falling from 3rd to 6th. One structure actually gets cleaner at equal cell numbers: all three myeloid lineages sit above all four lymphoid lineages after downsampling, which is not true in the full data. Seed-to-seed variation is small next to the depth effect, so the instability is caused by 15-cell libraries rather than by sampling noise.

A limma-voom refit on all seven lineages, same design and same normalised object, agrees on effect almost exactly, at Spearman 0.987 to 0.999 on fold changes, and disagrees systematically on significance. voom calls more genes in every single lineage: 96.7% to 99.5% of edgeR's calls are recovered by voom, but only 80.9% to 98.0% of voom's are recovered by edgeR. The direction of that gap is reported rather than averaged away, with edgeR-QLF named as the more conservative caller and kept as primary.

CCL8 tops both monocyte subsets, CXCL10 tops all three of CD4 T, CD8 T and NK, CXCL11 tops B cells, and IFIT2 tops dendritic cells. Six genes are top-25 drivers in all seven lineages: CXCL10, IFIT1, IFIT2, IFIT3, RSAD2 and ISG15. That is the first hint of what the interaction test will later make formal, which is that the lineages share a repertoire and differ in gain.

Two of those confidence intervals read "degenerate", and that is the run's own flag rather than a rendering failure. Genes that are near-absent in control produce quasi-complete separation, which leaves the deviance-based FDR valid and the Wald interval meaningless. 17 of the 291 driver rows carry the flag. They are marked in the tables rather than dropped or shipped silently, and ranking by absolute fold change is noted as being systematically enriched for exactly this failure mode.

The consistency check is the part that would be hard to fake. 285 of the 291 driver rows move in the model's direction in every single donor. The largest single donor's share of the summed per-donor effect has a median of 0.147 to 0.193 per lineage, against 0.125 for a perfectly even eight-donor split. The leading drivers are carried by all the donors at roughly equal weight, not by one enthusiastic outlier.

Exactly one driver is disqualified. CXCL13 in FCGR3A+ monocytes has an FDR of 7.0e-04 and a log2 fold change of 5.03, and would have read as a finding. But only 5 of 8 donors move in the model's direction, donor 1488 supplies 55% of the summed effect, and the leave-one-donor-out mean collapses to 0.42. It trips two independent flags, and the summary states plainly that it should not be reported as an IFN-β driver.

The cross-lineage heatmap is where a reader is most likely to over-read, so the figure argues against itself in its own title:

The step gives two concrete reasons the eye will mislead here. Missingness is structural: CXCL11 is untested in CD8 T cells, LAMP3 in CD8 T and NK, and none of those cells carries a zero. And topping one lineage's list is not specificity. Of the 25 genes that appear in exactly one lineage's top-25, 17 are not even tested in all seven lineages and 22 are significant with a fold change above 1 in more than one. ZBP1 tops only the CD14+ monocyte list yet is significant in all six lineages where it is tested. Deciding specificity, the step says, requires the interaction test in the next step.

The interaction outnumbers the uniform core by roughly thirty to one. The core genes are the ones you would nominate for a compartment-agnostic readout: MX1, IFI44L, IFI27, ETV7, IRF7, STAT2, LGALS3BP, SPATS2L, PARP10, PARP11, DTX3L, MOV10, C19orf66. The genes that are not in it are the surprise: CXCL10, IFIT2, TNFSF10, OASL, CCL8, CXCL11, SOCS1, CD38 and IFITM2 all classify as lineage-specific. They are induced everywhere. They are not induced equally, and the test can tell.

The interaction has a direction, and it is one axis:

Every myeloid lineage is dominated by stronger-than-average contrasts and every lymphoid lineage by weaker. This reproduces the lineage ranking from step 5 using a completely different statistic. There it was a per-donor ISG panel effect size; here it is a coefficient-equality test that never sees the panel. Two independent instruments landing on the same axis is what licenses the biological claim.

Three sensitivity analyses sit behind those numbers, and one of them moves them. Donor handling does not: refitting with donor as a random block instead of a fixed factor gives Jaccard 0.89 to 0.93 on the interaction set, and the fitted intra-donor correlation with donor removed from the design is only 0.048, so a random donor effect would recover almost nothing.

The choice of F statistic does move them. limma's moderated omnibus F rotates moderated t-statistics by one shared correlation matrix, while voom weights make the true coefficient correlation gene-specific. Against an exact per-gene Wald F the two agree at Spearman 0.976, but 274 genes, 6.40% of the total, change class (core 61 versus 68, lineage-specific 1,906 versus 1,846), and the switchers sit at the decision boundary, with a median limma FDR of 0.0437 against an exact FDR of 0.0602. Both statistics ship per gene, limma's is named primary and why, and the summary instructs that genes near the boundary are provisional and the counts should not be quoted without the range.

One more check earns its place. A lineage can look like a weak responder because it started high, so every weaker-than-average call was tested against its own control baseline:

SSB encodes La, one of the classical lupus and Sjögren autoantigens. The donors here are lupus patients, and the compartment where the type I interferon signature is measured clinically is the monocyte compartment. So a result showing that IFN-β induces the autoantigen specifically in myeloid cells is exactly the shape of a mechanism, which is why the more interesting thing is what the run does with it.

Three limits are attached to the claim at the point of making it, not in a footnote. The only prior mechanistic precedent for IFN-β induction of SSB is in salivary gland epithelium, a non-haematopoietic system, at around two-fold and with no protein change. A 2.3 to 3.3 log2 unit effect in human PBMCs far exceeds that, and wants protein-level confirmation. This is mRNA at six hours in cryopreserved cells. And because no donor clinical covariates ship with the data, the hypothesis cannot be tested against autoantibody titre or disease activity anywhere inside this run.

The word the run uses for it is "hypothesis". Not "finding", not "novel target". The report's own section heading is SSB/La: the sharpest lineage split, and a hypothesis rather than a result, and that is the correct altitude for a six-hour ex-vivo stimulation with one epithelial precedent behind it.

The interferon-alpha and interferon-gamma sets are the top two by enrichment score in all eight rankings, spanning just 2.29 to 2.76. Myeloid lineages suppress oxidative phosphorylation and MYC targets while mounting that response: CD14+ monocytes at NES -2.111 for OXPHOS, FCGR3A+ monocytes at -2.147. B cells' second-largest GO cluster is cytoplasmic translation, in the down direction. What individuates a lineage is what it gives up to respond, not which antiviral effectors it deploys.

The over-representation test makes the same point from the opposite direction, and this is the single most informative result in the step. The interferon sets are top-enriched in both the myeloid stronger-than-average foregrounds and the lymphoid weaker-than-average ones: FCGR3A+ monocytes' weaker set at odds ratio 55.2, NK cells' weaker set at 16.3, CD14+ monocytes' stronger set at 6.3. ISGs dominating both directions is the signature of one shared programme running at different gains.

Two pieces of hygiene are worth copying. GO biological process is a nested hierarchy, so the raw significant counts are not counts of processes: 1,570 raw hits collapse to 205 distinct ones under semantic clustering, an inflation factor of 5.8 to 8.9 across rankings. The step reports the collapsed numbers and says the raw ones overstate. And the enrichment background is the 4,274 symbols actually tested in the joint model, not the genome, with the 0.23% shrinkage from unmapped identifiers quantified rather than waved at.

One more piece of housekeeping the step insists on: "weaker" means below the seven-lineage average, not downregulated. All seven lineages induce ISGs strongly. The high interferon odds ratios in the weaker foregrounds do not mean ISGs fall anywhere.

GSVA agrees with the model-based ranking at Spearman 0.857, recovering both the top-1 and the top-2. singscore, on the identical matrix and the identical gene set, agrees at only 0.643 and disagrees on both. Two rank-based scoring methods, same input, different answer.

The step diagnoses this rather than picking the one it likes. singscore's paired difference correlates negatively with the control baseline (Spearman -0.464); GSVA's does not (+0.107). CD14+ monocytes start from a high control baseline with the second-lowest headroom on a bounded scale, and they nonetheless reach the highest stimulated score of any lineage, 0.3958 against a range of 0.2400 to 0.3769. They end furthest along and travel a shorter distance, and singscore reports the distance.

The instruction the step leaves behind is specific: do not report the singscore rank order as an independent lineage ranking, because it is scale-confounded. Report GSVA's, or better, report only the CD14+ monocyte and dendritic cell top-2 claim that all three statistics agree on.

The same bounded-scale logic recurs one section earlier. In the interaction model, 63% to 68% of the myeloid weaker-than-average calls occur where that monocyte subset already holds the highest control baseline. It is the same ceiling appearing in two unrelated statistics, and noticing that is what turns a methodological annoyance into a finding.

Three habits from this run transfer to any single-cell dataset.

State the pairing in the prompt, because it is the one fact about your design that the files cannot tell you and it changes every model downstream. Ask for a standard rather than a method. "Defend to a reviewer" is what produced the ban on per-cell tests, the interaction test and the depth-independent ranking, and none of those three were named. And when two ways of measuring the same thing disagree, treat the disagreement as the result: the singscore and GSVA split on this run was worth more than either ranking on its own.

If you want to start with this exact dataset, it is nine files and 74.4 MB, and the first prompt is on this page.