diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 0000000..1b7c5a6 --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,43 @@ +# Builds the docs on every pull request, so a broken page fails there, and +# publishes them to GitHub Pages on every merge to main. +name: docs + +on: + push: + branches: [main] + pull_request: + +permissions: + contents: read + +jobs: + build: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v7 + - uses: actions/setup-python@v7 + with: + python-version: '3.12' + - run: pip install -r docs/requirements.txt + - run: sphinx-build -W --keep-going docs _site + - if: github.event_name == 'push' + uses: actions/upload-pages-artifact@v5 + with: + path: _site + + deploy: + if: github.event_name == 'push' + needs: build + runs-on: ubuntu-latest + permissions: + pages: write + id-token: write + concurrency: + group: pages + cancel-in-progress: false + environment: + name: github-pages + url: ${{ steps.deployment.outputs.page_url }} + steps: + - id: deployment + uses: actions/deploy-pages@v5 diff --git a/README.md b/README.md index c45b4a2..58a89db 100644 --- a/README.md +++ b/README.md @@ -11,38 +11,28 @@ - **Automated archetype annotation** — contrastive image analysis and standardized reports generated with multimodal language models and user review. - **Morphological program definition** — integration of co-regulated archetypes into interpretable programs of coordinated cellular change. -## Requirements +## Documentation -Most of CHART is Python and arrives with the package. The archetype -analysis is the exception: every step of it but the first and last is an -R script that `chart archetypes` runs through `Rscript`, so that stage -needs R and eight R packages. `environment.yml` describes them: +The [documentation](https://genentech.github.io/CHART/) describes the input data, each stage's usage, outputs and method, the configuration, and running on a cluster. + +## Installation ``` -mamba env create -f environment.yml -conda activate chart-r -Rscript -e 'install.packages("archetypes", repos="https://cloud.r-project.org")' +pip install . ``` -The archetype annotation is optional. It needs two more packages and a key -for the language model, read from the environment: +The archetype analysis also needs R and the packages in `environment.yml`: ``` -pip install ".[annotation]" -export OPENAI_API_KEY="sk-..." +mamba env create -f environment.yml +conda activate chart-r +Rscript -e 'install.packages("archetypes", repos="https://cloud.r-project.org")' ``` -`src/chart/annotation/llm.py` lists the gateway and Azure alternatives. +The optional annotation stage needs `pip install ".[annotation]"` and credentials for a language model in the environment. [Installation](https://genentech.github.io/CHART/installation.html) covers these, scallops for reading images, and remote storage. ## Quickstart -Install the package, plus the environment above if you want the archetype -analysis: - -``` -pip install . -``` - Write a config saying where the input is and where output should go: ```yaml @@ -67,86 +57,4 @@ chart effects --config my_screen.yaml --all-channels chart programs --config my_screen.yaml ``` -`chart effects` tests each gene against the controls on every archetype -of the weights `chart archetypes` left, and writes one table per channel -to `effects/` under `data_output_dir`. Once every channel in -`archetypes.channels` has its table, it corrects all their p-values -together and writes them, with q-values, to `effects/pooled.parquet`. - -`chart programs` finds programs, sets of archetypes that genes move -together, by sparse PCA on the genes and archetypes in -`effects/pooled.parquet` with any q-value below 0.01. It writes to -`programs/` under `data_output_dir`: each program's loadings and each -gene's scores; the genes' modules, from how often they cluster together -across 150 refits on part of them, which takes about 25 minutes; and every -cell's scores, with 20 example cells per program. `--steps` runs only some -of `matrix`, `fit`, `modules` and `cells`. `examples/programs.ipynb` shows -what each table holds and how to read them; Jupyter is not installed with -CHART. - -The last step of `chart archetypes`, `crops`, draws each archetype's most -heavily weighted cells from the images, five to a PNG, under -`archetypes//crops/` in `local_output_dir`. It reads the images -through scallops, from the locations under `preprocessing_bywell`. A -channel that `archetypes.image_channels` does not list, or a machine -without scallops, gets a warning and no crops, or an error if the step -was asked for with `--steps crops`; `example_config.yaml` lists its -settings. - -`chart annotation` describes each archetype from those crops with a -language model, leaving out the archetypes `archetypes_to_remove` lists, -and writes the reports and a dashboard to `annotation//` under -`local_output_dir`. Nothing downstream reads them. - -Wells are read off the input filenames. Name wells -individually with `--well` and channels with `--channel`, and use -`--steps` on every stage but `chart effects` to run only part of a pipeline. -`chart archetypes --dry-run` prints every command the stage would run -without running any of them. - -`chart archetypes` runs one channel at a time. -`archetypes.parallel_channels`, or `--parallel-channels`, runs several -at once, each in its own process and logging to -`archetypes//` under `local_output_dir`; the screen says when -each starts and ends. Memory sets how many fit, not cores: the metacells -step needed about 62 GB for one channel of 875,000 cells, and uses one -core. The other stages use every core the run may use, which under SLURM -is the job's share. `n_jobs` at the top of the config, or `--n-jobs`, -sets a number instead. - -## Configuration - -Every step reads one YAML file, given with `--config`: - -``` -chart preprocessing --config my_screen.yaml --all-wells -``` - -Use `example_config.yaml` as a starting point: copy it and edit the paths. Every key in it except the two output directories is optional, and the values shown are the defaults. - -The `schema` block maps CHART variables to the input screen, accounting for differences in input column names and feature naming conventions: - -```yaml -schema: - cell_column: object_number - guide_column: sgrna_id - gene_column: target_gene - feature_patterns: ['^Nuclei_', '^Cells_', '^Cytoplasm_'] - control_gene: safe_harbour - control_prefix: null - centroid_columns: [Nuclei_Location_Center_Y, Nuclei_Location_Center_X] -``` - -Every column in `feature_patterns` must be numeric. - -Input filenames are configurable in the same file, under -`preprocessing_bywell`. Only `.parquet` input files are accepted. - -```yaml -preprocessing_bywell: - objects_pattern: 'plate1_{well}.objects.parquet' - features_pattern: 'plate1_{well}.features.parquet' -``` - -When one file per well holds objects and features together, set -`premerged: true` and name it with `merged_pattern` instead. \ No newline at end of file +Before `chart effects`, inspect each archetype's crops under `archetypes//crops/` in `local_output_dir`, and list those reflecting segmentation errors or staining artefacts under `archetypes.archetypes_to_remove`, as described in [Curating archetypes](https://genentech.github.io/CHART/archetypes.html#curating-archetypes). `example_config.yaml` lists every setting, and [Configuration](https://genentech.github.io/CHART/configuration.html) explains which matter most. diff --git a/docs/annotation.md b/docs/annotation.md new file mode 100644 index 0000000..d9b2323 --- /dev/null +++ b/docs/annotation.md @@ -0,0 +1,85 @@ +# Annotation + +`chart annotation` describes each archetype from its crops with a multimodal language model. For every archetype it writes a morphological summary, three proposed names and a functional hypothesis, and collects them with the images in one HTML dashboard per channel. The stage is optional: nothing downstream reads its output. + +It reads the grids drawn by the `crops` step of `chart archetypes`, leaving out the archetypes listed under `archetypes.archetypes_to_remove`. The packages and credentials it needs are described under [Installation](installation.md). + +## Usage + +Each channel needs one or more normal archetypes, chosen from the crops, against which every archetype is compared: + +```yaml +annotation: + normal_archetypes: + TOM20: [10] + FAK: [5, 8, 11] +``` + +``` +chart annotation --config my_screen.yaml --all-channels +chart annotation --config my_screen.yaml --channel TOM20 --steps annotate baseline +``` + +| Step | Output | +|---|---| +| `annotate` | An observational description of each of the archetype's first two grids | +| `baseline` | The archetype's first grid compared with the normal archetypes' | +| `synthesize` | A summary from the first grid's description and the baseline note; three names and a functional hypothesis | +| `refine` | Summaries that read alike, contrasted pairwise and rewritten | +| `dashboard` | The final reports beside their grids, in one HTML file | + +The steps run in order, each reading the previous step's output. `--resume` keeps the descriptions `annotate` has already written, so an interrupted run continues where it stopped. `--continue-on-error` carries on with the remaining channels after a failure. `dashboard` on its own needs no model or credentials. + +> [!NOTE] +> Calls are made one at a time. For a channel with $K$ archetypes, `annotate` makes up to $2K$ image calls, `baseline` $K$, `synthesize` $2K$ text calls, and `refine` two image calls per similar pair and one rewrite per archetype in a pair. `refine` compares summaries with the sentence-transformers model `all-MiniLM-L6-v2`, downloaded from Hugging Face on first use; on a machine without internet access, download it beforehand. + +## Outputs + +All files are written to `annotation//` under `local_output_dir`. Archetype numbers are those of the crops, `_Archetype`. + +| File | Contents | +|---|---| +| `archetype_annot/archetype_grid__annotation.md` | Description of one grid | +| `contrastive_baselines/contrastive_baseline_notes.json` | Each archetype against the normal archetypes; `contrastive_baseline_report.md` holds the same notes | +| `archetype_reviews/archetype_synthesis.md` | Summary, names and hypothesis | +| `archetype_reviews/archetype_contrast_vs.md` | Traits distinguishing archetype N from M | +| `archetype_reviews/archetype_final.md` | Final summary, names, hypothesis, and the archetypes it was contrasted with | +| `archetype_reviews/combined_archetype_synthesis.md`, `combined_archetype_final.md` | Every archetype's report in one file | +| `archetype_dashboard.html` | The final reports with their grids embedded, readable anywhere | + +## Configuration + +```yaml +annotation: + reports_dir: annotation/ # under local_output_dir + normal_archetypes: {} # per channel, by archetype number + model: gpt-5.5 + similarity_threshold: 0.08 # cosine distance below which summaries are contrasted + max_pairs_per_archetype: 0 # closest partners kept per archetype; 0 keeps all +``` + +The values shown are the defaults. `normal_archetypes` has to be set for every channel before `baseline` runs. The stage also reads `archetypes.channels`, `archetypes.archetypes_to_remove`, and `archetypes.reports_dir`, under which the crops are. + +## Method + +**Grids.** The crops step draws each archetype's most heavily weighted cells, `archetypes.top_cells` (10 by default), five to a grid: the stain in inverted grey, the cell outlined in red and its nucleus in blue. The first two grids of each archetype are described; only the first enters the later steps. + +**Observation.** The descriptions, the baseline comparisons and the pairwise contrasts are restricted to visible features: intensity, localisation relative to the two outlines, texture, cell and nuclear size and shape, and artefacts. The model is instructed not to infer biology or to assume the stain's usual pattern. The only interpretive call is the one in `synthesize` that proposes names and a hypothesis, from the first grid's description. + +**Baseline.** Each archetype's first grid is sent together with the normal archetypes' first grids in one request, at reduced resolution (longest edge at most 1024 pixels). + +**Refinement.** The summaries are embedded with `all-MiniLM-L6-v2`, and each pair closer than `similarity_threshold` in cosine distance is contrasted in both directions, each call showing both first grids. Every archetype in a pair then has its summary rewritten to include what distinguishes it. With `max_pairs_per_archetype` above 0, only each archetype's closest partners are kept. Archetypes in no pair keep their synthesis summary. + +**Validation.** Every reply must contain the sections the next step reads; a reply without them is requested once more, then the step fails. Empty and truncated replies are errors, and rate-limited calls are retried with increasing delays. A step therefore either writes every archetype's report or stops. + +## Troubleshooting + +| Message | Cause and remedy | +|---|---| +| `No top_cells.parquet in ...; run the crops step of chart archetypes for ...` | The channel has no crops; run `chart archetypes --steps crops` for it | +| `top_cells.parquet lists archetype ..., but ... is not in ...; rerun the crops step` | The grid images are incomplete; rerun the crops step | +| `'annotation.normal_archetypes' has to name archetypes of ... that the crops step drew and 'archetypes.archetypes_to_remove' keeps` | No normal archetype is set for the channel, or one named was removed or not drawn | +| `Missing langchain_openai, sentence_transformers. ...` | The annotation extra is not installed | +| `No LLM provider is configured ...` | No credentials in the environment; see [Installation](installation.md) | +| `... model returned no text ...`, `... reply was cut off at the token ceiling ...` | The model spent its token budget; rerun the step, with `--resume` for `annotate` | +| `... reply is missing required section header(s) ...` | The model twice ignored the output format; rerun the step | diff --git a/docs/archetypes.md b/docs/archetypes.md new file mode 100644 index 0000000..798b6da --- /dev/null +++ b/docs/archetypes.md @@ -0,0 +1,119 @@ +# Archetypes + +`chart archetypes` describes each channel's single-cell morphology as a mixture of archetypes, extreme morphological states between which every cell lies. It builds a PCA space from the channel's features, groups the cells into metacells, fits archetypal analysis to the metacells, and weights every cell on each archetype. The weights are the input to `chart effects` and `chart programs`; the last step crops each archetype's most heavily weighted cells from the images, for inspection and for `chart annotation`. + +## Usage + +``` +chart archetypes --config my_screen.yaml --all-channels +chart archetypes --config my_screen.yaml --channel TOM20 --steps selectk weights crops +``` + +| Step | Runs | Output | +|---|---|---| +| `pca` | Python | The channel's PCA space, centred on the controls | +| `npc` | R | The number of principal components kept | +| `metacells` | R, per parameter pair | Metacells, and their coordinates | +| `evaluate` | R | A comparison of the metacell partitions; optional, as nothing downstream reads it | +| `archetypes` | R, per parameter pair | Repeated archetype fits for each number of archetypes from 6 to 20 | +| `selectk` | R | The number of archetypes, K | +| `weights` | R | Every cell's weight on each of the K archetypes | +| `crops` | Python | Each archetype's top cells, cropped from the images | + +The R steps need the environment described under [Installation](installation.md); the stage checks that R can load every package before the first step. `--dry-run` prints every command the stage would run, and needs no R. `--continue-on-error` carries on with the remaining steps and channels after a failure. + +Channels run one at a time. `archetypes.parallel_channels`, or `--parallel-channels`, runs several at once, each in its own process and logging to `archetypes//` under `local_output_dir`. Memory, not cores, limits how many fit: the `metacells` step needed about 62 GB for one channel of 875,000 cells. Channels can also run as separate jobs with `--channel`. + +> [!IMPORTANT] +> `selectk` chooses K automatically and says so in the log. Check its choice against `evalK_recovery_plot.png` before going further; to overrule it, set `archetypes.k` for the channel and rerun `--steps selectk weights crops`. + +`crops` needs scallops, the image locations described under [Your data](data.md), and the channel's position in the images under `archetypes.image_channels`. Without them the step is skipped with a warning, or stops the run under `--steps crops`. Once the crops are drawn, curate the archetypes, as described under Curating archetypes below, before running `chart effects`. + +## Outputs + +Each channel's data are written to `archetypes//` under `data_output_dir`, and its reports to `archetypes//` under `local_output_dir`. File names carry the parameters, `gamma`, `knn`, `npc` and `k`. + +| File | Under | Contents | +|---|---|---| +| `cell_pca_centered.parquet` | data | The cells' principal components, with `Well`, `Label`, `Guide` and `Gene` | +| `npc.txt`, `k.txt` | data | The number of components and of archetypes in use | +| `supercells-_gamma…_knn…_npc….rds`, `.parquet` | data | Metacell membership, and metacell coordinates | +| `archetype_robust_runs_gamma…_knn…_npc….rds` | data | Every archetype fit, with its errors and stability measures | +| `robust_archweights-_k…_gamma…_knn…_npc….parquet` | data | One row per cell, a column per archetype, `Archetype1` to `ArchetypeK`, read by `chart effects` and `chart programs` | +| `npc_elbow.pdf` | local | The variance curve and the components kept | +| `space_evals.csv` | local | Size, dispersion and separation of each partition into metacells | +| `evalK_recovery_curves.csv`, `evalK_recovery_plot.png` | local | Reconstruction error against K, and why any K was excluded | +| `robust_archweights-….pdf` | local | The distribution of each archetype's weights over cells | +| `crops/archetype__grid_.png` | local | Each archetype's top cells, five to an image | +| `crops/top_cells.parquet` | local | The cells shown, with their weights and bounding boxes | + +## Curating archetypes + +Some archetypes capture technical variation rather than morphology. These should be identified from the crops and excluded before testing. Each grid shows five of an archetype's most heavily weighted cells, `top_cells` in all: the channel's stain in inverted grey, the cell outlined in red and its nucleus in blue, and above each the weight, well and cell label. Archetypes to exclude are typically those whose top cells reflect segmentation errors, for instance outlines enclosing several cells or a fragment of one, or staining artefacts such as debris, precipitates or out-of-focus regions. `top_cells.parquet` gives each cell's well, which shows whether an archetype is drawn from a single well. + +Excluded archetypes are listed under `archetypes.archetypes_to_remove`, grouped under category names of your choice; only the archetype names, `_ArchetypeN`, are read: + +```yaml +archetypes: + archetypes_to_remove: + seg_err: [TOM20_Archetype1] + staining_artifacts: [TOM20_Archetype4] +``` + +`chart effects` does not test them, and since all tests are corrected together, removing one also changes the q-values of the rest; `chart annotation` does not describe them. The numbering belongs to one fit: changing K, `npc` or the selected parameter pair renumbers the archetypes, and the list has to be made again. + +## Configuration + +```yaml +archetypes: + channels: [TOM20, Golgin97] + all_channels: [DAPI1, TOM20, Golgin97, CALR, EEA1] + outputs_dir: archetypes/ # under data_output_dir + reports_dir: archetypes/ # under local_output_dir + params: [[30, 10]] # gamma and k.knn pairs + # selected_params: [30, 10] # the pair selectk and weights use; by default the first + # npc: {TOM20: 14} # by default the elbow of the variance curve + # k: {TOM20: 12} # by default the knee of the recovery curve + coverage_sample_size: 12000 + eval_sample_size: 8000 + fit_sample_size: 12000 + image_channels: {TOM20: 8, Golgin97: 9} # position among the image channels, from 0 + top_cells: 10 + # display_percentiles: {TOM20: [0, 99.9]} # by default 0 and 99.5 + # display_ranges: {TOM20: {A1: [105, 4000]}} + parallel_channels: 1 +``` + +Where a key has a default, the value shown is the default and that of the published analysis. `channels` are the channels analysed; `all_channels`, every channel in the feature table, is required by `pca`: a channel's features are those that name no other channel, so a channel left out of the list keeps its features in every other channel's space. `display_percentiles` sets, per channel, the pixel percentiles between which the stain is shown, and `display_ranges` fixes the range for particular wells. The stage also reads `schema`, the guide-filtered cell table from [Preprocessing](preprocessing.md), and, for `crops`, the combined objects table and the image locations. + +## Method + +**Channel space.** A channel's features are the columns of the guide-filtered cell table that name no other channel in `all_channels` and measure no shape. They are standardised, reduced to the principal components covering 90% of the variance, and centred on the centroid of the control cells. + +**Components.** `npc` is the elbow of the log variance per component over the first 30 components, the point furthest from the chord joining the curve's ends, held between 10 and 30. + +**Metacells.** SuperCell groups the cells on a k-nearest-neighbour graph of their first `npc` components, with `k.knn` neighbours and a graining level `gamma`, the ratio of cells to metacells. A metacell's coordinates are the mean of its cells'. With several pairs under `params`, `evaluate` tabulates for each partition the metacell sizes, the dispersion within metacells, and the distance between neighbouring metacells relative to their radius, as a basis for `selected_params`. + +**Archetype fits.** On the standardised metacell coordinates, archetypal analysis is fitted eight times for each K from 6 to 20, each fit the best of 10 random starts. Each fit uses `fit_sample_size` metacells: 70% from a coverage set of `coverage_sample_size` metacells chosen by farthest-point sampling, the rest drawn with probability proportional to rarity, the distance to the 30th nearest metacell, raised to 1.5. Every fit is scored by the residual sum of squares of a shared evaluation sample of `eval_sample_size` metacells, half from the coverage set and half rarity-weighted, projected onto the simplex of its archetypes. For each K, the fit with the lowest error is the reference; drift is the 90th percentile of the distances between its archetypes and those of the other fits, matched by Hungarian assignment, and support is the smallest number of evaluation metacells weighted above 0.3 on any one archetype. + +**Choice of K.** For each K, the evaluation metacells' weights on the reference archetypes reconstruct their coordinates in all components, with the archetypes' positions in that space fitted by least squares; the root-mean-square error, discarding the worst 10% of metacells, forms the recovery curve. K values with an unsupported archetype, or with drift above twice the median drift, are excluded, and K is the knee of the remaining curve with both axes scaled to [0, 1]. The error decreases with K, so it has no minimum to choose. + +**Weights.** K archetypes are fitted once to all metacells, unstandardised, as the best of 20 random starts with a fixed seed, so a rerun on the same metacells gives the same archetypes. Each cell's weights are the convex coefficients, non-negative and summing to one, that best reconstruct its first `npc` components from the archetypes. + +**Crops.** Each archetype's `top_cells` most heavily weighted cells are cropped at their bounding box with an 8-pixel margin. A well's display range lies between two percentiles of its pixels, sampled in five 5000-pixel patches, one central and one per quadrant, inside a 10% margin. + +## Troubleshooting + +| Message | Cause and remedy | +|---|---| +| `Rscript is not on PATH ...` or `R cannot load ...` | Activate the R environment described under [Installation](installation.md) | +| `The pca step needs every channel the feature table holds ...` | Set `archetypes.all_channels` | +| `None of the ... columns belong to ... alone ...` | Every feature names another channel or measures a shape; check the channel names against the feature names | +| `There are only ... metacells, so every fit uses all of them ...` (warning) | The fits differ only in their starts; lower `coverage_sample_size` and `fit_sample_size` to about 40% of the metacells | +| `Only ... of ... K values survive the drift and support checks ...` | `evalK_recovery_curves.csv` gives the reason for each exclusion; set `archetypes.k` | +| `'archetypes.npc' gives ... but its npc.txt holds ...` | The weights were built with another `npc`; rerun from the `npc` step | +| `No K given and no k.txt in ...` or `No npc.txt ...` | Run the `selectk` or `npc` step first, or set `archetypes.k` | +| `The process running ... stopped before it finished ...` | Usually out of memory; lower `parallel_channels` | +| `Skipping the crops for ...: ...` (warning) | `archetypes.image_channels` lacks the channel, or scallops is not installed | +| `No '...' is configured for well ... under 'preprocessing_bywell' ...` | Set `pheno_dir`, `cell_labels_dir` and `nuclei_labels_dir`, described under [Your data](data.md) | +| `... cells shown for ... have no outline in the cell label image ...` (warning) | The label images are not the segmentation the objects table came from, or `cell_column` does not hold the label values | diff --git a/docs/cluster.md b/docs/cluster.md new file mode 100644 index 0000000..33f8b91 --- /dev/null +++ b/docs/cluster.md @@ -0,0 +1,58 @@ +# Running on a cluster + +Each stage is a single command, so on a cluster it becomes one batch job, or one per well or channel. This page describes how the stages divide into jobs, what they use of a job's cores and memory, and what the compute nodes need. The examples use SLURM. + +`data_output_dir` and `local_output_dir` have to be reachable from every node that runs a job, on a shared filesystem or on [remote storage](installation.md). + +## Dividing a run into jobs + +| Order | Command, each with `--config` | Jobs | +|---|---|---| +| 1 | `chart preprocessing --well --steps qc filter` | One per well | +| 2 | `chart preprocessing --all-wells --steps combine normalize outliers guides` | One, after every job of 1 | +| 3 | `chart archetypes --channel ` | One per channel | +| 4 | `chart annotation --channel ` | Optional; one per channel, after 3 | +| 5 | `chart effects --channel ` | One per channel, after 3 and curation | +| 6 | `chart programs` | One, after every job of 5 | + +Each row can also run as a single job, with `--all-wells` or `--all-channels`. Between 3 and 5 the archetypes are curated by hand, as described under [Archetypes](archetypes.md). Each effects job pools the q-values once every channel has a table, so the last to finish writes `pooled.parquet`. Every command exits with a non-zero status when a step fails, so dependent jobs are held back. + +A job array runs the per-well steps, and a dependency starts the across-wells steps once all of them have succeeded: + +```bash +#!/bin/bash +#SBATCH --array=0-5 +wells=(A1 A2 A3 B1 B2 B3) +chart preprocessing --config my_screen.yaml --well "${wells[$SLURM_ARRAY_TASK_ID]}" --steps qc filter +``` + +``` +jobid=$(sbatch --parsable bywell.sh) +sbatch --dependency=afterok:$jobid allwells.sh +``` + +`chart archetypes --dry-run` prints every command a job would run, without running any and without R, which checks a configuration before submission. + +## Cores and memory + +`n_jobs: -1`, the default, uses the cores the process may run on, which under SLURM are the job's allocation rather than the node's. + +| Stage | Parallelised over | Set by | +|---|---|---| +| `chart preprocessing` | Threads in `normalize`, and the local outlier factor in `guides` | `n_jobs`, `--n-jobs` | +| `chart archetypes` | Channels, as separate processes; the R scripts are not parallelised | `archetypes.parallel_channels`, `--parallel-channels` | +| `chart annotation` | Nothing; one model call at a time | | +| `chart effects` | Perturbations, over worker processes | `effects.n_jobs`, `n_jobs`, `--n-jobs` | +| `chart programs` | The sparse PCA | `n_jobs`, `--n-jobs` | + +Memory is the limiting resource for `chart archetypes`. For one channel of 875,000 cells, the `metacells` step took about 5 hours and held 62 GB, and the `weights` step took 10 minutes. A job running several channels at once with `parallel_channels` needs that memory for each, and shares its cores evenly between them. `normalize` holds two copies of the screen's feature matrix as 64-bit floats, about 1.2 GB each for 300,000 cells by 500 features. The `modules` step of `chart programs` took about 20 minutes on 96 cores for 471 perturbations and 188 archetypes. + +## Compute nodes + +`chart archetypes` needs `Rscript` on `PATH`, so the job script activates the R environment described under [Installation](installation.md); the stage checks that R can load every package before its first step. The `pheno` and `pheno_sbs` components and the `crops` step need scallops on the node. + +`chart annotation` needs network access to the model's endpoint, and the credentials in environment variables. SLURM passes the submitting shell's environment to the job by default, so they can be exported before `sbatch`. The `refine` step loads the all-MiniLM-L6-v2 model from Hugging Face on first use; where compute nodes have no internet access, download it beforehand into a Hugging Face cache (`HF_HOME`) that the nodes share. + +## Logs + +`chart preprocessing` writes its log to `preprocessing/bywell/qc_reports/` under `local_output_dir`, one file per run named after its well, the number of wells, or `all_wells`, and a timestamp, as well as to standard output. The other stages log to standard error, which SLURM writes to the job's output file, together with the output of the R scripts. With `parallel_channels`, each channel instead logs to `archetypes_log_.txt` in `archetypes//` under `local_output_dir`, and the job's output records when each channel starts and ends. diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..1228b91 --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,15 @@ +project = 'CHART' +author = 'Genentech, Inc.' +copyright = '2025, Genentech, Inc.' + +extensions = ['myst_parser'] +myst_enable_extensions = ['alert', 'colon_fence', 'dollarmath'] +exclude_patterns = ['_build'] +highlight_language = 'none' + +html_theme = 'pydata_sphinx_theme' +html_title = 'CHART' +html_theme_options = {'github_url': 'https://github.com/Genentech/CHART', + 'header_links_before_dropdown': 9} +html_sidebars = {'**': []} +html_show_sourcelink = False diff --git a/docs/configuration.md b/docs/configuration.md new file mode 100644 index 0000000..92045a8 --- /dev/null +++ b/docs/configuration.md @@ -0,0 +1,50 @@ +# Configuration + +Every stage reads one YAML file, given with `--config`. [`example_config.yaml`](https://github.com/Genentech/CHART/blob/main/example_config.yaml), reproduced at the end of this page, is the starting point: copy it, set the two output directories, and adjust what the screen requires. Each key is commented there, and the values shown are the defaults, except the channel lists, which are examples. + +## How the file is read + +| Section | Read by | +|---|---| +| top level: `data_output_dir`, `local_output_dir`, `n_jobs` | every stage | +| `schema` | every stage | +| `preprocessing_bywell` | `chart preprocessing`; the image locations also by the `crops` step | +| `preprocessing_allwells`, `guide_filtering` | `chart preprocessing` | +| `archetypes` | `chart archetypes`, and the channels and removed archetypes by every later stage | +| `annotation`, `effects`, `programs` | the stage of the same name | + +Only the two output directories are required. Tables are written under `data_output_dir`, and reports, plots and filter lists under `local_output_dir`. A relative path in the file is taken under the corresponding output directory, and image directories under `scallops_dir` where it is set; absolute paths and URLs such as `s3://` are used as given. A key CHART does not recognise is reported with a warning, `Unknown config key '...' in ...; it will be ignored`, so a misspelt key does not pass silently. Command-line options such as `--n-jobs` and `--parallel-channels` override the corresponding keys for one run. + +## Settings by importance + +**Required.** `data_output_dir` and `local_output_dir`. + +**Describing the screen.** These depend on the data, and their defaults fit only the reference dataset. They are described under [Your data](data.md), and the channels under [Archetypes](archetypes.md). + +- `schema`: the identifier columns, `feature_patterns`, the controls `control_gene` and `control_prefix`, and the coordinate, bounding-box and area columns that quality control and the crops read +- `preprocessing_bywell`: the input location and file names, `merged_dir`, `objects_pattern`, `features_pattern`, or `premerged` and `merged_pattern`; and the image locations under `global` and `wells` +- `archetypes`: `channels`, `all_channels` and `image_channels` + +**Decided during the analysis.** These are set after inspecting the outputs of an earlier step. + +- `archetypes.k`, and `archetypes.npc`, after checking the recovery curve and the variance elbow ([Archetypes](archetypes.md)) +- `archetypes.archetypes_to_remove`, after inspecting the crops ([Archetypes](archetypes.md)) +- `annotation.normal_archetypes`, before annotating ([Annotation](annotation.md)) + +**Analysis parameters.** The defaults are those of the published analysis. Changing them changes the results, so they are best changed deliberately and recorded. + +- `preprocessing_bywell.thresholds` and `guide_filtering`: the quality-control and guide-filtering thresholds ([Preprocessing](preprocessing.md)) +- `archetypes.params`: the metacell graining, `gamma` and `k.knn` ([Archetypes](archetypes.md)) +- `effects.n_perm` and `effects.trimmed_alpha` ([Effects](effects.md)) +- `programs.k`, `programs.n_modules` and `programs.qval_threshold` ([Programs](programs.md)) +- `annotation.model` ([Annotation](annotation.md)) + +**Resources.** `n_jobs`, `effects.n_jobs` and `archetypes.parallel_channels` set how many cores and channels a run uses, and do not change the results ([Running on a cluster](cluster.md)). + +**Rarely changed.** The remaining keys have defaults that suit most screens: the output subdirectories (`outputs_dir`, `reports_dir`), `drop_unassigned` and `exclude_patterns`, the normalisation and outlier settings in `preprocessing_allwells`, the archetype sample sizes and `selected_params`, the crop display (`top_cells`, `display_percentiles`, `display_ranges`), the seeds, and the annotation's `similarity_threshold` and `max_pairs_per_archetype`. + +## The example file + +```{literalinclude} ../example_config.yaml +:language: yaml +``` diff --git a/docs/data.md b/docs/data.md new file mode 100644 index 0000000..d69c190 --- /dev/null +++ b/docs/data.md @@ -0,0 +1,96 @@ +# Your data + +CHART starts from single-cell image feature tables (e.g. generated by CellProfiler or scallops): one row per segmented cell, with its guide assignment and its image-derived features. Images are needed only by two quality-control components and by the archetype crops. + +## Tables + +Each well contributes an objects table and a features table, both parquet. Both hold one row per cell and three identifier columns: the cell's label, unique within the well; its guide; and the perturbation the guide targets. The two tables are joined on these, and cells missing from the features table are dropped with a warning. The objects table carries the geometry that quality control and the crops read (areas, centroids, bounding boxes); the features table carries the measurements. A screen that writes one table per well holding both sets `premerged: true` under `preprocessing_bywell`. + +Cells without a guide or perturbation assignment are kept unless `preprocessing_bywell.drop_unassigned` is `true`. + +## File names and wells + +The tables are read from `preprocessing_bywell.merged_dir`, by default `preprocessing/bywell/` under `data_output_dir`; an absolute path or a URL points elsewhere. + +```yaml +preprocessing_bywell: + merged_dir: /data/screen1/objects + objects_pattern: 'plate1_{well}.objects.parquet' + features_pattern: 'plate1_{well}.features.parquet' +``` + +Without patterns the files are `{well}-objects.parquet` and `{well}-features.parquet`, and a premerged table is `{well}.parquet` unless `merged_pattern` names it. `{well}` is the only placeholder. + +The wells are the values `{well}` takes among the objects files, or the premerged ones. That name identifies the well in every output and must be the key of its images. Once wells are listed under `preprocessing_bywell.wells`, which per-well image settings require, only those are processed. + +## Columns + +The `schema` block says which columns hold each role. The defaults are the published screen's names. + +| Key | Role | Default | +|---|---|---| +| `cell_column` | cell label, unique within a well | `label` | +| `guide_column` | guide | `sgRNA` | +| `gene_column` | perturbation the guide targets | `gene_symbol` | +| `feature_patterns` | regular expressions selecting the feature columns | `[cell_, nuclei_, cytosol_]` | +| `control_gene` | perturbation marking control cells | `NTC` | +| `control_prefix` | prefix marking further control perturbations; `null` for none | `OR` | +| `centroid_columns` | cell coordinates, y then x, for the quality-control and outlier position plots | `[nuclei_centroid-0, nuclei_centroid-1]` | +| `bbox_columns` | bounding box for the crops: first row, first column, then one past the last of each | `[cell_bbox-0, cell_bbox-1, cell_bbox-2, cell_bbox-3]` | +| `nuclei_area_column` | nuclear area, for the `plots` and `segmentation` quality-control components | `nuclei_area` | +| `cell_area_column` | cell area, for `plots` | `cell_area` | +| `cytosol_area_column` | cytoplasmic area, for `segmentation` | `cytosol_area` | +| `nuclei_bbox_columns` | nuclear bounding box, ordered as `bbox_columns`, for `pheno` | `[nuclei_bbox-0, nuclei_bbox-1, nuclei_bbox-2, nuclei_bbox-3]` | + +A quality-control component whose columns are missing stops the run, naming the key; `--components` leaves it out. + +For a CellProfiler screen with its own identifiers: + +```yaml +schema: + cell_column: object_number + guide_column: sgrna_id + gene_column: target_gene + feature_patterns: ['^Nuclei_', '^Cells_', '^Cytoplasm_'] + control_gene: safe_harbour + control_prefix: null + centroid_columns: [Nuclei_Location_Center_Y, Nuclei_Location_Center_X] + nuclei_area_column: Nuclei_AreaShape_Area + cell_area_column: Cells_AreaShape_Area + cytosol_area_column: Cytoplasm_AreaShape_Area +``` + +### Features + +`feature_patterns` match anywhere in a column name: `nuclei_` matches `cell_nuclei_ratio`, `^nuclei_` only names that start with it. Of the matching columns, those matching `preprocessing_bywell.exclude_patterns` are dropped. These are case-insensitive and by default `bbox`, `centroid`, `location`, `centermass`, `to_origin`, `_x`, `_y` and `orientation`, which mark position and orientation rather than measurements. Every remaining column must be numeric; a text column stops the run at the `normalize` step. + +### Channels + +The archetype analysis works in one feature space per channel, so feature names must contain their channel, as in the published screen's `nuclei_mean_TOM20`. `archetypes.all_channels` lists every channel in the table. A channel's space is every feature naming no other channel, less shape measures such as eccentricity, solidity and perimeter. Consequently: + +- colocalisation features, such as `cell_corr_TOM20_CALR`, belong to neither channel; +- features naming no channel, such as `cell_area`, enter every channel's space; +- a channel missing from `all_channels` stays in every other channel's space; +- no channel name may contain another, as `DAPI` is contained in `DAPI1`. + +### Controls + +Control cells are those whose perturbation equals `control_gene` or starts with `control_prefix`: in the published screen, non-targeting guides (`NTC`) and guides against olfactory receptors (`OR…`), which these cells do not express. They are the reference for the guide filter, the origin of each channel's PCA space and the comparison group of every test in `chart effects`. A step that finds none stops. + +## Images + +Images are read through scallops. Their locations go under `preprocessing_bywell`, in `global` for every well or under `wells` per well; a relative path is taken under `scallops_dir`. + +```yaml +preprocessing_bywell: + global: + cell_labels_dir: /data/screen1/cell-seg.zarr + nuclei_labels_dir: /data/screen1/nuclei-seg.zarr + wells: + A1: + pheno_dir: /data/screen1/pheno/A1 + A2: + pheno_dir: /data/screen1/pheno/A2 +``` + +The crops read each well's phenotype image, `images/` in `pheno_dir`, and its segmentation, `labels/-cell` and `labels/-nuclei` in the two label zarrs. `archetypes.image_channels` gives each stain's position among the image's channels. A cell is outlined where the label image holds its label, so `cell_column` must hold the segmentation's label values. The `pheno` and `pheno_sbs` components read further keys, described under [Preprocessing](preprocessing.md). diff --git a/docs/effects.md b/docs/effects.md new file mode 100644 index 0000000..70c1401 --- /dev/null +++ b/docs/effects.md @@ -0,0 +1,86 @@ +# Effects + +`chart effects` tests every perturbation against the controls on every archetype of each channel, using the archetype weights from `chart archetypes`. Once every channel has been tested, the p-values of all channels are corrected together, and the result, `effects/pooled.parquet`, is the input to `chart programs`. + +## Usage + +> [!IMPORTANT] +> Run `chart archetypes` through its last step, `crops`, before testing. Inspect each channel's crops, under `archetypes//crops/` in `local_output_dir`, and list the archetypes that reflect segmentation errors or staining artefacts under `archetypes.archetypes_to_remove`. These archetypes are not tested, and since all tests are corrected together, removing one also changes the q-values of the rest. + +```yaml +archetypes: + archetypes_to_remove: + seg_err: [TOM20_Archetype1] + staining_artifacts: [TOM20_Archetype4] +``` + +The category names are free; only the archetype names, `_ArchetypeN`, are read. Then: + +``` +chart effects --config my_screen.yaml --all-channels +chart effects --config my_screen.yaml --channel TOM20 CALR +``` + +Channels can be tested in separate runs, for instance as separate cluster jobs. Each run ends by pooling, which happens only once every channel in `archetypes.channels` has a table, from that run or an earlier one; until then it is skipped with a warning. A channel that fails in a run blocks pooling in that run, since any table of its on disk is from an earlier one. `--continue-on-error` carries on with the remaining channels after a failure. + +Perturbations are distributed over worker processes, by default one per core available to the run, which under SLURM is the job's allocation. `effects.n_jobs` sets their number for this stage, the top-level `n_jobs` for every stage, and `--n-jobs` overrides both. + +## Outputs + +All tables are written to `effects/` under `data_output_dir`. + +| File | Rows | Contents | +|---|---|---| +| `.parquet` | perturbation × archetype | effect size, p-value and cell counts | +| `.json` | | weights read, test settings, archetypes removed, perturbations not tested | +| `pooled.parquet` | perturbation × archetype | every channel's rows, with q-values | +| `pooled.json` | | correction method, number of tests, tables pooled | + +| Column | Contents | +|---|---| +| `Gene` | the perturbation | +| `Archetype` | `_ArchetypeN` | +| `w1` | Wasserstein-1 distance between the perturbation's and the controls' distributions of the archetype's CLR-transformed weight | +| `signed_w1` | `w1`, positive when the perturbation's median exceeds the controls' | +| `w1_trimmed` | `w1` after discarding the outer `trimmed_alpha` of each distribution at both ends | +| `pval` | permutation p-value for `w1` | +| `qval` | Benjamini–Hochberg q-value, in `pooled.parquet` only | +| `n_gene`, `n_ctrl` | perturbation and control cells tested | +| `n_wells` | wells in which the perturbation appears | + +`w1` and `signed_w1` are in centred log-ratio (CLR) units. `chart programs` uses `signed_w1` and `qval`; `w1_trimmed` indicates how much of an effect is carried by the tails. + +## Configuration + +```yaml +effects: + outputs_dir: effects/ + n_perm: 1000 # permutations per test + trimmed_alpha: 0.02 # tail fraction discarded for w1_trimmed + seed: 42 + # n_jobs: -1 # worker processes; the top-level n_jobs unless set +``` + +The values shown are the defaults, and those of the published analysis. The stage also reads `archetypes.channels`, `archetypes.archetypes_to_remove`, and the controls defined by `schema.control_gene` and `schema.control_prefix`. + +A permutation p-value cannot fall below $1/(n_\mathrm{perm} + 1)$, and since Benjamini–Hochberg never lowers a p-value, neither can a q-value. With 1000 permutations the floor is about $10^{-3}$, which leaves room below the default `programs.qval_threshold` of 0.01; a stricter threshold needs more permutations. + +## Method + +**Log-ratio transform.** Each cell's weights in a channel are smoothed away from zero, $(w + \alpha)/(1 + K\alpha)$ with $\alpha = 10^{-6}$ and $K$ the channel's number of archetypes, log-transformed, and centred on the cell's mean log weight. All of the channel's archetypes enter the transform, including those removed from testing. + +**Test.** Controls are the cells whose perturbation is `schema.control_gene` or starts with `schema.control_prefix`. Each perturbation is compared with the control cells pooled from the wells in which it appears. For each archetype, `w1` is the exact one-dimensional Wasserstein-1 distance between the two empirical distributions of the archetype's CLR value, and `signed_w1` carries the sign of the difference in medians. The p-value is $(b + 1)/(n_\mathrm{perm} + 1)$, where $b$ is the number of permutations, shuffling perturbation and control labels within each well, whose distance is at least the observed one. The distance responds to any difference between the two distributions, not only a shift in location, and in either direction. Permutations are seeded from `seed`, so a rerun on the same inputs reproduces the p-values exactly. + +**Correction.** The p-values of every channel in `archetypes.channels` are corrected together by Benjamini–Hochberg, over every perturbation × archetype test. + +## Troubleshooting + +| Message | Cause and remedy | +|---|---| +| `No archetype weights for ... at ...; run the weights step` | `chart archetypes` has not produced weights for that channel with the configured `k` and parameters | +| `None of the ... cells for ... carry a control guide` | No perturbation matches `schema.control_gene` or `schema.control_prefix` | +| `'archetypes.archetypes_to_remove' names ..., but ... has only ...` | An archetype listed for removal does not exist in that channel | +| `... genes have no controls in their wells and were not tested` (warning) | Those perturbations appear only in wells without control cells; they are listed under `skipped_genes` in `.json` | +| `Nothing in 'archetypes.archetypes_to_remove' belongs to ...` (warning) | No archetype of that channel is excluded; inspect its crops and list artefactual ones before testing | +| `Not pooling the q-values: ...` (warning) | A channel in `archetypes.channels` has no table yet, or failed in this run; test it, and pooling follows | +| `The pooled q-values cover only the channels in 'archetypes.channels' ...` (warning) | A channel was tested that `archetypes.channels` does not list; its table is written but not pooled | diff --git a/docs/index.md b/docs/index.md new file mode 100644 index 0000000..23d5a5d --- /dev/null +++ b/docs/index.md @@ -0,0 +1,26 @@ +# CHART + +**CHART** (**Cellular High-content imaging Archetype Response Toolkit**) translates high-content imaging features into interpretable morphological phenotypes. CHART identifies marker-specific cellular archetypes, annotates representative images using multimodal large language models, and integrates co-regulated archetypes into morphological programs that describe coordinated cellular responses to perturbation. + +![Overview of the CHART framework](../images/chart_diagram.png) + +- **OPS preprocessing** — quality control, filtering, normalization, outlier detection, guide filtering and single-cell/aggregated embeddings. +- **Compartment archetype analysis** — identification of distinctive reference morphological states for each imaged marker. +- **Automated archetype annotation** — contrastive image analysis and standardized reports generated with multimodal language models and user review. +- **Morphological program definition** — integration of co-regulated archetypes into interpretable programs of coordinated cellular change. + +Start with [Installation](installation.md) and [Your data](data.md); a quickstart is in the [README](https://github.com/Genentech/CHART#readme). + +:::{toctree} +:hidden: + +installation +data +preprocessing +archetypes +annotation +effects +programs +Cluster +configuration +::: diff --git a/docs/installation.md b/docs/installation.md new file mode 100644 index 0000000..f689c49 --- /dev/null +++ b/docs/installation.md @@ -0,0 +1,56 @@ +# Installation + +CHART is a Python package (Python ≥ 3.9) that installs the `chart` command. Some parts of the pipeline need more than the package itself. + +| Component | Needed for | Install | +|---|---|---| +| CHART | every stage | `pip install .` | +| R 4.5 and its packages | `chart archetypes`, every step but `pca` and `crops` | `environment.yml` | +| scallops | the `crops` step of `chart archetypes`; the `pheno` and `pheno_sbs` QC components | `pip install scallops` | +| annotation extra | `chart annotation` | `pip install ".[annotation]"` | +| fsspec and a filesystem driver | inputs or outputs on remote storage | `pip install fsspec s3fs` for S3 | + +## CHART + +``` +git clone https://github.com/Genentech/CHART.git +cd CHART +pip install . +``` + +## R + +All steps of `chart archetypes` except `pca` and `crops` are R scripts, which the stage runs with `Rscript`. `environment.yml` describes a conda environment, `chart-r`, with R 4.5 and the packages the scripts load. The `archetypes` package is on CRAN but not conda-forge, so it is installed separately: + +``` +mamba env create -f environment.yml +conda activate chart-r +Rscript -e 'install.packages("archetypes", repos="https://cloud.r-project.org")' +``` + +`Rscript` must be on `PATH` when `chart archetypes` runs, for instance with `chart-r` activated alongside the Python environment holding CHART. Before its first step, the stage checks that R can load every package the scripts need, and names any it cannot. The compiled packages, `r-supercell` in particular, are built for one R minor version, which is why the environment fixes R to 4.5. + +## scallops + +[scallops](https://github.com/Genentech/scallops) reads the registered images. Without it, CHART skips the components that need images, with a warning: the crops of archetypal cells, and the phenotype-correlation QC components. An explicit request for one of them, `--steps crops` or `--components pheno`, fails instead. scallops requires Python ≥ 3.12. + +## Annotation + +`chart annotation` describes archetypes from their crops with a multimodal language model. Its packages are an optional extra: + +``` +pip install ".[annotation]" +``` + +The model is set by `annotation.model` in the config, `gpt-5.5` by default. Credentials are read from the environment only. Two providers are supported: + +| Provider | Variables | +|---|---| +| OpenAI-compatible gateway, such as LiteLLM | `LITELLM_API_KEY` and `LITELLM_BASE_URL`; `OPENAI_API_KEY` and `OPENAI_BASE_URL` are also accepted, and `OPENAI_API_KEY` alone addresses OpenAI itself | +| Azure OpenAI | `AZURE_OPENAI_API_KEY` and `AZURE_OPENAI_ENDPOINT`; optionally `AZURE_OPENAI_API_VERSION`, and `AZURE_DEPLOYMENT_MAP` where deployment names differ from model names, e.g. `'{"gpt-5.5": "gpt-5-5-deploy"}'` | + +`LLM_PROVIDER` selects `gateway`, `azure` or `auto`; `auto`, the default, uses the gateway when its variables are set and Azure otherwise. + +## Remote storage + +Paths in the config may point to remote storage, such as `s3://` URLs. Reading and writing them needs `fsspec` and the driver for the filesystem, `s3fs` for S3, which are not installed with CHART. diff --git a/docs/preprocessing.md b/docs/preprocessing.md new file mode 100644 index 0000000..f694fc5 --- /dev/null +++ b/docs/preprocessing.md @@ -0,0 +1,104 @@ +# Preprocessing + +`chart preprocessing` turns the per-well single-cell tables into the screen-wide cell table that the archetype analysis starts from. Within each well it removes cells that fail quality control; across wells it normalises the features, removes outlier cells within each perturbation, and removes guides that do not behave like the rest of their perturbation or, for controls, like the other controls. Its inputs are described under [Your data](data.md). + +## Usage + +``` +chart preprocessing --config my_screen.yaml --all-wells +chart preprocessing --config my_screen.yaml --well A1 --steps qc filter +``` + +| Step | Runs | Output | +|---|---|---| +| `qc` | per well | Lists of cells failing quality control | +| `filter` | per well | The well's tables without those cells, reduced to the feature columns | +| `combine` | across wells | One objects table and one feature table for the screen | +| `normalize` | across wells | Features quantile-normalised across wells | +| `outliers` | across wells | Outlier cells removed within each perturbation | +| `guides` | across wells | Unreliable guides removed; cell, guide and perturbation tables | + +The steps run in order. With `--well`, only the per-well steps run unless the across-wells steps are named in `--steps`, since those write one table for the whole screen. The per-well steps can therefore run as separate jobs, one per well, followed by a single `--all-wells --steps combine normalize outliers guides`. When steps run together, a step whose input was not written earlier in the same run stops rather than read an older table; naming it alone in `--steps` uses the existing table. + +`--components` selects the quality-control components: `plots`, `segmentation`, `pheno`, `pheno_sbs`, and `filters`, which writes the lists. A component whose prerequisites are missing, such as image directories or scallops, is skipped with a warning, or stops the run when named in `--components`. + +Plots are written only with `--save-plots`; `--validate` adds before-and-after distributions of the normalisation. `--continue-on-error` carries on with the remaining wells after a failure. `normalize` and `guides` use every core available to the run; the top-level `n_jobs`, or `--n-jobs`, sets a number instead. The log is written to `preprocessing/bywell/qc_reports/` under `local_output_dir`. + +## Outputs + +Tables are written under `data_output_dir`; filter lists, records and plots under `local_output_dir`. + +| File | Under | Contents | +|---|---|---| +| `preprocessing/bywell/filters/_*.parquet` | local | The well's filter lists, and `_filters.json`, the components run and the lists written | +| `preprocessing/bywell/filtered/-objects.parquet`, `-features.parquet`, `.parquet` | data | The well's filtered objects, features, and both joined | +| `preprocessing/allwells/combined/allwells-objects.parquet`, `allwells-features.parquet` | data | All wells stacked | +| `preprocessing/allwells/normalized/allwells-features_cell.parquet` | data | Normalised features | +| `preprocessing/allwells/outlier_filtered/allwells-features_cell_outlier_filtered.parquet` | data | Outlier cells removed; `_guide_outlier_filtered` holds the per-guide medians | +| `preprocessing/allwells/filters/outlier_filters.parquet` | local | Every cell removed and why; `outlier_filters.json` records the settings | +| `guide_filtering/guide_filtered/allwells-features_cell.parquet` | data | The final cell table, read by `chart archetypes` | +| `guide_filtering/guide_filtered/allwells-features_guide.parquet`, `allwells-features_gene.parquet` | data | Per-guide medians, and per-perturbation medians of those | +| `guide_filtering/filters/guide_filters.parquet` | local | Every guide removed and why; `guide_filters.json` records the settings | + +Cells are indexed by `Label`, `Guide`, `Gene` and `Well` from `combine` on. + +## Configuration + +```yaml +preprocessing_bywell: + thresholds: + segmentation_error_threshold: 15 # cytoplasm-to-nucleus area ratio + pheno_correlation: 0.8 # DAPI correlation across phenotyping rounds + pheno_sbs_correlation: 0.8 # DAPI correlation, phenotyping against sequencing + +preprocessing_allwells: + normalization_type: quantile + contamination: auto # or the proportion of cells to flag + random_state: 42 + n_estimators: 100 + missing_value_threshold: 100000 # a count of cells, not a fraction + rcv_threshold: 0.01 + +guide_filtering: + cosine_similarity_threshold: 0.2 + min_cells_per_guide: 50 +``` + +The values shown are the defaults. The stage also reads `schema`, and the input settings described under [Your data](data.md). + +The `pheno` and `pheno_sbs` components read their images through scallops, from these keys under `preprocessing_bywell.global` or `preprocessing_bywell.wells.`: + +| Key | Contents | +|---|---| +| `pheno_dir` | Phenotyping images | +| `dapi_channel_names` | The DAPI channel of each phenotyping round; `pheno` needs at least two | +| `pheno_sbs_registered_dir` | Phenotyping images registered to the sequencing images | +| `pheno_sbs_dapi_channel` | The DAPI channel among those, by name or position; by default the last of `dapi_channel_names` | +| `sbs_dir` | Sequencing images; DAPI is read from the first channel of the first cycle | +| `sbs_feature_dir` | Holds `nuclei/-objects.parquet`, the nuclei in sequencing coordinates | + +Filters computed elsewhere are added under `preprocessing_bywell.precomputed_filters`, one entry per filter, keyed by its name. Each entry gives `path`, a directory; `file_pattern`, default `{well}-objects.parquet`; optionally `threshold_column` and `threshold`, keeping the rows at or above the threshold; and `role`, `inclusion` (default) or `exclusion`. The table's index has to hold the cell labels. + +## Method + +**Quality control.** `segmentation` excludes cells whose cytoplasmic area exceeds `segmentation_error_threshold` times their nuclear area. `pheno` correlates, within each nucleus's bounding box, the first phenotyping round's DAPI image with each later round's, and keeps the cells whose mean Pearson correlation reaches `pheno_correlation`; a poor correlation marks a cell misregistered between rounds. `pheno_sbs` does the same between the phenotyping DAPI image registered to the sequencing images and the sequencing DAPI image, which ties a cell's phenotype to its guide assignment. A cell is kept if it appears in every inclusion list and in no exclusion list. + +**Normalisation.** Each feature is quantile-normalised across wells: a well's values are ranked, ties taking their average rank, and mapped by linear interpolation onto the feature's distribution pooled over the screen. Missing values take no part and stay missing. + +**Outliers.** Three prefilters run first: features with more than `missing_value_threshold` missing values are dropped, then features whose robust coefficient of variation (normal-scaled median absolute deviation over the median) is below `rcv_threshold`, then cells with any missing value. An isolation forest (`n_estimators` trees, seeded by `random_state`) is then fitted to the cells of each perturbation separately, and the cells it flags are removed; perturbations with fewer than two cells are not tested. `contamination: auto` uses the forest's own score threshold. + +**Guides.** The step works on per-guide medians of the outlier-filtered cells. Control guides are flagged when at least two of three detectors agree: an isolation forest on the features, and a local outlier factor (20 neighbours) and an elliptic envelope in a PCA space. The space is fitted to all guides after standardisation and keeps the components up to 90% of the variance. A second PCA space is fitted to the remaining guides and centred on the control guides, and each perturbation guide is scored by its mean cosine similarity to the other guides of its perturbation; guides below `cosine_similarity_threshold` are removed. Perturbations with a single guide cannot be scored and are kept. Last, guides with fewer than `min_cells_per_guide` cells left are removed. A guide's profile is the median of its cells, and a perturbation's the median of its guides', so that a guide with many cells does not dominate. + +## Troubleshooting + +| Message | Cause and remedy | +|---|---| +| `well ... has no QC filter lists in ...; run the qc step first` | QC wrote no list for the well; without images only `segmentation` produces one, and only if it flags a cell. Without `--well` the well is skipped with a warning instead | +| `Every one of the ... objects was filtered out ...` | The filter lists and the tables identify cells differently, often by label type, or the lists belong to another well | +| `No feature columns: none of the ... columns match any of the patterns ...` | `schema.feature_patterns` matches no column | +| `... of the ... columns to normalise hold something other than numbers ...` | `schema.feature_patterns` selects a text column; narrow the patterns | +| `Not normalising: ...`, `Not filtering outliers: ...`, `Not filtering guides: ...` | The previous step wrote nothing in this run; fix it, or name the step alone to use an existing table | +| `Prefiltering left ... cells and ... columns ...` | The prefilters removed everything; check `missing_value_threshold` and `rcv_threshold`, or `--skip-prefilters` | +| `No control guides among ... rows ...` | No perturbation matches `schema.control_gene` or `schema.control_prefix` | +| `Guide filtering removed every cell ...` | `cosine_similarity_threshold` or `min_cells_per_guide` is too strict for the screen | +| `'...' was requested for well ... but ...` | A component named in `--components` lacks its images, settings or scallops | diff --git a/docs/programs.md b/docs/programs.md new file mode 100644 index 0000000..44e0354 --- /dev/null +++ b/docs/programs.md @@ -0,0 +1,83 @@ +# Programs + +`chart programs` decomposes the perturbation × archetype effect matrix into morphological programs, sparse sets of archetypes across channels that perturbations shift together. Perturbations are grouped into modules by consensus clustering of their program scores, and every cell is scored on the same programs, which links each program back to individual cells and their images. + +The stage reads `effects/pooled.parquet` from `chart effects`, the per-channel archetype weights from `chart archetypes`, and the objects table from `chart preprocessing`. + +## Usage + +``` +chart programs --config my_screen.yaml +``` + +The four steps run in order, each reading the previous step's output; `--steps` selects a subset: + +``` +chart programs --config my_screen.yaml --steps modules cells +``` + +| Step | Output | +|---|---| +| `matrix` | Signed effects of the significant perturbations on the significant archetypes | +| `fit` | Program loadings and perturbation scores | +| `modules` | Consensus matrix and module assignments | +| `cells` | Cell scores and example cells per program | + +> [!NOTE] +> `modules` refits the decomposition 150 times; for 471 perturbations and 188 archetypes it takes about 20 minutes on 96 cores. + +The sparse PCA uses every core available to the run, which under SLURM is the job's allocation; the top-level `n_jobs`, or `--n-jobs`, sets a number instead. + +## Outputs + +All tables are written to `programs/` under `data_output_dir`. + +| File | Rows | Contents | +|---|---|---| +| `matrix.parquet` | perturbation | signed effect on each retained archetype | +| `loadings.parquet` | archetype | loading on each program | +| `scores.parquet` | perturbation | score on each program | +| `consensus.parquet` | perturbation | co-clustering frequency with every other perturbation | +| `modules.parquet` | perturbation | module, nearest other module, margin, core or weak | +| `cell_scores.parquet` | cell | score on each program | +| `top_cells.parquet` | program × rank | example cells with score, purity and bounding box | + +Positive loadings mark archetypes in which perturbations scoring high on the program gain weight relative to controls; negative loadings, archetypes in which they lose it. Perturbation and cell scores are both uncentred projections in centred log-ratio (CLR) units. A perturbation scores zero when it does not differ from controls; a cell scores zero when its weights are uniform across every channel's archetypes. + +A perturbation's `margin` is its mean consensus with its own module minus its highest mean consensus with any other (`best_other_module`). Perturbations with a margin below 0.15 are `weak`, the rest `core`. + +`matrix`, `scores`, `modules` and `cell_scores` each come with a `.json` record of the inputs and settings that produced them. + +## Configuration + +```yaml +programs: + outputs_dir: programs/ + k: 18 # number of programs + qval_threshold: 0.01 # significance for retaining perturbations and archetypes + n_modules: 25 # number of modules + seed: 42 +``` + +The defaults are those of the published analysis, 471 perturbations and 188 archetypes; `k` and `n_modules` should be chosen for the screen at hand. + +## Method + +**Effect matrix.** `chart effects` measures each perturbation's effect on each archetype as the Wasserstein-1 distance between the distributions of that archetype's CLR-transformed weight in its cells and in the control cells of the same wells, signed by the difference in medians (`signed_w1`). P-values come from permuting labels within wells, and Benjamini–Hochberg q-values are computed jointly across channels. The matrix keeps every perturbation with q < `qval_threshold` on at least one archetype, and every archetype with at least one such perturbation. + +**Decomposition.** Sparse PCA (scikit-learn, `k` components, at most 500 iterations, tolerance $10^{-4}$) factorises the matrix; in the published analysis about a quarter of the loadings are non-zero. Each program's sign is fixed so that its largest-magnitude loading is positive. Scores are the ridge projection $S = XL\,(L^\top L + 0.01\,I)^{-1}$, as in scikit-learn's `transform` but without centring $X$, so that zero means no effect. + +**Consensus modules.** The decomposition is refitted on 150 random subsamples of 85% of the perturbations (10 seeds × 15 subsamples). Each refit's programs are matched to the reference by Hungarian assignment on the absolute cosine similarity of their loadings and sign-aligned. Each refit's scores are then partitioned 38 ways: Ward and average-linkage hierarchical clustering of the z-scored scores, cut at 12, 14, …, 28 clusters, and Leiden clustering (RB configuration, resolution 0.8, 1.0, 1.2, 1.5, 1.8) of symmetrised nearest-neighbour graphs with 5, 8, 10 and 12 neighbours. The consensus of two perturbations is the fraction of partitions containing both that place them together. Modules are obtained by average-linkage clustering of 1 − consensus, cut into `n_modules` clusters. + +**Cell scores.** Each cell's weights are CLR-transformed per channel, as in `chart effects`, concatenated across channels and projected with the same ridge projection. A cell absent from a channel contributes zero there, as if its weights were uniform. For each program, the example cells are the 20 with the highest purity among its 100 top-scoring cells, purity being the cell's percentile on the program minus its mean percentile on the other programs. + +## Troubleshooting + +| Message | Cause and remedy | +|---|---| +| `No pooled.parquet in ...; run 'chart effects' first` | Pooled effects are written only once every channel in `archetypes.channels` has been tested; run `chart effects` on the rest | +| `No gene has a q-value below ...` | No effect passes `qval_threshold` | +| `... genes were not tested on every archetype kept ...` | Those perturbations' wells contain no control cells in some channel, so they have no effect there | +| `The loadings in ... were fitted on another matrix; rerun the fit step` | The matrix changed after the decomposition; rerun from `fit` | +| `The programs load on ..., which 'archetypes.archetypes_to_remove' now leaves out` | An archetype was excluded after the programs were fitted; rerun `chart effects`, then this stage from `matrix` | +| `The programs load on ..., which the weights do not hold` | The archetypes were refitted, or a channel removed from `archetypes.channels`, since the programs were fitted; rerun from `matrix` | diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100644 index 0000000..a2be2e2 --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1,3 @@ +sphinx +myst-parser>=5.1 +pydata-sphinx-theme diff --git a/pyproject.toml b/pyproject.toml index 3fd7552..2ca103e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,8 +34,8 @@ dependencies = [ # range of numpy versions than anything else here, so it is the dependency # most likely to constrain an install; see OPEN_ISSUES.md for the way out. -# scallops supplies the readers for registered images. It is not published -# on PyPI, so it has to be installed separately; without it CHART skips the +# scallops supplies the readers for registered images. It needs Python 3.12 +# or later, so it has to be installed separately; without it CHART skips the # components that need image data rather than failing. # The annotation stage is optional, so its packages are an extra.