cyRAVEN 1.0.0
Explore mode, the statistics catalogue and the input format, together with parametric tests beside the rank tests, a feasibility check before either, the channel control that stops a cohort splitting into panels, per-marker UMAPs, parallel reading, and one fix for a failure that was silent and severe.
Nothing an earlier pre-release produced is renamed or removed, with the single deliberate exception in “A sample that fails staining QC is no longer assumed to be a control” below. That one changes numbers, because the numbers it changes were wrong.
Clinical variables are analysed instead of only carried
read_samplesheet() recognises a fixed list of subject attributes, patient_id, sex, age_years, height_cm, weight_kg, infection_focus, cohort, wbc_per_ul, and carried everything else as a “study column”, usable only to group samples. A SOFA score cannot group anything and a 28-day survival flag grouped as a cohort loses the fact that it is an outcome, so those columns were read, written into the manifest, and never analysed.
-
--clinical-columns(orsamples: clinical_columns:in the config) names any sheet column as a clinical variable in its own right. Both metadata tables are searched, because the sheet’s own columns end up split between them. - The test follows the variable’s type: Spearman for numeric, Wilcoxon with Cliff’s delta for two levels, Kruskal-Wallis with epsilon-squared for more. Multiplicity is corrected within each variable, not across all of them.
- Every signed effect carries a 95% percentile bootstrap interval. On ten patients most intervals span zero, and that is the finding rather than a shortcoming of the figure.
-
stats_clinical_association(),fig_clinical_heatmap(),fig_clinical_detail(),fig_clinical_forest(),fig_clinical_correlogram(),fig_clinical_landscape()andfig_clinical_trajectory()are exported. -
fig_clinical_correlogram()shows the clinical variables against each other. Correcting within each variable assumes they are separate questions; where two of them describe one gradient, a population associated with both is one finding counted twice, and nothing else in a run would show that. - This is association and not survival analysis. A 28-day flag is tested as the two-group comparison it is, and no time-to-event model is fitted, because the sheet carries no follow-up time.
Between-group differences, all populations on one figure
-
group_differences.png(fig_group_volcano()): Cliff’s delta against -log10 p, one point per population.group_comparison.pngis the figure for reading one population carefully; this is the one for deciding which to read. - It draws the smallest p-value the design can reach. A rank-sum test has
choose(n1 + n2, n1)rank arrangements under the null, so nothing below2 / choose(n1 + n2, n1)is attainable, 0.016 at 4 against 5. Where that floor sits above 0.05, an empty upper region of the figure is a property of the design and the figure now says so. -
population_trajectories.png: with--paired-columnand--condition-column, one line per patient across the conditions with the median over them, coloured by the first two-level clinical column. A box per timepoint averages away the direction each patient moves, which is usually the signal.
Report and figure fixes
- Every table in
report.htmlis its own toggle, closed, with a one-line description of what it holds. A section with a dozen tables of several hundred rows no longer separates its figures by yards of scrollable grid, and table bodies are built on first open rather than all at load. -
populations_by_batch.pnghas a key.facet_wrapdraws x-axis labels on the bottom panel of each column only, so on a three-row grid two thirds of the panels had coloured boxes with nothing naming them. - Categorical strata with exactly two levels are drawn blue and amber rather than green and red. Levels arrive alphabetically, so a survival flag came out “no” green and “yes” red, a value judgement the wrong way round, asserted by the palette rather than by anything in the data.
- The pooled marker UMAP in
marker_umaps_by_group/gets the same border and grey strip as the split ones, from the same code path. -
paired_comparison_stats.csv,covariate_adjusted_stats.csvandsubcluster_marker_shifts.csvare named by the between-group section instead of reaching the catch-all, whose purpose is to be empty. - Captions on figures whose canvas is sized from the data are wrapped to that canvas. A caption is drawn at one size with no wrapping of its own, so a line longer than the canvas was silently clipped.
- The licence in the documentation sidebar links to the full text. pkgdown prints the
Licensefield verbatim, andGPL (>= 3)is not one of the abbreviations it autolinks, so the site named a licence and gave no way to read it.
A sample that fails staining QC is no longer assumed to be a control
-
staining_verdict()reads an absentis_controlcolumn as “might be a control”, because onlyis_control = FALSEpositively asserts that a file is a biological sample. On a sheet without that column, any sample whose CD45 gate found no density minimum was promoted to the unstained reference for its whole panel. - Every threshold in that panel then came from the 99.5th percentile of that one sample. When it was really a stained sample, every cut landed at the top of a real distribution and every population collapsed. Measured on a real cohort: T cells at 0.034% of CD45+ against 4.4% once corrected, with no error raised and no warning printed. The frequency table looked like a finding.
- A control is now believed only when the sheet declares one. Absence of evidence is not evidence of an unstained tube, and a cohort with no controls is the ordinary case.
- When samples would have been promoted, the run now names them, says they are being treated as failed samples instead, and points at both remedies.
-
--discover-controlsrestores the previous behaviour. - Where a control reference IS used, the log now states that every threshold in the panel is a 99.5th percentile of it, and what a wrong reference looks like in
thresholds_used.csv. - If you have results from an earlier version on a cohort with no declared control, check
thresholds_used.csv. Asourceofcontrol_q995on most markers means those frequencies came from this failure and should be re-run.
Knowing what to write in the config
-
--list-channelsprints one row per acquisition parameter: index,$PnN,$PnS, and the symbol the run resolves it to, which is the name a population specification has to use. On a spectral panel those differ, because$PnSreadsCD45 : SparkUV-387 - Areaand the run usesCD45. Reads headers only. - It reads every file rather than the first, and reports whether they all carry the same panel, naming the markers that differ.
A channel that splits the cohort
-
--ignore-channelsdrops named channels before anything else, including before the panel fingerprint. - This matters more than it sounds. The fingerprint is the set of resolved marker names, so a channel present in some files and absent in others makes them separate panels, each with its own cofactor, embedding and thresholds. Autofluorescence channels from spectral unmixing do this routinely: they are not stains, and how many are written depends on the sample. Twelve comparable files can become seven panels of one or two files each.
- Names match whole and case-insensitively; a name containing
*is a glob, so[AF color*catches all of them andCD16does not catchCD161.
Testing what is testable
-
design_feasibility.csvis written whenever a group column resolves. One row per group withn_samples,n_donors,will_be_tested,valid_unpairedand a reason. - It catches two failures that are otherwise invisible. A group below the minimum is skipped silently, and a skipped comparison is indistinguishable from one that ran and found nothing. Worse, a donor contributing to more than one group makes an unpaired test treat repeated measures as independent, and that test runs and produces a well-formed p-value answering a question nobody asked.
-
--no-group-testskeeps the grouping and skips the between-group tests. Figures stay split and coloured by group and every per-sample quantity is still reported. For a cohort worth describing but too small, or too confounded, to test. Diagnostics are unaffected. -
--min-group-nexposes the minimum group size, which was fixed at 3.
The t-test and ANOVA, with their assumptions attached
-
parametric_tests.csvholds Welch’s t-test and Student’s t-test for two groups, Welch’s ANOVA and the classical one-way ANOVA for more, with Cohen’s d and eta squared. -
posthoc_tests.csvholds Games-Howell, Tukey HSD and Dunn for every pair. - All of it is computed on arcsine square root transformed percentages, the standard variance-stabilising transform for proportion data, and both scales are reported because a difference in arcsine units is uninterpretable.
- Every row carries
shapiro_p,brown_forsythe_pand arecommendedcolumn. Normality is tested on within-group residuals rather than pooled values, which are bimodal whenever the groups genuinely differ. - These sit beside the rank tests and never replace them. No existing table changes.
--no-parametricturns them off.
Figures are written on white, not on transparency
- Every figure now carries an opaque white background.
ggsave()takes its background from the theme, and the theme leftplot.backgroundunset, so PNGs were written RGBA with a transparent ground. Measured on one run: 11 of 20 figures. - That is fine inside a white report and unreadable anywhere else. Expand a figure in a dark viewer, or drop it on a dark slide, and black axis text sits on whatever is behind it.
- Set in three places:
safe_ggsave(), which every one of the 43 figure writes passes through;theme_cyto(), so the plot object itself carries a white ground however it is rendered; and the rawpng()device the multigraph overlay draws on, which does not go throughggsave()at all. A caller passing its ownbgstill wins.
Legends sit top-left on the large multi-panel figures
-
functional_markers.png,population_frequencies.png,population_ratios.png,group_comparison.png,marker_state.pngandumap_multigraph_overlay.pngput their key at the top left instead of the right. - patchwork collects the guide once for a whole composition and centres it vertically on the right. On a tall grid that is thousands of pixels from the top:
functional_markers.pngruns to 26 panels and six thousand pixels, so finding out what a colour means meant scrolling away from the panels and back. -
This did not work on the first attempt, and the figures shipped with the key still on the right. Setting
legend.positionon a panel has no effect onceguides = "collect"lifts the guide out of that panel: the guide is drawn against the composition, which falls back to its own default. The position has to be set on the patchwork theme, and now is. -
umap_multigraph_overlay.pnghad a second, worse form of the same problem. Its cluster key was built withncol = 1, so a twelve-cluster cohort stacked twelve rows into the height reserved for the UMAP row and collapsed those panels to a sliver. The key now flows across in rows of six, and the row’s height grows with it. -
threshold_drift.pngalso moves to a top-left key. An earlier note here said it kept a right-hand one “for the same reason” as the UMAP figures; that was wrong, because it is a facet grid whose height grows with the marker count and not a single panel. - The standalone UMAP figures are unchanged. They are single panels where the key already sits beside what it describes.
- Captions on those grids are wrapped to the figure’s own width. The line naming the groups reached 205 characters on a three-cohort study and was silently clipped mid-name at the canvas edge.
Group colours are separated hues, and no two groups share one
- Study, cohort and timepoint groups take green, red, purple, amber, blue, pink in that order, with green anchoring the reference. Each is roughly a sixth of the colour wheel from the last.
-
A fourth group used to be drawn in the reference’s green. The palette held three colours and was recycled with
rep_len(), so with four or more groups two of them came out identical, with nothing said about it. Past the palette the run now generates additional hues instead of repeating. The same bug and the same fix applied to the multigraph overlay, where overlaid curves sharing a colour are unreadable rather than merely ambiguous. - Small non-study category sets no longer get red, orange and yellow. Those are the first three entries of the population palette, which is ordered by hue, so taking the first n handed two or three categories three neighbouring colours: the hardest triple to separate at bar-chart size and nearly identical under red-green colour vision deficiency. Up to six categories the picks are now spread across the palette instead.
-
fig_group_comparison()passes the reference it resolved rather than relying on the global default, so a figure’s colours no longer depend on what an earlier call happened to leave set.
Long panel keys no longer overrun the panel beside them
- A panel’s key is its y-axis title, and a y-axis title is rotated, so the space it has is the panel’s height rather than its width. Population names fit. The keys
functional_markers.pnguses areblock: population, marker, reaching 60 characters, and they ran out of the panel onto the neighbour’s axis numbers and tag letter. - Long titles now wrap onto two or three lines, with the left margin growing to match, so the tag letters stay legible. Short titles are untouched, so
population_frequencies.pngandpopulation_ratios.pngare unchanged. - Lines are balanced rather than wrapped at a fixed width: wrapping a 47-character title at 30 leaves one full line and a stub, which reads badly rotated.
A full-size UMAP for every marker
-
marker_umaps_by_group/holds full-size UMAPs per marker, written by every run.umap_<marker>.pngpools every sample, andumap_<marker>_by_<category>.pngis written beside it for every category present, not only the one named by--group-column. - Naming a group column decides what is tested. It says nothing about what is worth seeing, and a cohort usually carries several categories: timepoint, infection focus, phenotype, sex. Before this, a run grouped by timepoint drew a single merged panel for infection focus even though the column was in the sheet and already panelled on
umap_overview.png, and an ungrouped run drew no split at all. - Pooled and split are both kept, because they answer different questions. Splitting asks whether a marker sits differently between categories; pooling asks where the marker is at all, which no split figure can show, because each facet holds only a subset of the cells.
- What is never faceted: numeric columns, which would give one panel per distinct value; identifiers such as
sample_idandpatient_id, which are not categories and would multiply the folder by the donor count; and single-level columns, since one facet is not a comparison. Past four categories the extras are named in the log rather than dropped silently, and--group-columnmoves one to the front. - Neither existing figure covers this.
umap_markers.pngcolours by intensity but pools the groups and shrinks each marker into a grid cell.umap_overview_by_group.pngsplits by group but colours by population or covariate, never by intensity. “Is this marker brighter, or somewhere else, in one group” needs both at once. - Within one file the colour scale is shared across the facets and clipped to the 1st-99th percentile of the cells drawn, so a colour difference between two panels is a real intensity difference. It is not shared between files, for the reason
fig_marker_grid()gives. -
--no-marker-group-umapsturns it off. Two PNGs per marker on a grouped run is a folder you may not want on a panel of forty.
Batch correction per cell type, not per file
-
--batch-method clusterfits one alignment map per marker per cell type instead of one per marker over the whole file.cytonormis accepted as a synonym, because that is the name the method is known by.quantile, the whole-file map, stays the default and is unchanged. - Why it matters: a detector shift moves a bright population and a dim one by different amounts, so a single map fitted to the pooled distribution over-corrects one and under-corrects the other. A per-file map also moves with the biology, so it can remove the difference it was meant to preserve.
- It is fitted in two stages, and the second stage is not optional. Clustering the raw matrix lets a large batch shift become the dominant variance, so the clusters turn out to be the batches, every cluster holds one batch and the correction silently does almost nothing. Measured on a synthetic three-batch shift: per-cluster alignment fitted that way left a mean between-batch gap of 1.296 against whole-file alignment’s 0.003. Fitting the clustering on a whole-file-aligned copy and the per-cluster maps on the original values brings it to 0.014, while leaving the true between-cell-type separation at 3.996 against a true 3.994, where whole-file alignment inflates it to 4.078.
-
--batch-cluster-ksets the number of cell types, default 10. - No new dependency. The clustering reuses
run_unsupervised_clusters(), which is already seeded, already stream-safe, and already falls back to a built-in SOM when FlowSOM is absent. -
The Cramer’s V refusal is evaluated before the method is consulted. Both methods are refused on identical evidence, and
--force-batch-correctionremains the only way past. A better alignment algorithm does not make a confounded design correctable. -
batch_correction.csvis new. The decision was previously held in memory and never written, so no output said whether a correction had run, which method fitted it, or why it was refused. Written whenever--correct-batchis set, including on a refusal.
--list-channels now names the fix, not just the symptom
- When files carry different panels, it reports which channels vary and, where the cause is mechanical, prints the exact
--ignore-channelspattern that merges the cohort. - Two causes are detected: autofluorescence and unmixing artefacts, which are not stains and whose count varies per acquisition, and a real marker present in a minority of files. Anything else is reported as a genuine panel difference and left alone, because that is a decision about the experiment.
- This was previously something you had to know to look for. A fragmented cohort raises no error: every file loads and every table is written, while each panel quietly gets its own cofactor, thresholds and embedding.
Speed
-
--read-threads Nreads N files at once. Reading is the slowest stage on a large cohort and files are independent, so this scales close to linearly. Peak memory scales with it too, since each worker holds one expression matrix. - Everything after reading stays sequential, because cofactors, thresholds and the embedding are derived across samples from a seeded generator that has to be consumed in one order for a run to reproduce. Results do not depend on the thread count.
Fixes
- A UMAP marker panel whose marker is constant over the cells drawn no longer fails the run. The 1st and 99th percentile came back equal, the colour scale was zero width, ggplot collapsed its two identical breaks into one while both labels remained, and the run died with “
breaksandlabelshave different lengths” after most of its output had been written. -
--paired-columnwithout--condition-columnwas documented as sufficient in places. It is not, and never was: the pipeline logspaired test skippedand writes no paired output.
Explore mode
-
--exploreruns unsupervised discovery over every eligible channel – scatter and viability included. ignoring the population specification and the parent gate. It is what the declared path structurally cannot do: find a population nobody declared, and look outside the parent gate. - Everything it writes lands in
<outdir>/explore/, including its own self-containedexplore_report.html. No existing output changes. That is verified rather than asserted: the demonstration cohort is run with and without the flag and every top-level deliverable compared byte-for-byte. -
--explore-onlyneeds no specification at all.--diris the whole input. For a panel you have not written a config for yet. - A cluster-level quality gate, the one the unsupervised family uses: cluster coarsely, then judge whole clusters by their marker profile, so no per-event cutoff is invented. Calls debris, dead and saturated, and records the basis for each in
explore_qc_clusters.csv. - Cells are equalised per sample after the gate, so a faceted panel’s density is comparable between groups rather than reflecting group size.
-
explore_suggested_spec.yamlturns the clusters into a draft specification to curate and run back through the supervised path.
--maybe-learn
- The only thing that lets the two analyses see each other, and it is off by default. An explore run that used the declared thresholds is no longer a blind unsupervised run, and a check on a specification has to be independent of the specification to be worth anything.
- With it, the declared run lends explore its per-sample thresholds, so a cluster gets a phenotype measured against each sample’s own cut –
CD19+ HLA-DR+ CD3-, with a percentage behind every call. instead of one read off a colour scale. Its staining QC verdicts and its batch/group confounding verdict travel too, the latter into the cluster statistics, so a q-value cannot be read without it. - In the other direction, explore writes
spec_gaps.csv: declared populations spanning several clusters, and clusters no population covers. The first is the one worth the flag. a population whose total is flat while a subset inside it moves is invisible to any amount of testing on the declared table.
The statistics catalogue
-
statistical_methods.csv, on every run, names every commonly reported method in the immunophenotyping and cytometry literature. Student’s and Welch’s t, one-way and two-way ANOVA with Tukey, repeated measures, Mann-Whitney, Kruskal-Wallis with Dunn, chi-squared, Bonferroni, and diffcyt’s edgeR, limma-voom and GLMM. and says whether this run computed it and why. -
normality_tests.csvis the evidence for the choice: Shapiro-Wilk per population per group, Brown-Forsythe for equal variances across groups. Itsinterpretationcolumn states the thing a p-value alone hides. that Shapiro-Wilk on 4 to 10 donors has almost no power, so a non-significant result is not evidence of normality, which is itself the argument for the rank tests. - What deliberately does not happen: running every test and reporting all the p-values. At these group sizes that is p-hacking with extra steps.
Input format and the run report
Two changes to how a run is specified and how its result is read, and a version number that says the interface is now settled. Every option, output file, column name and value that 0.4.0 produced is produced identically; nothing is renamed and nothing is removed.
One CSV and one YAML
- A run needed up to three tables: a sample map keyed by filename, a patient table keyed by patient identifier, and an absolute-count export keyed by whichever of the two the counting instrument happened to write. Each is a different shape with a different key, and the joins between them are where cohorts get mislabelled: a patient present in one and absent from another produces no error, only an empty covariate panel or a silently dropped count.
-
--samplestakes one sheet with a row per acquired file, carrying every fact about that file: what it is, whose it is, which group and batch it belongs to, and any externally measured counts incount.<Population>columns. One key, one file, no joins.--configcontinues to carry the analysis. Between them the two files specify a run completely, and the config’s newsamples:section can name which column is the group and which the batch so that no further flag is needed. - The split is deliberate. Anything varying per sample belongs in the CSV; anything that is one decision for the whole study belongs in the YAML. Putting an analysis choice in the CSV would repeat it on every row and invite the rows to disagree about it.
- The sheet is a different way to supply the same facts, not a different analysis. The reader splits it into the same three structures the pipeline always consumed and applies the same coercions by calling the same code:
normalise_patient_columns()was factored out ofload_patient_table()so the two routes cannot diverge in how they parse a date or translate a value. The claim is measured rather than asserted: one study run both ways produced 31 of 31 tables and 22 of 22 figures byte-identical. -
--sample-map,--patient-tableand--absolute-countsstill work exactly as before and are not deprecated. They cannot be combined with--samples, because two sources of truth for one fact is what the format removes. - One hazard the three-file shape did not have. A subject attribute is a property of a patient, but the sheet has a row per FILE, so a patient with several acquisitions repeats it. Two rows can therefore disagree about that patient’s sex. A reader taking the first value would silently pick one, so
read_samplesheet()reports every disagreement assubject / column: value vs valueand stops. A blank is not a disagreement and is filled from the rows that carry a value. -
--write-sampleswrites a sheet covering every file in the input directory, with the filename-derived identifiers filled in. Identifiers are never inferred from plate order at run time: a wrong guess mislabels a patient and nothing downstream reveals it, so a file the sheet does not cover is a fatal error. - Annotated templates ship as
inst/examples/samples_template.csvandinst/examples/analysis_template.yaml.
Validation before the run rather than during it
-
--checkreads FCS headers and the two input files, reports what the run would do, and exits without analysing anything. It names every marker resolved from the files, every specification entry matching none of them, whether the sheet covers every file, the levels of the group column and their sizes, the number of batches, and the study variables available. - Everything it catches was already knowable before the first event was read, and was already being reported: partway through a run that costs many minutes on a real cohort. A marker name not matching
$PnSexactly is the most common cause of an empty frequency table, and it was surfacing after the reading stage rather than before it. - It reads keyword blocks rather than event matrices, so it costs seconds regardless of cohort size.
The report is one self-contained file
-
report.htmlembedded nothing. It referenced the figures beside it, so moving it produced a page of broken images, and a result that cannot survive being moved is not a record. - Every figure is now embedded at full resolution and every table in full. The file references nothing, loads no font or script from a network, and works from a
file://path. It can be attached to an email, put in a supplement or archived on its own. - Base64 encoding is written out in
R/base64.Rrather than taken from a package: the CRAN options would each add a dependency to an import list that is deliberately short and a container image that is pinned, for fifteen lines with no edge cases beyond padding. It is vectorised over the raw vector, because a per-byte loop over a 6 MB PNG is minutes in R. - Sections are collapsible, with a sidebar indexing every section, figure and table. Figures share one display box whatever their native aspect ratio, zoom on click with keyboard control, and download at full resolution from the same bytes that are already embedded, so nothing is stored twice. Tables are searchable, sortable by column, shown 10, 50, 100 or all rows, and export to CSV exactly as filtered and sorted.
- Tables are carried as JSON and rendered in the browser rather than written as markup. Over pre-rendered rows the export would re-parse the DOM and the search would hide rather than filter; carrying the data means all three operations read one array, so what is exported is what was filtered.
- Section headings state what the section reports rather than posing a question, and each carries a description of what its figures and tables show and how to read them.
- No output can be silently omitted. The named sections put files in reading order; anything a section does not name is collected under “Further outputs” rather than dropped. A table above 8 MB is named with its row count instead of being embedded, which in practice means only the per-cell exports, whose row count is the number of cells; the limit is
options(cyRAVEN.report_table_max_mb = ). - The file is large, tens of megabytes on a full run, and its size is logged. That is the cost of self-containment and is stated rather than discovered.
A failed run explains itself
- A run that failed left the outputs written before it stopped, a manifest marked
failed, and an R error on stderr that whoever opens the results directory afterwards may never have seen. Reconstructing what happened meant reading a stack trace out of a container log. - A failed run now writes
report.htmltoo, with the diagnosis first: what stopped it, which stage it reached, the log leading up to that point, what the message means in the vocabulary of the analysis rather than of R, and the specific next action. Everything produced before the failure is embedded below it, because the partial output is usually where the evidence is: agating_qc.pngwritten before the failure often shows the cause directly. - Twelve failure modes carry an interpretation, including a sheet that does not cover every file, subject rows that disagree, a path that does not exist inside the container, memory exhaustion, an absent optional package, and a marker name that resolves to nothing. “Subscript out of bounds” is not a useful thing to tell a cytometrist; “the specification names a marker this panel does not contain” is the same fact in terms they can act on.
- The report is written by the entry point rather than by an exit handler, because an exit handler cannot see the condition that ended the call and “something failed” is not a diagnosis. The original error is re-raised immediately afterwards, so the exit status and the message a caller sees are unchanged, and a failure inside the reporting cannot replace the real error with one about reporting.
- The run log is now kept in a bounded in-memory buffer as well as being written to stderr, which is what lets the report show what the run was doing when it stopped.
Documentation
- New Inputs article covering both routes, every reserved column, the count columns, the config sections, the templates, and the errors the format can raise with what each means.
- The Get started article is organised by use case, each with the full command sequence: a first run on a new cohort, the full analysis, a cohort that must fit in memory, comparing against an accepted baseline, and driving it from R.
- The diagnostics article opens with
--check, since every other check in it applies to a run that has already happened. - README, the options reference, the output article and the Claude skill all document the new flags, the report’s behaviour and the failure report. # cyRAVEN 0.4.0
Seven additions, completing the feature backlog in dev/TODO.md. Every file the previous version wrote is still written with the same name, the same columns and the same values, verified by running both versions on the same cohort and comparing byte for byte: 34 of 35 outputs identical, none lost, and the one that differs is miflowcyt.md, differing only in its own generation timestamp.
Height-only acquisitions
- The reading stage kept area channels only, because a marker recorded as both area and height would otherwise be counted twice in the embedding. Applied to a file that records height and no area, which older instruments and several clinical archives do, that rule resolved zero markers. The run then failed several stages later with every population reported UNAVAILABLE, which points at the specification rather than at the panel.
-
read_fcs_resolved()now falls back to pulse height when a file contains no area channel, and says so. The rule exists to avoid double-weighting a marker that appears twice; where nothing appears twice there is nothing to avoid. The substitution is announced rather than silent because height understates a bright wide event, so thresholds derived from such a file are not interchangeable with those from an area acquisition of the same panel. - The scatter gate is defined on area and stops without it, so a height channel now stands in there too. That alone would have made
derive_singlet_band()divide a channel by itself: a ratio identically one, a MAD of zero, a band of zero width, and every event discarded. The singlet gate now requires FSC-H and FSC-A to be distinct columns and skips with a log line when they are not.
Demonstration cohort
- The worked example now uses the graft-versus-host disease dataset distributed with flowCore (Brinkman et al. 2007, Biol Blood Marrow Transplant 13:691), replacing two CytoTrol control acquisitions partitioned into pseudo-samples with randomised group labels.
- The previous cohort could demonstrate calibration and nothing else: its groups were randomised, so the correct answer was always that nothing differs, and its panel carried neither CD45 nor a viability dye. The replacement has five patients, seven visits as genuine acquisition batches, a real clinical contrast between GvHD grades, CD45 in the panel, a Time channel, and 2,205 to 66,105 events per file.
- It requires no download. The data ship inside a package cyRAVEN already depends on, so the example is reproducible offline and cannot break when a repository moves or a certificate expires.
- Both groups are transplant recipients, so the contrast is disease severity rather than disease against health, and the grades are unbalanced across patients.
inst/scripts/demo_data.Rstates both, and the confounding diagnostic is the check to read before interpreting the contrast.
Documentation
- The demonstration configuration now declares two
functional_blocks:and one entry underratios:in addition to its five populations, so a default run writes all 22 figures the package can produce from a cohort with no volumetric counting, rather than the 20 that thepopulations:block alone reaches.absolute_counts.pngandabsolute_counts_qc.pngare the two it cannot, and the worked example says why and documents the measurement they need instead of filling the gap with invented counts. -
functional_blocks:andratios:had no syntax reference. Both are now documented in the Gating article and in the skill, including the requirement that decides whether a block is meaningful: a marker must not be read inside a gate its own threshold helped draw, because a percent positive pinned at 100 in every sample has zero variance and measures the definition rather than the biology. - The
--absolute-countsinput format is documented in the Output article: the sheet layout, how the unit label is read, and the two-key sample matching with its ambiguity rule. - The worked example carries every figure the run produces, each with what it measures and what this cohort shows, and each with alt text.
- New
options.Rmd, a reference for all 83 command-line options with their defaults and consequences. An audit found 26 of them documented nowhere in the README, the articles or the skill. It also states the convention that decides which options are on by default: analyses that only add output are, and the five that change numbers a previous run reported are not. - Maintainer working files moved to
dev/, which is excluded from the build, the container context and version control.dev/README.mdrecords which documents are authoritative where they disagree.
Acquisition-time quality control
- Every number a run reports for a sample is derived from that sample’s pooled events, which is correct only if the instrument was doing the same thing throughout the acquisition. The Time channel was read and then discarded, so nothing checked. A partial clog, a bubble or a drift in laser power makes a file two instruments over its run, and one threshold then suits neither half.
- New
acquisition_qc.csv,acquisition_qc_bins.csvandacquisition_qc.png. The Time channel is binned into equal-width intervals and two quantities are tracked: the event rate, which a clog lowers and a bubble spikes, and each channel’s median, which catches a shift that leaves the rate untouched. Both are judged by robust z against the file’s own bins. - Bins are equal in TIME rather than equal in count. A bin defined to hold a fixed number of events cannot show that the rate changed.
- Nothing is removed.
acquisition_qc_impact.csvstates how far each population would move if the flagged intervals were excluded, which is what turns a flag into a decision: a file with a visible clog that moves no population by more than its own gate uncertainty does not need re-acquiring.--drop-unstable-eventsperforms the removal and is recorded in the manifest. - The detector is written directly rather than wrapping PeacoQC, which would add a Bioconductor dependency for a binned robust location test. Emmaneel et al. 2022, Cytometry A 101:325 is the reference for the approach.
- The figure is drawn at the end of the run rather than where it is computed. Drawing it early changed the rendering of nine existing figures; the data behind them was provably identical either way, but a new diagnostic that silently re-renders every published figure is not an addition.
Fluorescence-minus-one controls
- Control handling was one unstained tube per panel, used as the reference for every marker. That is the right control for where autofluorescence ends and the wrong one for where a marker’s background ends in a panel, because it cannot show spillover. An FMO is the same panel with one reagent left out, so its distribution in that channel is the negative population under the spreading the real samples experience.
- Two optional sample-map columns,
fmo_forandcontrol_group. The first names the markers a file controls for; the second confines a control to the batch it was acquired in, since a reagent lot changes between batches. - New threshold source
fmo_q995, taking the same place in the resolution chain as the unstained control and named apart from it because the two are different experiments supporting different claims. - New
fmo_agreement.csv, which is the reason to have the feature rather than just the control. It reports the distance between the derived cut and its FMO-anchored equivalent in units of that threshold’s own uncertainty. Within about one they agree to the precision either can claim, and the derived cut is corroborated by an independent experiment. Beyond about three one of them is wrong: far above the FMO the cut is discarding signal, far below it is calling spillover positive.
Per-sample threshold overrides
-
--configcould pin a threshold for the whole run, which corrects one tube by applying one number to every sample and so reintroduces the fixed-coordinate bias the package exists to remove. Asample_overrides:block now corrects one sample and one marker. - New source
manual, distinct fromconfig: the first says a named person moved this one cut for a stated reason, the second says the assay declares this cut everywhere. -
thresholds_used.csvgainsoverride_reasonandoverride_by, and the run manifest lists every override. Both are absent on a run that declares none, so that run writes the table it always wrote. - An override matching no sample or marker in the cohort is reported rather than ignored. Silence there means the analyst believes a cut was corrected, the run says nothing, and the uncorrected number is published.
-
specification_conformance.csvreports a hand-set marker asmanually setrather than as agreement. A marker made to match the baseline by hand did not conform.
Spillover spreading
- New
spreading_pairs.csvandspreading_receivers.csv. Spreading is the one optical effect that removes a density minimum without moving anything else, so it is the most common reason a marker that should resolve does not, and the run could previously only reportquantile_fallbackwith no cause. - For each ordered channel pair, the spread of the receiver’s negative population is compared between cells negative and positive for the source. Restricting to the receiver’s own negatives is what makes it spreading rather than biology: co-expression moves the positive cells, not the width of the negatives.
- The actionable join is to the fallback rate. A marker that both fails to resolve in most samples and receives substantial spreading is reported as a panel design problem, which no gating strategy fixes. On the demonstration cohort CD38 and CCR7 are both flagged that way.
- This is a ranking from the samples in hand, not a spillover spreading matrix: the published SSM is computed from single-stain controls, which this package is not given.
Run report
- New
report.html, written by default. A run writes several dozen files and the reading order is the product: each check can invalidate the ones after it, and a directory listing presents a failed staining QC and a headline p-value as equals. The report puts them in order, states what each section answers, and links every claim to the file it came from. A run that excluded samples says so above every result. - Written as HTML directly rather than through rmarkdown, which would need pandoc: a report the documented execution path cannot produce would be worse than none. Figures are referenced rather than embedded, so the report is read in the results directory beside them.
Opt-in changes to reported numbers
Both of these change existing numbers, so both default to the previous behaviour and are recorded in the manifest.
-
--subsample rareweights the embedding draw by inverse local density. A population at 0.3% contributes about 60 cells out of 20,000 under a uniform draw, which is enough to be a smudge and not enough to be a cluster, so the unsupervised cross-check cannot recover it. That makes the sampling the binding constraint on this package’s own falsification claim for rare populations.cells_umap.csvgainssampling_weightso anything computed from the embedded cells can be weighted back to the true composition. -
--calibration-beadswith--calibration-valuesconverts channel units to MESF or ERF before any threshold is derived. Arbitrary units are comparable within an instrument and never between two, which is the ceiling on what--baselinecan establish: it can say a cohort moved and not whether the assay or the cytometer did. A channel whose fit does not hold is left in instrument units and said to be, because a calibration nobody checked converts an honest arbitrary number into a dishonest absolute one.
cyRAVEN 0.3.0
Counting uncertainty and detection limits
- Every frequency now carries a second uncertainty beside the gate one: what it gets from the number of events behind it. The budget in 0.2.0 perturbed thresholds and nothing else, so a population of twelve events and one of twelve thousand, both sitting behind the same clean valley, were reported with the same uncertainty. Displacing a well-separated cut moves neither of them much. One of those numbers is worth acting on.
- New
counting_uncertainty(). The interval is Wilson rather than the textbooksqrt(p(1-p)/n), which goes to zero as p does and so claims perfect knowledge of a population nobody observed. The two agree to three figures once a population has a few hundred events, so nothing is given up in the common case to fix the boundary. -
population_frequencies.csvgainsn_parent_events,u_counting_pct_points,u_total_pct_points,pct_lo_total,pct_hi_total,lod_pct,loq_pctanddetection.u_pct_pointsstill means gate placement alone and still holds the value it held before, andpct_lo/pct_hiare still built from it. Only the new columns are new. -
group_comparison_stats.csvgainstotal_u_pct_pointsanddifference_over_total_u.difference_over_gate_uis untouched. The two separate for a rare population, where the cut can be well placed and the frequency still rest on too few events: the first ratio passes and the second does not. - The limits follow the clinical event-count convention,
--lod-events(20) and--loq-events(50), against the parent-gate events THIS run saw. That means--max-events-per-fileraises every limit in proportion, which is the honest reading: a population under the limit of a subsample may be well resolved in the whole file, and the answer is to raise the cap rather than to believe the limit. - New
detection_limits.pngcounts, per population, how many samples clear each limit. A population mostly below the limit of quantification cannot be fixed by changing the gating, because the numerator is small for want of cells. One split across the limit is the dangerous case, since its group difference can be driven by which samples happened to clear it. - Counting is deliberately absent from
uncertainty_budget.csv, which answers which THRESHOLD a population’s uncertainty comes from. Counting is not a threshold, and adding a term to that table would have changed a published figure to say something the file does not mean. - The two components are not strictly independent, since displacing a cut also changes the count. Quadrature treats them as though they were, which is the same approximation the GUM makes for the marker terms; it is stated in the documentation rather than left to be discovered.
MIFlowCyt report
- New
miflowcyt.md, written by default beside the run manifest. MIFlowCyt is the ISAC reporting standard and is checked at submission by Cytometry A, Nature and PLOS. Its hardest section to write by hand is the one this package already knows in full – the gating specification, the transform and its derived parameters, and how every threshold was obtained – and its instrument section was sitting unused in memory, becauseread_fcs_resolved()keeps the whole keyword block and the pipeline read four keys out of it. - Sections 3 and 4, instrument and data analysis, are completed from the keyword block and from the run: cytometer, serial, acquisition software, date, operator and event counts per file, then per panel a detector table carrying
$PnN,$PnS, voltage, range and whether amplification was linear or logarithmic. - Sections 1 and 2, experiment intent and specimen biology, are marked TO BE COMPLETED rather than omitted. An absent section reads as one that did not apply; a marked one reads as work outstanding, which is what it is. The same rule governs individual keywords: one the file does not carry is written as “not recorded in the FCS file”, because its absence is a fact about the acquisition and is often the reason an analysis could not be done.
- Compensation is reported as applied, absent, or mixed. Mixed is a real and dangerous state – some files carrying a
$SPILLOVERmatrix and others not – and both of the other two answers would be false for it. -
--no-miflowcytskips it.
Per-marker batch drift
-
--batch-columnnow also writesmarker_batch_drift.csvandthreshold_batch_drift.csv. iLISI says whether the batches mix in the shared embedding, which is the right question and names nothing anyone can act on. A flagged marker names a reagent lot, a detector voltage, or a laser that was serviced between runs. -
threshold_batch_drift.csvcost nothing to add:stats_threshold_drift()was always generic in its grouping vector, and the pipeline had only ever handed it the study group. It is the same test, grouped by batch. -
marker_batch_drift.csvis the part that is new. A threshold is one number per sample, so the test above only sees drift that moves the CUT. A marker can change its spread, grow a tail, or lose the separation between its modes while the density minimum between them stays exactly where it was, and no threshold statistic can detect that. Newmarker_batch_drift()compares the distributions themselves. - The statistic is the one-dimensional Earth Mover’s distance, new
emd_1d(): the L1 distance between two empirical quantile functions. No binning choice, no dependency, and it is in the units of the analysis scale, so it is reported as a distance rather than as a score. Since those units differ per marker it is also divided by the marker’s own pooled MAD, which is what makes two markers comparable;emd_over_madat or above 0.5 is flagged. - Its subsampling borrows and returns the RNG stream, for the reason in 0.2.0.
Fixes
- The RNG regression test was testing
withr, not cyRAVEN. Its fixture usedwithr::local_seed(), and withr 3.0.2 restores the generator’s position but not its state vector, so the fixture leaked a changed stream into a test whose entire purpose is to prove nothing leaks. The fixture now saves and restores by hand, the same way the package does. The package’s own guard was correct throughout and is unchanged; verified directly by callingthreshold_uncertainty()on a pre-computed vector and comparing.Random.seedand the next five draws either side.
cyRAVEN 0.2.0
Three additions. All three are additive in the strict sense: every file the previous version wrote is still written, with the same name, the same columns in the same order, and the same values. Verified by running both versions on the same cohort and comparing byte for byte, which is also how the RNG hazard below was ruled out.
Gate placement uncertainty
- Every reported frequency now carries a standard uncertainty saying how far it moves when the thresholds behind it move. Two components, combined in quadrature per the GUM convention: resampling the parent-gate events, and sweeping the settings
density_valley()itself takes. The second is the honest form of the claim that automation removed analyst discretion. It did not, it fixed it, and this is how much was fixed. - The parent gate is a term in every population’s budget, because displacing the CD45 cut moves cells into and out of the denominator all of them are expressed against. On the development cohort it is the largest median contributor, which is what the operator-variability literature reports for manually gated data as well.
-
group_comparison_stats.csvgainsdifference_over_gate_u. Below 1, the groups differ by less than the distance the cut itself moves. A p-value does not disclose that. - New:
threshold_uncertainty.csv,uncertainty_budget.csv,frequency_uncertainty.png,uncertainty_budget.png. On by default;--no-uncertaintyrestores the previous output exactly. - The reported threshold is never touched. Perturbation runs on copies, so the numbers are identical with the analysis on or off.
- The RNG is borrowed and returned.
run_cyraven()seeds once and STEP 6’s UMAP cell selection draws from that one stream, so a new step consuming draws would have changed which cells are embedded and silently redrawn every UMAP in the run. Every entry point saves and restores.Random.seed; the embedding is bit-identical either way, and a test asserts it. -
bootstrap_valley_raterecords how often a resample finds the cut at all. A threshold that comes and goes is not imprecise, it is unresolved, and it looks identical to a solid one inthresholds_used.csv.
Conformance against a baseline
- New
--write-baselineand--baseline. The existing peer check inthresholds_used.csvcompares each sample against the others in its own run, which catches one bad tube. It cannot catch a cohort that moved as a whole, because the leave-one-out peer median moves with it. A baseline written from an accepted run is the fixed reference that can. -
threshold_scale_qc.csvkeeps its within-run meaning. Repointing it at a baseline would have kept the filename and the column names while changing what the numbers mean, so the baseline comparison goes tospecification_conformance.csvinstead. - A changed transform is reported once as “not comparable” rather than as every marker having drifted, and a redefined population is reported as a different measurement rather than a moving one.
-
--fail-on-driftturns the verdict into an exit code for scheduled runs. It raises after the run is marked completed, because every output was still produced.
Labels from another tool
- New
--external-labels, which takes a CSV of cell labels this package did not produce – a cyCONDOR clustering, for instance – and learns a gating strategy for each one.learn_convex_gate()already accepted an arbitrary label; what was missing was a way to supply one, a way to find out whether the result transfers, and a file the instrument can read. - The join is on sample and event index, never position. Tools subsample independently, so row i of one table is not row i of another and a positional join mislabels every cell without erroring.
- New
gate_transferability()refits the whole strategy with one donor withheld and scores it on that donor, once per donor. The existing held-out metric reserves CELLS, which come from the same donors acquired in the same tubes, so it can look excellent on a gate that fails on the next patient. Report the minimum across donors, not the mean. - New
--export-gateswrites ISAC Gating-ML 2.0 plus a plain table of vertices, both in the linear units the FCS file stores. A gate that stays in R cannot be sorted on. - Polygon edges are subdivided before inversion. An edge that is straight on the analysis scale is a curve in linear units, so inverse-transforming only the corners would emit a different region from the one that was fitted and validated.
- No new dependencies. The XML is written directly rather than through a
GatingSet, which keeps the export available in the container this package tests in.
cyRAVEN 0.1.0
First packaged release. Previously a set of command-line R scripts; the analysis is unchanged, the delivery is not.
Standalone coverage
- New
--transform logiclealongside the existing arcsinh path. Logicle is the flow-cytometry convention: linear near zero, so compensated negative values stay on scale rather than piling against an axis limit. Parameters are pooled per PANEL, not fitted per file, so every sample stays on one ruler and cross-sample medians remain comparable. - The transform choice is not cosmetic. On a seven-colour T-cell panel the top-level populations agree to 0.1% between methods while the CCR7/CD45RA memory subsets move by tens of percent, because under arcsinh the CCR7 threshold resolves to a density valley in some samples and a quantile fallback in others. Logicle finds a valley in every sample.
thresholds_used.csvrecords which happened. - New
--correct-batch, which aligns each marker across batches by monotone quantile mapping. Being monotone means it cannot reorder cells within a batch, so it corrects location and spread without inventing structure. - Batch correction REFUSES when Cramer’s V between batch and group exceeds
--batch-max-cramers-v. Where the two overlap that far, removing the batch and removing the finding are the same operation.--force-batch-correctionoverrides and is recorded in the run manifest.
Learned gating strategies
- New:
--explain-clustersderives a gating strategy for any cluster the population specification does not describe.cluster_gate_agreement()could already report that a cluster matched nothing; it could not say what the cluster was, which left the finding as a cluster number and a cell count – nothing anyone could act on at the instrument.explain_cluster()learns a short sequence of two-marker gates that selects those cells and reports how well each performs. - Gates are convex polygons – intersections of half-planes – rather than the axis-aligned rectangles a conjunction of one-dimensional thresholds can express. Diagonal boundaries (CD4/CD8, CD14/CD16, FSC-A/FSC-H) cannot be captured by a rectangle without either admitting contaminants or rejecting real cells;
derive_singlet_band()was already a hand-written special case of exactly this, and a polygon is its general form. - No new dependencies. The objective’s gradient is derived in closed form and the fit is 24 parameters, so
stats::optim’s L-BFGS-B solves it directly with no automatic-differentiation package;stats::prcompinitialises andgrDevices::chulltightens, both already imported. - Every reported metric is computed on cells held out of the fit, with the in-sample value printed beside it. Eight free half-planes fitted to a few thousand cells can memorise them, so a resubstitution F1 is optimistic by construction and the size of that gap is information the reader needs.
- Strictly descriptive. A proposed gate says where cells sit in marker space, not that they are a real population; it carries no p-value, and the scored populations, frequencies and tested results are identical whether this runs or not.
Packaging
- Reorganised into a standard R package: documented, namespaced, tested, and installable with
R CMD INSTALL. The command-line entry point survives assystem.file("scripts", "cyraven.R", package = "cyRAVEN")and calls the samerun_cyraven()that is now exported. - Mutable globals (
COLORS,REFERENCE_GROUP) replaced by a private package environment with accessors:fcs_colors()/set_fcs_colors()andreference_group()/set_reference_group(). A namespace is locked after loading, so the previous pattern could not work in a package; assigning into the global environment instead is something a package must never do. - All progress reporting moved from
cat()tomessage(), so it can be suppressed, captured, and redirected. Verbosity is set once per session withoptions(cyRAVEN.verbose = )rather than threaded through every call. - The figure raster ceiling is now
options(cyRAVEN.max_raster_px = )instead of a hard-coded constant.
Performance
All three changes below are exact: they compute the same numbers by a cheaper route, and are covered by tests that assert bit-identity rather than closeness.
-
Cofactor derivation is roughly two orders of magnitude faster. The bisection evaluated
asinh()over the entire background vector on each of its 40 iterations, around 264 million transcendental calls per sample at 300,000 events. Becauseasinh(./m)is strictly increasing it maps order statistics to order statistics, and R’s type-7 quantile interpolates between positions that depend only onnandp, never onm. The four order statistics the interquartile range needs are therefore selected once, and each iteration costs four scalarasinh()calls. -
SOM node updates are O(n) per epoch instead of O(nodes x n). The per-node
colSums(X[mapping == j, ])loop is replaced by a singlerowsum()pass, removing about 100 million redundant comparisons per--clusterrun. -
LISI neighbour selection uses a partial sort. Only the identity of the k nearest neighbours affects the Simpson index, so the full
order()of every distance row was unnecessary.
Memory
- Each sample’s raw expression matrix is released as soon as its transformed matrix exists, rather than being held for the whole run. Both were previously alive simultaneously during population scoring, which is exactly where an uncapped run was killed by the OOM reaper. Unstained controls are held back until every sample that could reference them has been scored, and
--keep-exprsstill retains everything.
Fixes
- Long plot subtitles and captions are wrapped to the figure width.
ggplot2does not wrap them, so a long one was drawn on a single line and silently clipped at the device edge, mid-sentence, with no warning. - Heatmap tile labels choose black or white per tile by WCAG contrast ratio against the fill they sit on. A single fixed grey was unreadable on the dark end of a sequential palette, and which tiles were affected depended on the data.
-
--cofactor-from-first-sampleno longer collides with--cofactorthrough R’s partial matching of$on lists.optparseomitsNULL-default options from the list it returns, soopt$cofactormatched the longer name and returnedFALSE, which was then used as a cofactor of zero: every marker becameasinh(x/0) = Infand every sample failed staining QC with a message that pointed nowhere near the cause.
