diff --git a/.planning/codebase/ARCHITECTURE.md b/.planning/codebase/ARCHITECTURE.md new file mode 100644 index 00000000..af9dc806 --- /dev/null +++ b/.planning/codebase/ARCHITECTURE.md @@ -0,0 +1,227 @@ +--- +last_mapped_commit: 3c83c5b8aa313077e0ce43239a7b818b281b32ea +last_mapped_at: 2026-10-09 +--- + + +# Architecture + +**Analysis Date:** 2026-10-09 + +## System Overview + +deepCSA (`bbglab/deepCSA`) is a **Nextflow DSL2 pipeline** for analysis of the clonal structure of tissues using duplex-sequencing data. It is a workflow-orchestration architecture: Groovy/Nextflow processes orchestrate containerized Python/R scripts. + +```text +┌─────────────────────────────────────────────────────────────────────┐ +│ ENTRY POINT / ORCHESTRATION │ +│ `main.nf` → workflow BBGTOOLS → DEEPCSA │ +│ `workflows/deepcsa.nf` (798 lines — the pipeline "spine") │ +├──────────────────┬──────────────────┬───────────────────────────────┤ +│ SUBWORKFLOWS │ CONFIG LAYER │ PARAM VALIDATION │ +│ `subworkflows/ │ `nextflow.config`│ nf-schema plugin │ +│ local/*` (23) │ `conf/*.config` │ `nextflow_schema.json` │ +└────────┬─────────┴──────────────────┴───────────────────────────────┘ + │ + ▼ +┌─────────────────────────────────────────────────────────────────────┐ +│ PROCESS MODULES │ +│ `modules/local/**/main.nf` (97 process definitions, 48 top dirs) │ +│ `modules/nf-core/**` (multiqc, custom/dumpsoftwareversions) │ +│ `subworkflows/nf-core/**` (utils, vep annotation) │ +└────────┬────────────────────────────────────────────────────────────┘ + │ invokes via PATH (bin/ auto-mounted by Nextflow) + ▼ +┌─────────────────────────────────────────────────────────────────────┐ +│ ANALYSIS SCRIPTS (bin/) │ +│ ~90 Python/R scripts + shared utils (`bin/utils*.py`) │ +│ Run inside containers: docker.io/bbglab/*, quay.io registry │ +└────────┬────────────────────────────────────────────────────────────┘ + │ + ▼ +┌─────────────────────────────────────────────────────────────────────┐ +│ OUTPUT: ${params.outdir}/ (mutdensity/, mutational_profile/, │ +│ plots/, depths/, pipeline_info/ — publish paths in conf/*.config) │ +└─────────────────────────────────────────────────────────────────────┘ +``` + +## Component Responsibilities + +| Component | Responsibility | File | +|-----------|----------------|------| +| `main.nf` | Entry point; declares `BBGTOOLS` wrapper workflow and runs `PIPELINE_INITIALISATION` | `main.nf` | +| `DEEPCSA` workflow | The entire pipeline logic: channel wiring, feature flags, fan-out of analyses | `workflows/deepcsa.nf` | +| `PIPELINE_INITIALISATION` | Banner, version validation, params summary log | `subworkflows/local/utils_nfcore_deepcsa/main.nf` | +| Local subworkflows | Compose modules into logical analysis stages (depths, panels, omega, signatures...) | `subworkflows/local/*/main.nf` | +| Local modules | Single `process` definitions with container/conda, script invocation, versions.yml | `modules/local/**/main.nf` | +| nf-core modules/subworkflows | Vendored community components (MultiQC, VEP annotation, utils) | `modules/nf-core/`, `subworkflows/nf-core/` | +| `bin/` scripts | Actual analysis logic (Python with click/pandas, some R) | `bin/*.py`, `bin/*.R` | +| Config layer | Params defaults, resources, publish paths, tool/mode presets | `nextflow.config`, `conf/*.config` | + +## Pattern Overview + +**Overall:** nf-core-style DSL2 pipeline with local (non-nf-core) module tree. + +**Key Characteristics:** +- **Single mega-workflow:** All pipeline logic lives in `workflow DEEPCSA` inside `workflows/deepcsa.nf`. `main.nf` is a thin wrapper (`BBGTOOLS` just calls `DEEPCSA`). +- **Feature-flag fan-out:** Analyses are toggled by boolean params (`params.omega`, `params.oncodrivefml`, `params.dnds`, `params.signatures`, `params.mutationdensity`, `params.profileall`, ...) checked with `if` blocks in the workflow body. +- **Subworkflow reuse via aliasing:** The same subworkflow is included multiple times with aliases, e.g. `MUTATION_DENSITY as MUTDENSITYALL/PROT/NONPROT/SYNONYMOUS` (`workflows/deepcsa.nf:33-37`). +- **Meta-map tuples:** Data flows as `tuple val(meta), path(file)` where `meta.id` is the sample/group key. +- **Channel aggregation with `collectFile`:** Per-sample outputs are flattened and concatenated into cohort-level files (`all_mutdensities.tsv`, `all_profile_stabilities.tsv`) directly in the workflow. +- **Result accumulation channel:** `positive_selection_results` is progressively `.join(..., remainder: true)`-ed with each positive-selection tool's output, then filtered and fed to `PLOTTINGSUMMARY`. + +## Layers + +**Orchestration layer:** +- Purpose: Entry point and pipeline initialization +- Location: `main.nf` +- Contains: `BBGTOOLS` wrapper workflow, `PIPELINE_INITIALISATION` invocation +- Depends on: `workflows/deepcsa.nf`, `subworkflows/local/utils_nfcore_deepcsa` +- Used by: `nextflow run` CLI + +**Workflow layer (the spine):** +- Purpose: All channel wiring, branching, and analysis fan-out +- Location: `workflows/deepcsa.nf` +- Contains: `workflow DEEPCSA` (~700 lines of channel logic), imports of all subworkflows/modules +- Depends on: every subworkflow in `subworkflows/local/`, many modules in `modules/local/` and `modules/nf-core/` +- Used by: `main.nf` + +**Subworkflow layer:** +- Purpose: Reusable multi-process analysis stages +- Location: `subworkflows/local//main.nf` (23 workflows) +- Contains: `include` of modules + `workflow X { take: main: emit: }` blocks +- Key stages: `depthanalysis`, `createpanels`, `mutationpreprocessing`, `mutationdensity`, `mutationprofile`, `mutability`, `omega`, `oncodrivefml`, `oncodrive3d`, `oncodriveclustl`, `dnds`, `indels`, `signatures`, `signatures_hdp`, `mutatedcells`, `regressions`, `plotdepths`, `plotting_qc`, `plottingsummary`, `enrichpanels`, `adjmutdensity`, `input_check.nf` +- Depends on: `modules/local/**`, `modules/nf-core/**` +- Used by: `workflows/deepcsa.nf` + +**Module (process) layer:** +- Purpose: One Nextflow `process` per file, wrapping a script/tool invocation +- Location: `modules/local///main.nf` (97 `main.nf` files) +- Contains: process definition with `conda`, `container`, `input/output` (meta tuples, named `emit:`), `script:` heredoc, `stub:` block, `versions.yml` emission +- Depends on: `bin/` scripts (resolved via Nextflow's automatic `bin/` PATH) and container images +- Used by: subworkflows and directly by `workflows/deepcsa.nf` + +**Script layer:** +- Purpose: Actual computation logic +- Location: `bin/*.py`, `bin/*.R`, `bin/saturation_mutagenesis/` +- Contains: Python 3 CLI scripts (click + pandas + pybedtools/polars), R scripts (dNdScv, mutrate), shared helpers (`bin/utils.py`, `bin/utils_context.py`, `bin/utils_filter.py`, `bin/utils_impacts.py`, `bin/utils_plot.py`, `bin/read_utils.py`) +- Depends on: conda envs declared per-process (`python=3.10`, `pybedtools`, `polars`, `click`) +- Used by: module processes + +**Config layer:** +- Purpose: Defaults, resources, publishing, tool presets +- Location: `nextflow.config` (params defaults + profiles), `conf/base.config` (resource limits, labels, retry strategy), `conf/modules.config` (78 `withName` overrides: `ext.args`, `publishDir`), `conf/results_outputs.config` (final publish paths), `conf/tools/*.config` (per-tool: omega, oncodrive3d, oncodrivefml, mutdensity, panels, regressions, hdp), `conf/modes/*.config` (presets: `basic`, `clonal_structure`, `get_signatures`) +- Depends on: `nextflow_schema.json` (nf-schema validation), `assets/schema_input.json` (samplesheet schema) +- Used by: Nextflow launcher + +## Data Flow + +### Primary Request Path + +1. **Input validation** — samplesheet CSV validated against `assets/schema_input.json` via nf-schema; `INPUT_CHECK` → `SAMPLESHEET_CHECK` builds `sample_inputs_ch` of `[meta, vcf, bam]` (`subworkflows/local/input_check.nf`). Alternative entry: `--input_maf` + `--use_custom_depths` converts MAF→VCF via `INPUTMAF2VCF` (`workflows/deepcsa.nf:224-238`). +2. **Grouping definition** — `TABLE2GROUP` parses the features table into JSON group definitions (`modules/local/table2groups/main.nf`); group keys are extracted with `JsonSlurper` in the workflow (`workflows/deepcsa.nf:247-260`). +3. **Depth computation** — `DEPTHANALYSIS` computes per-position depths from BAMs (`COMPUTEDEPTHS`) or accepts a custom depths table, then filters by minimum depth (`subworkflows/local/depthanalysis/main.nf`). +4. **Panel creation** — `CREATEPANELS` builds consensus panels (all/exons/protein-coding/non-coding/synonymous BEDs + VEP-annotated panels) from depths + WGS trinucleotide counts (`subworkflows/local/createpanels/main.nf`, `modules/local/createpanels/{captured,consensus,compare,custombedfile}`). +5. **Mutation preprocessing** — `MUT_PREPROCESSING` annotates VCFs with Ensembl VEP (`subworkflows/nf-core/vcf_annotate_ensemblvep*`), blacklists artifacts, builds a mask matrix; emits `somatic_mafs` (`subworkflows/local/mutationpreprocessing/main.nf`). +6. **Depth annotation & enrichment** — `ANNOTATEDEPTHS` merges depths with panels; optional `DOWNSAMPLEDEPTHS`; `ENRICHPANELS` expands panels with subgenic regions and DNA→protein/domain mappings (`subworkflows/local/enrichpanels/main.nf`). +7. **Parallel analyses (feature-flag gated)** — each consumes `somatic_mutations` + relevant consensus BED/panel + depths: + - Mutational profiles: `MUTPROFILEALL/NONPROT/EXONS/INTRONS` + - Mutation density: `MUTDENSITYALL/PROT/NONPROT/SYNONYMOUS` + adjusted variant `MUTDENSITYADJUSTED` → `DNDSPROXY` + - Mutability: `MUTABILITYALL/NONPROT` (feeds oncodrivefml/3d/clustl) + - Positive selection: `ONCODRIVEFMLALL`, `ONCODRIVE3D`, `ONCODRIVECLUSTL`, `DNDS`, `OMEGA` (+ `OMEGAMULTI`, `OMEGANONPROT` variants), `INDELSSELECTION` + - Clonal structure: `EXPECTEDMUTATEDCELLS`, `MUTATEDCELLSVAF`, `VAFSMOOTHING` + - Signatures: `MAF2VCF` → `SIGPROMATRIXGENERATOR` → `SIGNATURESALL/...` → `MUTS2SIGS` +8. **QC & summary** — `PLOTTINGQC` flags failing omegas/QC metrics; results are combined into `positive_selection_results_ready` and passed to `PLOTTINGSUMMARY` (`workflows/deepcsa.nf:640-712`). +9. **Regressions (optional)** — `REGRESSIONSMUTDENSITY/OMEGA/OMEGAGLOB` correlate metrics with metadata (`subworkflows/local/regressions/main.nf`). +10. **Publishing** — `publishDir` paths defined per-process in `conf/modules.config` and `conf/results_outputs.config`; `CUSTOM_DUMPSOFTWAREVERSIONS` collates `versions.yml` topics; MultiQC report at the end. + +### Aggregation Pattern (secondary flow) + +1. Per-sample process emits `tuple val(meta), path(result.tsv)`. +2. Workflow maps to file only: `.out.mutdensities.map{ it -> it[1] }.flatten()`. +3. `.collectFile(name: "all_*.tsv", storeDir: "${params.outdir}/...", skip: 1, keepHeader: true)` writes the cohort file directly to the outdir (`workflows/deepcsa.nf:316-322, 395-410`). + +**State Management:** +- No persistent state; everything is Nextflow channels + task work dirs. +- Placeholder channels (`channel.value(file("${projectDir}/assets/placeholder_no_file.tsv"))`) initialize optional result channels so downstream `.join(..., remainder: true)` never blocks (`workflows/deepcsa.nf:141-148`). +- `scratchhhh/` is a gitignored developer scratch area (not part of the pipeline). + +## Key Abstractions + +**Meta map (`meta`):** +- Purpose: Sample/group identity carried with every file tuple +- Examples: `workflows/deepcsa.nf` (`meta.id`, `meta.batch`), `subworkflows/local/input_check.nf` (`create_input_channel`) +- Pattern: `tuple val(meta), path(file)`; `meta.id` used for file prefixes via `task.ext.prefix` + +**Process module:** +- Purpose: One tool invocation, fully self-contained +- Examples: `modules/local/createpanels/consensus/main.nf`, `modules/local/bbgtools/omega/estimator/main.nf` +- Pattern: `tag "$meta.id"`, `conda "..."`, `container 'docker://bbglab/...'`, `label 'cpu_medium'`, named `emit:` outputs, `versions.yml` with `topic: versions`, `stub:` block for testing + +**Feature-flag params:** +- Purpose: Enable/disable analysis branches +- Examples: `nextflow.config` params block (`omega`, `omega_multi`, `omega_globalloc`, `oncodrivefml`, `oncodrive3d`, `dnds`, `signatures`, `mutationdensity`, `profileall`, `regressions`, `downsample`, ...) +- Pattern: derived booleans in the workflow (`run_mutabilities`, `run_mutdensity`, `run_profile_all` at `workflows/deepcsa.nf:206-208`) + +**Grouping JSONs:** +- Purpose: Define sample/gene groupings used by omega, regressions, plotting +- Examples: `modules/local/table2groups/main.nf`, `assets/omega_consequences_groupings.json` +- Pattern: `TABLE2GROUP` emits `json_samples`/`json_groups`/`json_allgroups`; keys parsed in workflow to build `samples_keys_ch`/`group_keys_ch` + +## Entry Points + +**`nextflow run main.nf` (or repo root):** +- Location: `main.nf` +- Triggers: CLI with `--input` samplesheet (or `--input_maf`), profile (`docker`/`singularity`/`conda`), optional `-c conf/modes/.config` +- Responsibilities: init (banner, version check, params validation via `PIPELINE_INITIALISATION`), then run `DEEPCSA` + +**Tests:** +- Location: `tests/deepcsa.nf.test` (nf-test with snapshot `tests/deepcsa.nf.test.snap`), `bin/test/test_*.py` (pytest) +- Triggers: `nf-test test` (config in `nf-test.config`, profile `test,singularity`) + +## Architectural Constraints + +- **Threading:** Nextflow manages parallelism; per-process resources via labels `cpu_low` (2 cpus), `cpu_medium` (4), `cpu_high` (8), `mem_low` (1 GB) in `conf/base.config`; global caps `params.max_cpus/max_memory/max_time`; per-process `withName` overrides for heavy steps (e.g. `CREATEPANELS:SITESFROMPOSITIONS` 8 GB). +- **Error handling:** `errorStrategy` retries with exponential backoff for exit codes 130–145 and 104 (`conf/base.config`); labels `error_ignore` and `error_retry` opt processes out/in. +- **Containers:** Default registry `quay.io` (`docker.registry`/`singularity.registry` in `nextflow.config`); bbglab tools pinned images (e.g. `docker.io/bbglab/omega:0.2.1`, `bbglab/deepcsa_bed:latest`). Processes declare both `conda` and `container`. +- **bin/ on PATH:** Nextflow auto-prepends `bin/` to the task PATH — scripts in `bin/` are callable by name from any process script block. Do not duplicate scripts elsewhere. +- **Global state:** None beyond `params`; workflow-level Groovy variables (`positive_selection_results`, `all_compiled_omegas`, ...) are channel builders inside `DEEPCSA`. +- **Vendored nf-core code:** Only 2 nf-core modules (`multiqc`, `custom/dumpsoftwareversions`) and 5 nf-core subworkflows are vendored and pinned in `modules.json`; nf-core dirs are excluded from nf-test (`ignore 'modules/nf-core/**/*'` in `nf-test.config`). +- **Circular imports:** None observed; dependency direction is strictly `main.nf → workflows → subworkflows → modules → bin/`. + +## Anti-Patterns + +### Monolithic workflow file + +**What happens:** `workflows/deepcsa.nf` is ~800 lines containing all channel wiring, JSON parsing, aggregation, and branching. +**Why it's wrong:** Hard to test in isolation; any change risks breaking unrelated branches; merge conflicts likely. +**Do this instead:** Follow the existing subworkflow pattern — group related processes into `subworkflows/local//main.nf` with `take/main/emit` (see `subworkflows/local/depthanalysis/main.nf` as the model) and keep `workflows/deepcsa.nf` as wiring only. + +### Inline Groovy logic in the workflow + +**What happens:** `JsonSlurper` parsing of group JSONs and result-channel reshaping happen inline in `workflows/deepcsa.nf:247-260, 713-720`. +**Why it's wrong:** Untestable outside Nextflow; mixes orchestration with data transformation. +**Do this instead:** Push parsing into a `bin/` script invoked by a module process, or a small subworkflow. + +## Error Handling + +**Strategy:** Process-level retry with exponential backoff; validation fails fast. + +**Patterns:** +- `errorStrategy` retry on transient exit codes (130–145, 104) with `sleep(Math.pow(2, task.attempt) * 200)` (`conf/base.config`) +- Labels `error_ignore` / `error_retry` for per-process overrides (`conf/base.config`) +- Explicit `error "..."` guards in the workflow for invalid param combinations (e.g. `--input_maf` without `--use_custom_depths`, `--omega_covariates` without `--omega`) (`workflows/deepcsa.nf:210-218`) +- nf-schema validation of params and samplesheet (`nextflow_schema.json`, `assets/schema_input.json`) +- Python scripts raise/exit non-zero on bad input; shared helpers in `bin/utils.py`, `bin/utils_filter.py` + +## Cross-Cutting Concerns + +**Logging:** Nextflow task logs; `PIPELINE_INITIALISATION` prints banner + params summary (`paramsSummaryMap` from nf-schema plugin); MultiQC aggregates QC (`assets/multiqc_config.yml`). +**Validation:** nf-schema for params (`params.validate_params`), JSON schema for samplesheet (`assets/schema_input.json`), `SAMPLESHEET_CHECK` module. +**Authentication:** None (HPC/cluster execution; no external API auth). Seqera Platform integration via `tower.yml` (report display only). +**Versioning:** Every process emits `versions.yml` (topic `versions`), collated by `CUSTOM_DUMPSOFTWAREVERSIONS`; pipeline version in `params.version` checked at init. +**Publishing:** Two-tier publish config — default per-process paths in `conf/modules.config`, curated final outputs in `conf/results_outputs.config`; `versions.yml` never published (`saveAs` filter). + +--- + +*Architecture analysis: 2026-10-09* diff --git a/.planning/codebase/CONCERNS.md b/.planning/codebase/CONCERNS.md new file mode 100644 index 00000000..6cdb0ffd --- /dev/null +++ b/.planning/codebase/CONCERNS.md @@ -0,0 +1,210 @@ +--- +last_mapped_commit: 3c83c5b8aa313077e0ce43239a7b818b281b32ea +last_mapped_at: 2026-10-09 +--- +# Codebase Concerns + +**Analysis Date:** 2026-10-09 + +## Tech Debt + +**Unpinned container images (`latest` / untagged):** +- Issue: Many local modules reference floating or untagged Docker images, so rebuilds can silently pick up different tool versions and break reproducibility. +- Files: `modules/local/*/main.nf`, `subworkflows/local/*/main.nf`, `conf/modules.config` + - `docker://bbglab/deepcsa_bed:latest` (`modules/local/createpanels/consensus/main.nf`) + - `docker.io/ferriolcalvet/dnds:latest`, `docker.io/ferriolcalvet/oncodrivefml:latest`, `docker.io/ferriolcalvet/oncodriveclustl:latest`, `docker.io/ferriolcalvet/musical:latest`, `docker.io/ferriolcalvet/msighdp:latest`, `docker.io/axelrosendahlhuber/expected_mutrate:latest` + - Untagged: `docker.io/ferriolcalvet/hdp_wrapper`, `docker.io/ferranmuinos/test_mutated_genomes`, `docker.io/ferriolcalvet/sigprofilerassignment` + - Dev-tagged: `docker.io/rblancomi/bbgregressions:dev` (`conf/modules.config`) +- Impact: Non-reproducible results across runs/nodes; a pushed `latest` image can break the pipeline without any repo change. +- Fix approach: Pin every container to an immutable version tag (pattern already used by `docker.io/bbglab/omega:0.2.1`, `docker.io/ferriolcalvet/sigprofilermatrixgenerator:1.3.5`, `docker.io/spellegrini87/oncodrive3d:1.0.9-light`). + +**No Python dependency manifest for `bin/` scripts:** +- Issue: `pyproject.toml` only configures Black/isort; there is no `requirements.txt` or dependency list for the ~86 Python scripts in `bin/` (pandas, numpy, polars, click, matplotlib, etc.). +- Files: `pyproject.toml`, `bin/*.py` +- Impact: Local development and unit tests (`bin/test/`) depend on whatever is in the current conda env; version drift between the `deepcsa-core` container and local envs causes subtle failures. +- Fix approach: Add a pinned `requirements.txt` (or `[project.dependencies]` in `pyproject.toml`) matching the `bbglab/deepcsa-core` container contents. + +**Pending pandas 2.2.3 upgrade:** +- Issue: Two scripts carry `# TODO: bump pandas to 2.2.3`, indicating known incompatibilities with newer pandas. +- Files: `bin/concat_sbs_probs.py:3`, `bin/mut_density_simple.py:16` +- Impact: Blocks dependency modernization; mixed pandas/polars codebase (`bin/create_consensus_panel.py`, `bin/create_panel_versions.py`, `bin/merge_annotation_depths.py`, `bin/panel_custom_processing.py`, `bin/panels_computedna2protein.py` use polars) increases maintenance surface. +- Fix approach: Audit pandas-2.x breaking changes in these scripts, upgrade, and remove the TODOs. + +**Duplicated annotation post-processing logic:** +- Issue: `bin/postprocessing_annotation.py` (263 lines) and `bin/panel_postprocessing_annotation.py` (215 lines) are near-duplicates with the same TODOs repeated in both (`bin/postprocessing_annotation.py:85,190`, `bin/panel_postprocessing_annotation.py:70,164`). +- Impact: Bug fixes must be applied twice; the two files can drift (one already imports `utils.vartype`, the other does not). +- Fix approach: Extract shared consequence-mapping/muttype-conversion logic into `bin/utils_impacts.py` or a new shared module and parameterize the differences. + +**Stub CHANGELOG and incomplete nf-core template remnants:** +- Issue: `CHANGELOG.md` still contains the template placeholder `## v1.0dev - [date]` with empty Added/Fixed/Dependencies/Deprecated sections, while `nextflow.config` manifest declares `version = '1.0.1.dev'`. `subworkflows/local/utils_nfcore_deepcsa/main.nf:199-201` still has the placeholder Zenodo DOI TODO. +- Files: `CHANGELOG.md`, `nextflow.config:392`, `subworkflows/local/utils_nfcore_deepcsa/main.nf` +- Impact: Release history is not tracked; citation text may render incomplete. +- Fix approach: Backfill changelog per release; register Zenodo DOI and fill `toolCitationText`/`toolBibliographyText`. + +**Scratch/experiment artifacts in repo root:** +- Issue: A 1.1 GB `scratchhhh/` directory (SigProfiler outputs, VCFs) sits in the repo root. It is gitignored (`.gitignore` line `scratchhhh/`) but its name suggests it was renamed to dodge cleanup; stray files `nf-2eI0igpGfSq5RI-reports.tsv`, `.nextflow.log*` (10 rotated logs, ~2 MB total) also accumulate at root. `bin/explore_saturation.ipynb` (an exploration notebook) is git-tracked inside `bin/`. +- Files: `scratchhhh/`, `nf-2eI0igpGfSq5RI-reports.tsv`, `.nextflow.log*`, `bin/explore_saturation.ipynb` +- Impact: Confuses repo structure; risk of accidentally committing large outputs; notebooks in `bin/` mix exploration code with pipeline scripts. +- Fix approach: Move scratch data outside the repo or into a dedicated `scratch/` (already gitignored); delete rotated logs; move exploration notebooks to `assets/useful_scripts/` (where `.ipynb` files are already the convention and gitignored). + +**Unexplained/undocumented functions:** +- Issue: Several functions carry explicit "explain what this does" TODOs. +- Files: `bin/mutations_custom_processing.py:17`, `bin/mut_density_adjusted_dnds.py:12-13` (also requests a log file output) +- Impact: Onboarding cost and higher risk of misuse. +- Fix approach: Write docstrings and add the requested statistics log output. + +## Known Bugs + +**Five failing nf-test process tests (as of 2026-10-09 run):** +- Symptoms: `tests/2026-10-09_results.csv` records 5 FAILED / 20 PASSED: + - `modules/local/filtermaf/tests/main.nf.test` — "Should run and emit cohort filtered mutations" (FILTER_BATCH) + - `modules/local/group_genes/tests/main.nf.test` — GROUP_GENES genes-to-groups JSON + - `modules/local/mut_density/simple/tests/main.nf.test` — MUTATION_DENSITY all-panel densities + - `modules/local/sig_matrix_concat/tests/main.nf.test` — MATRIX_CONCAT concatenated WGS matrices + - `modules/local/sitesfrompositions/tests/main.nf.test` — SITESFROMPOSITIONS panel sites chunk +- Files: `modules/local/filtermaf/`, `modules/local/group_genes/`, `modules/local/mut_density/`, `modules/local/sig_matrix_concat/`, `modules/local/sitesfrompositions/` +- Trigger: Run `nf-test` suite (see `nf-test.config`); results CSVs in `tests/` (gitignored via `tests/2*results*`). +- Workaround: None — these indicate current regressions or stale snapshots (`tests/deepcsa.nf.test.snap` may also need regeneration). +- Note: The results CSV references test files (e.g. `modules/local/blacklistmuts/tests/main.nf.test`) that no longer exist in the working tree — only 2 `*.nf.test` files are git-tracked under `modules/local`. The suite was run against a different tree state; reconcile test files with the results. + +**Bedtools coordinate hack:** +- Symptoms: `bin/createcustombed.sh:21` contains `#HACK to handle issue arq5x/bedtools2#359` — an awk workaround mutating BED start/end coordinates. +- Files: `bin/createcustombed.sh` +- Trigger: Custom BED file creation path (`use_custom_bedfile = true`). +- Workaround: In place; will need revisiting if bedtools version changes (container `bbglab/deepcsa_bed:latest` is unpinned, compounding risk). + +**Dummy channel value in OMEGA subworkflow:** +- Symptoms: `subworkflows/local/omega/main.nf:68` — `// FIXME here I am using bedfile as a dummy value channel` passed into PREPROCESSING. +- Files: `subworkflows/local/omega/main.nf` +- Impact: Fragile channel wiring; refactoring the OMEGA inputs can silently misalign channels. +- Fix approach: Replace with an explicit value channel or remove the unused input. + +## Security Considerations + +**No secrets handling in repo (good), but no CI either:** +- Risk: No GitHub workflows exist (`.nf-core.yml` template `skip: [github, ci]`), so no automated linting/tests/secret scanning run on pushes. +- Files: `.nf-core.yml`, absent `.github/workflows/` +- Current mitigation: `.pre-commit-config.yaml` exists locally (prettier/markdownlint only). +- Recommendations: Add at minimum a CI workflow running `nf-test` and the `bin/test/` unittest suite; add nf-core `linting.yml` equivalent. + +**External subprocess calls:** +- Risk: `bin/merge_annotation_depths.py:161,175,183,191` shells out to `gzip` via `subprocess.run(..., check=True)` (no `shell=True`, arguments are list-form — low risk). +- Files: `bin/merge_annotation_depths.py` +- Current mitigation: List-argument invocation, `check=True`. +- Recommendations: None urgent; keep avoiding `shell=True` with user-derived paths. + +**Hardcoded reference/cache versions:** +- Risk: `vep_cache_version = 111` and `dnds_biomart_ref = "homo_sapiens.v111.canonical.biomart.tsv"` are hardcoded defaults (`nextflow.config:138,159`); annotation results silently depend on cache availability and version. +- Files: `nextflow.config`, `conf/` +- Current mitigation: Parameters are overridable via CLI/config. +- Recommendations: Document the required VEP cache version in `docs/usage.md` and validate cache presence at pipeline init. + +## Performance Bottlenecks + +**Memory-hungry annotation/mapping processes:** +- Problem: `DNA2PROTEINMAPPING` requests `30.GB * task.attempt` and is set to `errorStrategy = 'ignore'` after retries; `PLOTMUTATIONSPECIFIC` and `ONCODRIVEFMLSNVS` request `24.GB * task.attempt`. +- Files: `conf/base.config:211-216,220-222`, `conf/tmp_quick_fixes.config:9-12` +- Cause: Whole-panel tables loaded into pandas in one shot (e.g. `bin/panels_computedna2protein.py` uses polars, but several annotation scripts still use pandas `read_table` on full files). +- Improvement path: Extend the existing chunking mechanism — `params.panel_sites_chunk_size` (default 1,000,000, `nextflow.config:119`, consumed at `conf/modules.config:723` for `SITESFROMPOSITIONS`) — to other site-level processes; prefer polars/lazy evaluation in the largest scripts. + +**Silent output loss on OOM:** +- Problem: Because `DNA2PROTEINMAPPING` and several plot processes end in `errorStrategy = 'ignore'`, an OOM failure produces missing downstream outputs rather than a pipeline failure. +- Files: `conf/base.config:201-222`, `conf/modules.config:227` (`PLOTOMEGA`), `conf/tmp_quick_fixes.config` +- Improvement path: Convert to `retry` with a final `fail` (or emit an explicit empty-output marker) so missing results are detectable. + +## Fragile Areas + +**Positional column renaming in omega QC:** +- Files: `bin/omega_syn_qc.py:187` — `# TODO: dangerous column renaming here` — assigns `syn_muts_df.columns = [...]` positionally after a merge with suffixes `_loc`/`_gloc`. +- Why fragile: Any column-order change in the OMEGA container output (`bbglab/omega:0.2.1`) silently mislabels observed/estimated synonymous counts. +- Safe modification: Rename by explicit column-name mapping and assert expected columns before assignment. +- Test coverage: None (no unit test for this script). + +**Column-name robustness in annotation post-processing:** +- Files: `bin/postprocessing_annotation.py:162` (`# TODO: Is it robust enough to use columns names here?`), `bin/postprocessing_annotation.py:146` (bare `# TODO`), `bin/utils_impacts.py:196` (try/except around consequence handling flagged for revision), `bin/utils_context.py:33` (`# TODO remove this try-except`). +- Why fragile: VEP output format changes break parsing; bare `except:` clauses (`bin/postprocessing_annotation.py:178`, `bin/utils_context.py:42`, `bin/utils_impacts.py:202`, `bin/mutgenomes_driver_priority.py:30,43,134,259`, `bin/plot_depths.py:687`) swallow all errors including `KeyboardInterrupt`. +- Safe modification: Replace bare `except:` with narrow exception types; validate expected columns at read time. +- Test coverage: Only `bin/test/` covers `utils_filter`, `mask_matrix`, `check_samplesheet`, `check_contamination`, `plot_selectionsideplots` — none of the fragile files above. + +**Monolithic plotting scripts:** +- Files: `bin/plot_gene_saturation.py` (1546 lines, with `# TODO:`/`# FIXME` at lines 30, 1234, 1250, 1333), `bin/mutgenomes_expected_mutrisk.R` (836 lines, TODO at 49, FIXME at 218), `bin/plot_selectionsideplots.py` (796 lines), `bin/plot_depths.py` (771 lines). +- Why fragile: Single-file scripts mixing parsing, computation, and plotting; hard to test in isolation. +- Safe modification: Extract pure-computation helpers (pattern already proven by `bin/utils_filter.py` + `bin/test/test_utils_filter.py`). +- Test coverage: None for these scripts. + +**Giant central workflow:** +- Files: `workflows/deepcsa.nf` (798 lines, 89 `params.` references, ~40 subworkflow includes with repeated aliasing of the same subworkflow, e.g. `MUTATION_DENSITY` included 5×, `MUTATIONAL_PROFILE` 4×, `OMEGA_ANALYSIS` 4×, `SIGNATURES` 4×). +- Why fragile: Any channel rename ripples across all aliased invocations; conditional logic (`if (run_profile_all)` etc.) is hard to test. +- Safe modification: Keep changes local to one aliased block; verify with `tests/deepcsa.nf.test` pipeline-level test after edits. + +**Config-level failure masking:** +- Files: `conf/tmp_quick_fixes.config` (entire file is patches: `DNA2PROTEINMAPPING` retry-then-ignore, `COMPARE_SIGNATURES` ignore, `EXPECTEDMUTATEDCELLS` ignore), `conf/base.config:201-222`, `conf/exome.config:136-147`. +- Why fragile: "Quick fixes" have become permanent infrastructure; failures are invisible. +- Safe modification: Before removing an ignore, confirm the upstream bug is fixed; add explicit warning logs when a process is skipped. + +## Scaling Limits + +**Whole-cohort in-memory operations:** +- Current capacity: Cohort-level MAF/depth tables processed as single DataFrames; `DNA2PROTEINMAPPING` already needs 30 GB. +- Limit: Cohorts with large panels/WGS-scale site counts will exceed typical node memory in annotation and plotting steps. +- Scaling path: Roll out `panel_sites_chunk_size` chunking (implemented for `SITESFROMPOSITIONS`, `conf/modules.config:723`) to other processes; shard per-chromosome where possible. + +**Test runtime:** +- Current capacity: Process-level nf-tests take ~20-80 s each (per `tests/2026-10-09_results.csv`); the pipeline-level test (`tests/deepcsa.nf.test`) runs full `main.nf`. +- Limit: Adding tests for all 48 `modules/local/` processes at current per-test cost makes CI slow. +- Scaling path: Use minimal test data (`tests/test_data/` is only 20 KB — good), parallelize nf-test, and consider tagging slow tests. + +## Dependencies at Risk + +**`latest`-tagged method containers (see Tech Debt above):** +- Risk: `dnds`, `oncodrivefml`, `oncodriveclustl`, `musical`, `msighdp`, `deepcsa_bed`, `expected_mutrate`, `bbgregressions:dev`, plus three untagged images. +- Impact: Silent behavioral change of core selection/signature methods. +- Migration plan: Pin versions; publish versioned images for `hdp_wrapper`, `test_mutated_genomes`, `sigprofilerassignment`. + +**nf-core modules nearly absent:** +- Risk: `modules.json` tracks only 2 nf-core modules (`custom/dumpsoftwareversions`, `multiqc`) at a single git SHA (`911696ea0b62df80e900ef244d7867d177971f73`); everything else is local code in `modules/local/` (48 processes) and `subworkflows/local/` (24 subworkflows). +- Impact: No upstream fixes/updates flow in; local modules lack `meta.yml` in many cases (e.g. `modules/local/blacklistmuts/` has only `main.nf`). +- Migration plan: Incrementally port stable local processes to nf-core module standards (container pinning + `meta.yml` + tests). + +**Mixed pandas/polars codebase:** +- Risk: pandas pinned mentally to <2.2.3 (TODOs in `bin/concat_sbs_probs.py`, `bin/mut_density_simple.py`) while newer scripts use polars. +- Impact: Upgrade friction; two APIs to maintain. +- Migration plan: Complete pandas 2.2.3 bump, then standardize new development on polars for large tables. + +## Missing Critical Features + +**No CI pipeline:** +- Problem: `.nf-core.yml` skips all GitHub workflows; nothing runs tests or linting on push/PR. +- Blocks: Reliable regression detection (5 tests are already failing unnoticed in local runs). + +**No per-module tests for the vast majority of local modules:** +- Problem: 48 processes in `modules/local/` but only 1 `*.nf.test` file exists in the tree (`modules/local/expand_regions/tests/main.nf.test`) plus 2 tracked elsewhere; `subworkflows/local/` has none (only nf-core vendored subworkflows have tests). +- Blocks: Safe refactoring of `workflows/deepcsa.nf` and local modules. + +**Pipeline-level test coverage of feature flags:** +- Problem: `tests/deepcsa.nf.test` covers basic run, omega, and MAF-input validation; the many boolean feature switches in `nextflow.config` (e.g. `oncodrive3d`, `dnds`, `indels`, `signatures`, `omega_covariates`, `downsample`, `regressions`, `contamination`) lack dedicated pipeline tests. +- Blocks: Confidence that optional branches still work after workflow edits. + +## Test Coverage Gaps + +**`bin/` Python scripts:** +- What's not tested: ~81 of 86 scripts have no unit test; only `bin/test/test_check_samplesheet.py`, `bin/test/test_check_contamination.py`, `bin/test/test_mask_matrix.py`, `bin/test/test_plot_selectionsideplots.py`, `bin/test/test_utils_filter.py` exist (unittest style, `sys.path.insert` sibling imports). +- Files: `bin/*.py`, `bin/test/` +- Risk: Silent numeric errors in selection statistics (omega/dNdS QC scripts) and plotting code. +- Priority: High for `bin/omega_syn_qc.py`, `bin/postprocessing_annotation.py`, `bin/utils_impacts.py`, `bin/mutgenomes_driver_priority.py` (all contain bare excepts or flagged-fragile logic). + +**Local Nextflow modules:** +- What's not tested: 47 of 48 `modules/local/` processes; all 24 `subworkflows/local/` subworkflows. +- Files: `modules/local/`, `subworkflows/local/` +- Risk: Regressions like the 5 currently failing tests recur undetected. +- Priority: High for `filtermaf`, `group_genes`, `mut_density`, `sig_matrix_concat`, `sitesfrompositions` (currently failing); Medium for processes on the default execution path (`createpanels`, `mutationpreprocessing`). + +**R scripts:** +- What's not tested: `bin/dNdS_run.R`, `bin/mutgenomes_expected_mutrisk.R` (836 lines with FIXME), `bin/mutrate_genome_trinuc_corrected.R`, `bin/signatures_msighdp_run.R` have no test harness. +- Files: `bin/*.R` +- Risk: Statistical errors in expected-mutation-risk and dN/dS computations go unnoticed. +- Priority: Medium. + +--- + +*Concerns audit: 2026-10-09* diff --git a/.planning/codebase/CONVENTIONS.md b/.planning/codebase/CONVENTIONS.md new file mode 100644 index 00000000..f7066361 --- /dev/null +++ b/.planning/codebase/CONVENTIONS.md @@ -0,0 +1,156 @@ +--- +last_mapped_commit: 3c83c5b8aa313077e0ce43239a7b818b281b32ea +last_mapped_at: 2026-10-09 +--- +# Coding Conventions + +**Analysis Date:** 2026-10-09 + +## Language Mix + +The codebase is a Nextflow (DSL2) pipeline with three script languages in `bin/`: + +- **Python** (~76 scripts) — dominant language for analysis/plotting logic +- **R** (4 scripts) — `bin/dNdS_run.R`, `bin/mutgenomes_expected_mutrisk.R`, `bin/mutrate_genome_trinuc_corrected.R`, `bin/signatures_msighdp_run.R` +- **Groovy/Nextflow** — `main.nf`, `workflows/`, `subworkflows/local/`, `modules/local/` +- **Shell** (1 script) — `bin/createcustombed.sh` + +## Naming Patterns + +**Files (Python):** +- `snake_case.py` for all scripts: `check_samplesheet.py`, `mut_profile.py`, `plot_depths.py` +- Shared helper modules prefixed `utils_` or suffixed `_utils`: `bin/read_utils.py`, `bin/utils.py`, `bin/utils_plot.py`, `bin/utils_filter.py`, `bin/utils_context.py`, `bin/utils_impacts.py` +- Test files: `test_.py` in `bin/test/` + +**Files (Nextflow):** +- Processes live in `modules/local//main.nf` (e.g. `modules/local/compute_profile/main.nf`) +- Subworkflows live in `subworkflows/local//main.nf` (e.g. `subworkflows/local/omega/main.nf`); multi-file groups use `subworkflows/local//main` as entry +- One workflow file: `workflows/deepcsa.nf` + +**Functions:** +- `snake_case`: `compute_mutation_matrix()`, `validate_and_transform()`, `check_samplesheet()` +- Private methods prefixed `_`: `RowChecker._validate_sample()`, `_validate_vcf_format()` (`bin/check_samplesheet.py`) +- Test helpers prefixed `_`: `_make_df()`, `_make_maf()` (`bin/test/test_utils_filter.py`) + +**Classes:** +- `PascalCase`: `RowChecker` (`bin/check_samplesheet.py`), `TestSampleNameValidation` (`bin/test/test_check_samplesheet.py`) + +**Variables/Constants:** +- `snake_case` locals and module-level config: `custom_na_values` (`bin/read_utils.py`) +- `UPPER_SNAKE_CASE` for module constants: `MAX_CATEGORIES_PER_PLOT`, `MAX_HEATMAP_CELLS` (`bin/plot_depths.py`), `MIN_NONZERO_PVALUE` (`bin/utils.py`), `VALID_FORMATS_BAM`/`VALID_FORMATS_VCF` (`bin/check_samplesheet.py`) + +**Nextflow processes:** +- `UPPERCASE` process names: `COMPUTE_PROFILE`, `INPUT_CHECK` +- Subworkflows `UPPERCASE` with aliases on import: `MUTATIONAL_PROFILE as MUTPROFILEALL` (`workflows/deepcsa.nf`) + +## Code Style + +**Formatting:** +- Black configured in `pyproject.toml`: `line-length = 120`, `target_version = ["py37", "py38", "py39", "py310"]` +- isort with `profile = "black"` (`pyproject.toml`) +- `.editorconfig`: 4-space indent, LF line endings, UTF-8, final newline, trimmed trailing whitespace; 2-space for `*.md`, `*.yml`, `*.yaml`, `*.html`, `*.css`, `*.scss`, `*.js` +- **Reality check:** only `bin/check_samplesheet.py` is consistently Black-formatted (it is nf-core template code). Most other `bin/` scripts predate the config and use pandas-style `=` spacing (`sep = "\t"`), inconsistent quoting, and occasional 2-space indents. **Do not reformat existing scripts wholesale; match the style of the file you are editing. New nf-core-template-derived code must be Black-formatted.** + +**Linting:** +- No repo-level flake8/ruff/pylint config. `.devcontainer/devcontainer.json` enables pylint/flake8 in the dev container only. +- `markdownlint` config in `.markdownlint.json`: `MD013` (line length) off, `MD024` siblings-only. + +## Shebangs & Execution Model + +- Python: `#!/usr/bin/env python3` (nf-core template scripts) or `#!/usr/bin/env python` (older scripts) +- R: `#!/opt/conda/bin/Rscript --vanilla` (hard-coded conda path — matches the `bbglab/deepcsa-core` container) +- Scripts in `bin/` are invoked by Nextflow processes directly by name (Nextflow auto-stages `bin/` onto `PATH`), e.g. `mut_profile.py profile --sample_name ...` in `modules/local/compute_profile/main.nf` + +## CLI Argument Parsing + +**Two coexisting styles — prefer `click` for new scripts:** + +1. **click (dominant, newer scripts):** decorators with `click.Choice`, `click.Path(exists=True)` for input validation, `is_flag=True` for booleans. Example: `bin/mut_profile.py` (`@click.command()`, `@click.argument('mode', type=click.Choice(['matrix', 'profile']))`), `bin/plot_depths.py` +2. **argparse (nf-core template scripts):** `parse_args(argv=None)` function + `main(argv=None)` + `if __name__ == "__main__: sys.exit(main())`. Example: `bin/check_samplesheet.py` + +**Rules for new scripts:** +- Use `click.Path(exists=True)` so missing inputs fail fast +- Use `click.Choice` for enumerated modes +- Keep a `main()` entry point guarded by `if __name__ == '__main__':` (71 of 80 Python scripts do this) so functions stay importable by tests + +## Import Organization + +**Order (observed in `bin/check_samplesheet.py`, `bin/mut_profile.py`):** +1. Standard library (`sys`, `argparse`, `csv`, `logging`, `re`, `pathlib`) +2. Third-party (`click`, `pandas`, `numpy`, `matplotlib`, `seaborn`, `polars`) +3. Local sibling modules (`from utils import contexts_formatted`, `from read_utils import custom_na_values`) + +**Path Aliases:** +- None. Sibling imports rely on Nextflow staging all of `bin/` into the task working directory, so `import utils` works at runtime. Unit tests replicate this with `sys.path.insert(0, str(Path(__file__).parent.parent))` (`bin/test/test_utils_filter.py`). + +## Error Handling + +**Patterns:** +- **Validation scripts (nf-core template style):** raise `AssertionError` with descriptive f-string messages from validator methods; catch in the caller, log with `logger.critical(...)`, and `sys.exit(1)` (`bin/check_samplesheet.py`). Missing input file → `sys.exit(2)`. +- **click scripts:** rely on `click.Path(exists=True)` for input validation; domain errors often just `print(...)` + `exit(1)` (e.g. `bin/mut_profile.py` `compute_mutation_matrix`). +- **Fail-fast over recovery:** pipeline scripts generally do not catch exceptions; Nextflow's `errorStrategy = 'retry'` with `maxRetries = 2` (`tests/nextflow.config`) handles transient failures at the process level. +- **Defensive data handling:** explicit `fillna(0)` after reindexing, percentile clipping of outliers (99.5th percentile of `ALT_DEPTH` in `bin/mut_profile.py`), and `custom_na_values` lists passed to `pd.read_csv` (`bin/read_utils.py`). + +**Do this for new validation code:** raise `AssertionError`/`ValueError` with a message that names the offending value and the allowed set, and let the CLI wrapper convert to a non-zero exit. + +## Logging + +**Framework:** mixed — no single standard. + +- `logging` module with `logging.basicConfig(level=args.log_level, format="[%(levelname)s] %(message)s")` and a module-level `logger = logging.getLogger()` — `bin/check_samplesheet.py` (nf-core template pattern) +- `click.echo()` for user-facing progress in click scripts (`bin/mut_profile.py`) +- Plain `print()` in most analysis scripts, sometimes with a `[scriptname]` prefix: `print(f"[plot_depths] Skipping plot section '{section}': {reason}")` (`bin/plot_depths.py`) + +**When to log:** skip conditions with reasons, percentile/clip decisions, mode and parameter echo at startup. Keep messages one-line; Nextflow captures stdout/stderr per task. + +## Comments + +**When to Comment:** +- Section banners in test/Nextflow files using `/* ==== ... ==== */` blocks (`tests/deepcsa.nf.test`, `main.nf`) +- `# TODO` / `# FIXME` markers are used liberally (~20 occurrences across `bin/`), e.g. `bin/concat_sbs_probs.py:3` ("TODO: bump pandas to 2.2.3"), `bin/postprocessing_annotation.py:146` +- Explanatory comments for non-obvious numeric choices (plot size limits in `bin/plot_depths.py`) + +**Docstrings:** +- **Google style** in nf-core template code: `Args:`, `Returns:`, `Attributes:`, `Raises:` sections, plus `Example:` blocks (`bin/check_samplesheet.py`) +- **NumPy style** in bbglab-written code: `Parameters\n----------`, `Returns\n-------` (`bin/utils.py`, `bin/test/test_mask_matrix.py`) +- Module-level docstrings summarizing purpose (`"""Provide a command line tool to validate and transform tabular samplesheets."""`) +- **For new code:** match the host file's style; Google style for anything derived from nf-core templates. + +## Function Design + +**Size:** no enforced limit; legacy scripts contain very long functions and files (`bin/plot_gene_saturation.py` is 1546 lines). New code should keep functions single-purpose. + +**Parameters:** plain positional args for 2–3 params; keyword args with defaults for options (`def check_samplesheet(file_in, file_out, bam_required=False)`). R scripts use `optparse` long flags (`--inputfile`, `--outputprefix`). + +**Return Values:** DataFrames in/out for data-transform functions; `None` + file writes for plot/report functions. Functions that can fail return `None` and let the caller decide (`compute_mutation_profile` returning `None` in `bin/mut_profile.py`). + +## Module Design + +**Exports:** flat modules; no `__all__`. Tests import named functions directly: `from utils_filter import filter_maf, somatic_mask` (`bin/test/test_utils_filter.py`). + +**Barrel Files:** none. Shared constants live in dedicated modules (`bin/read_utils.py` for NA values / MAF reading, `bin/utils.py` for MAF filters and variant typing, `bin/utils_plot.py` for plotting helpers). + +## Nextflow-Specific Conventions + +**Process definition** (`modules/local/compute_profile/main.nf`): +- `tag "$meta.id"` for traceability +- Resource `label`s: `cpu_low`/`mem_low`/`process_high_memory` plus domain label `deepcsa_core` (container pinned in `conf/modules.config` via `withLabel: deepcsa_core { container = "docker.io/bbglab/deepcsa-core:0.1.0" }`) +- `input:`/`output:`/`script:`/`stub:` blocks; outputs use `emit:` names and `optional:true` where conditional +- `versions.yml` written via heredoc `cat <<-END_VERSIONS` and emitted `topic: versions` in every process +- `task.ext.args` / `task.ext.prefix` consumed with `?: ""` defaults for configurability from `conf/modules.config` +- `stub:` blocks present in 97 of 48 module dirs' main.nf files (widespread) for fast testing + +**Config layering:** `nextflow.config` → `conf/base.config`, `conf/modules.config`, `conf/results_outputs.config`, tool configs (`conf/tools/*.config`), profile configs (`conf/test.config`, `conf/exome.config`, `conf/mice.config`, ...). Site-specific paths isolated in `conf/general_files_IRB.config`. + +**Publishing:** output routing centralized in `conf/results_outputs.config` with `saveAs` filters that drop `versions.yml` from published dirs. + +## R Conventions + +- Shebang `#!/opt/conda/bin/Rscript --vanilla` +- `optparse` for CLI (`make_option(c("-n", "--samplename"), ...)`) — `bin/dNdS_run.R` +- `snake_case` function names (`is_SNV` is a legacy exception) +- Usage example in a header comment block + +--- + +*Convention analysis: 2026-10-09* diff --git a/.planning/codebase/INTEGRATIONS.md b/.planning/codebase/INTEGRATIONS.md new file mode 100644 index 00000000..ab1d268e --- /dev/null +++ b/.planning/codebase/INTEGRATIONS.md @@ -0,0 +1,133 @@ +--- +last_mapped_commit: 3c83c5b8aa313077e0ce43239a7b818b281b32ea +last_mapped_at: 2026-10-09 +--- +# External Integrations + +**Analysis Date:** 2026-10-09 + +## APIs & External Services + +**None at runtime.** The pipeline is fully offline: VEP runs with `--cache --offline` (`params.vep_params` in `nextflow.config`), and `BGDATA_OFFLINE=TRUE` is forced in the `env` block. No HTTP calls, no cloud SDKs, no database drivers anywhere in `bin/` or module scripts. + +**External data resources (consumed as files, not APIs):** + +| Resource | Purpose | Default param (`nextflow.config`) | Site override | +| --- | --- | --- | --- | +| Ensembl VEP cache (v111, GRCh38, homo_sapiens) | Variant annotation | `vep_cache = ".vep"` (staged into process) | `conf/general_files_IRB.config`: `/data/bbg/datasets/vep` | +| COSMIC v3.4 SBS signatures (GRCh38) | Signature fitting | `cosmic_ref_signatures` | IRB: COSMIC v3.5 at `/data/bbg/datasets/COSMIC_signatures/` | +| COSMIC v3.4 ID signatures (GRCh37) | Indel signature fitting | `indel_ref_signatures` | IRB: v3.5 | +| CADD v1.7 whole-genome SNV scores (+ .tbi) | OncodriveFML | `cadd_scores`, `cadd_scores_ind` | IRB: `/data/bbg/datasets/CADD/v1.7/hg38/` | +| Genome FASTA (GRCh38 no-alt masked) | Reference sequence | `fasta` (user-supplied) | IRB: `/data/bbg/datasets/genomes/GRCh38/...` | +| dNdScv Biomart reference (Ensembl v111) | dN/dS CDS mapping | `dnds_biomart_ref` | IRB: MANE version | +| dNdScv covariates (epigenome/PCAWG, .rda) | dN/dS covariates | `dnds_covariates` | IRB path | +| Omega covariates TSV (hg19/hg38 epigenome PCAWG) | Covariate-aware Omega | `omega_covariates_cov_file` (default ships in `assets/omega-covariates/`) | IRB path | +| Oncodrive3D datasets + annotations | 3D selection | `datasets3d`, `annotations3d` | IRB: dated snapshots (`datasets-260603`) | +| NanoSeq SNP/noise masks (BED) | Artifact filtering | `nanoseq_snp`, `nanoseq_noise` | IRB paths | +| GFF3 (Homo_sapiens.GRCh38.111) | DNA→protein/domain mapping | `gff3_file` | IRB path | +| Whole-genome trinucleotide counts | Mutability normalization | `wgs_trinuc_counts` (ships in `assets/trinucleotide_counts/`) | IRB path | +| bbgdomains annotated TSV | Subgenic/domain regions | `domains_file` (IRB only) | IRB path | +| gnomAD allele frequencies | Germline filtering | Embedded in VEP annotation (`--af_gnomadg --af_gnomade`), threshold `gnomad_af_threshold = 1e-3` | — | + +Consumption pattern: `channel.fromPath(params.X, checkIfExists: true)` in `workflows/deepcsa.nf` and `subworkflows/local/*/main.nf` — all file-based, no network fetch. + +## Data Storage + +**Databases:** +- None. All I/O is TSV/MAF/VCF/BED/PDF files on a shared filesystem. + +**File Storage:** +- Local/shared filesystem only. Outputs published under `params.outdir` with routing defined in `conf/modules.config` (`${params.outdir}/processing_files//`) and `conf/results_outputs.config` (curated outputs: `depths/`, `plots/`, `qc/`, `selection/`, `mutdensity/`, `signatures/`, …). `publish_dir_mode = 'copy'` by default. +- Some cohort-level aggregates written directly via `collectFile(storeDir: "${params.outdir}/...")` (e.g. `workflows/deepcsa.nf:341`, `subworkflows/local/dnds/main.nf:39-45`). + +**Caching:** +- Nextflow work dir (default `work/`); nf-test work dir `.nf-test` or `$DEEPCSA_TEST_WORKDIR` (`nf-test.config`). +- Singularity image cache: `singularity.cacheDir`/`libraryDir` = `/data/bbg/datasets/pipelines/nextflow_containers` (`conf/general_files_IRB.config`, `tests/nextflow.config`). +- No application-level cache service (no Redis/Memcached). + +## Authentication & Identity + +**Auth Provider:** None. No authentication anywhere — pipeline runs on trusted HPC infrastructure. Notification endpoints are the only externally-addressable values (see below). + +## Monitoring & Observability + +**Error Tracking:** +- None (no Sentry etc.). Failure handling is Nextflow-native: retry with exponential backoff on exit codes 104/130–145, `maxRetries 3` (`conf/base.config`); labels `error_ignore`/`error_retry`. + +**Logs / Reports:** +- Nextflow built-ins enabled in `nextflow.config`: `timeline`, `report`, `trace`, `dag` → timestamped files under `${params.outdir}/pipeline_info/`. +- **MultiQC** aggregates per-process `versions.yml` (emitted with `topic: versions` by ~76 local modules) plus custom `workflow_summary_mqc.yaml` and `methods_description_mqc.yaml` (`workflows/deepcsa.nf:762-790`, config `assets/multiqc_config.yml`). +- `CUSTOM_DUMPSOFTWAREVERSIONS` (`modules/nf-core/custom/dumpsoftwareversions`) captures the full software environment. + +**Notifications (outgoing):** +- **Email** — completion/failure email via local `sendmail`/`mail` binaries, HTML rendered from `assets/sendmail_template.txt` / `assets/email_template.html` (`subworkflows/nf-core/utils_nfcore_pipeline/main.nf:296-318`). Triggered by `params.email` / `params.email_on_fail`. +- **Slack / MS Teams** — `imNotification` posts to `params.hook_url`; Slack format (`assets/slackreport.json`) when URL contains `hooks.slack.com`, otherwise Adaptive Cards (`assets/adaptivecard.json`) for Teams (`subworkflows/nf-core/utils_nfcore_pipeline/main.nf:403-404`). + +## CI/CD & Deployment + +**Hosting:** +- Source: GitHub (`https://github.com/bbglab/deepCSA`, manifest in `nextflow.config`). DOI: `dx.doi.org/10.17504/protocols.io.dm6gp1jodgzp/v2`. +- No `.github/` directory — no GitHub Actions workflows detected in this checkout. + +**CI Pipeline:** +- **nf-test** is the test harness (`nf-test.config`, `tests/deepcsa.nf.test` + snapshot). Per `docs/test_data.md` and `tests/nextflow.config`, tests execute on the IRB SLURM cluster (queues `bbg_cpu_zen4,irb_cpu_zen4`) with Singularity; results CSVs (`tests/2026-10-03_results.csv`, `tests/2026-10-09_results.csv`) are committed. No hosted CI runner config found in-repo. +- Container Dockerfile recipes live out-of-repo: `https://github.com/bbglab/containers-recipes` (`docs/tools.md`). + +**Deployment model:** +- `nextflow run bbglab/deepCSA -profile ,` — no packaging, no container registry push automation in-repo. + +## Environment Configuration + +**Required env vars:** +- None mandatory. Optional: `DEEPCSA_TEST_WORKDIR` (nf-test work dir), `HOME=/tmp` and the isolation vars (`PYTHONNOUSERSITE`, `R_PROFILE_USER`, `R_ENVIRON_USER`, `JULIA_DEPOT_PATH`, `BGDATA_OFFLINE`) are set by the pipeline itself (`nextflow.config` `env` block). + +**Required params (schema-enforced, `nextflow_schema.json`):** +- `input` (samplesheet CSV, validated against `assets/schema_input.json`) and `outdir` are the only `required` entries. `fasta` required unless reference paths come from a site profile. `input_maf` alternative input mode requires `use_custom_depths = true` (`workflows/deepcsa.nf:189-201`). + +**Secrets location:** +- No secrets in repo. Notification hook URLs are passed at runtime via `--hook_url`; email addresses via `--email`. No `.env` files present. + +## Webhooks & Callbacks + +**Incoming:** +- None. + +**Outgoing:** +- Slack webhook (`hooks.slack.com`) or MS Teams Adaptive Card POST, only when `params.hook_url` is set (`subworkflows/nf-core/utils_nfcore_pipeline/main.nf`). +- Local `sendmail` invocation (not an HTTP callback). + +## Container Images (external registries) + +All pulled from Docker Hub / biocontainers at task launch; registry override via `docker.registry`/`singularity.registry` (default `quay.io` in `nextflow.config`, but module directives hard-code `docker.io`/`biocontainers`). + +| Image | Used by | Defined in | +| --- | --- | --- | +| `docker.io/bbglab/deepcsa-core:0.1.0` | ~60 `deepcsa_core`-labeled Python processes | `conf/modules.config:714` | +| `docker.io/bbglab/deepcsa_bed:latest` | Panel consensus (pybedtools/polars) | `modules/local/createpanels/consensus/main.nf` | +| `docker.io/bbglab/omega:0.2.1` | Omega preprocess/mutabilities/estimator | `modules/local/bbgtools/omega/*/main.nf` | +| `docker.io/ferriolcalvet/omegacovariates:v0.1.0` | Covariate-aware Omega | `modules/local/omega_covariates/run/main.nf` | +| `docker.io/spellegrini87/oncodrive3d:1.0.9-light` / `-chimerax` | Oncodrive3D run/plots | `modules/local/bbgtools/oncodrive3d/*/main.nf` | +| `docker.io/ferriolcalvet/oncodrivefml:latest` | OncodriveFML | `modules/local/bbgtools/oncodrivefml/main.nf` | +| `docker.io/ferriolcalvet/oncodriveclustl:latest` | OncodriveCLUSTL | `modules/local/bbgtools/oncodriveclustl/main.nf` | +| `docker.io/ferriolcalvet/dnds:latest` | dNdScv buildref/run | `modules/local/dnds/*/main.nf` | +| `docker.io/ferriolcalvet/sigprofiler_assignment:1.1.3` | Signature fitting | `modules/local/signatures/sigprofiler/assignment/*/main.nf` | +| `docker.io/ferriolcalvet/sigprofilermatrixgenerator:1.3.5` | SBS96 matrix generation | `modules/local/signatures/sigprofiler/matrixgenerator/main.nf` | +| `docker.io/ferriolcalvet/sigprofilerassignment` | SigProfilerExtractor | `modules/local/signatures/sigprofiler/extractor/main.nf` | +| `docker.io/ferriolcalvet/msighdp:latest`, `hdp_wrapper`, `musical:latest` | HDP/MUSICAL signature extraction | `modules/local/signatures/*/main.nf` | +| `docker.io/rblancomi/bbgregressions:dev` | bbgregressions (label-based) | `conf/modules.config:719` | +| `docker.io/ferriolcalvet/runningr:v1` | R mutation-density scaling | `modules/local/mut_density/wgscaled/main.nf` | +| `docker.io/ferriolcalvet/saturation:v0.1.0` | Saturation kinetics | `modules/local/saturation_kinetics/compute/main.nf` | +| `docker.io/axelrosendahlhuber/expected_mutrate:latest` | Expected mutated cells | `modules/local/mutated_cells_expected/main.nf` | +| `docker.io/ferranmuinos/test_mutated_genomes` | Mutated genomes from VAF | `modules/local/mutated_genomes_from_vaf/main.nf` | +| `biocontainers/ensembl-vep:111.0--pl5321h2a3209d_0` (also 102/108 variants) | VEP annotation | `modules/nf-core/ensemblvep/*/main.nf` | +| `biocontainers/samtools:1.18--h50ea8bc_1` | Depth computation | `modules/local/computedepths/main.nf` | +| `biocontainers/tabix:1.11--hdfd78af_0` | Indexed TSV queries | `modules/nf-core/tabix/*/main.nf` | +| `biocontainers/multiqc:1.20--pyhdfd78af_0` | QC report | `modules/nf-core/multiqc/main.nf` | +| `biocontainers/pybedtools:0.9.1--py38he0f268d_0` | Custom BED handling | `modules/local/createpanels/custombedfile/main.nf` | +| `biocontainers/python:3.8.3` | Samplesheet check | `modules/local/samplesheet_check.nf` | + +Conda fallbacks exist for a handful of processes (`-profile conda`/`mamba`): inline specs in `modules/local/createpanels/*/main.nf` (`python=3.10.17`, `pybedtools=0.12.0`, `polars=1.30.0`, `click=8.2.1`, `gcc_linux-64=15.1.0`) and `modules/local/computedepths/environment.yml` (`bioconda::samtools=1.18`). + +--- + +*Integration audit: 2026-10-09* diff --git a/.planning/codebase/STACK.md b/.planning/codebase/STACK.md new file mode 100644 index 00000000..d935b6eb --- /dev/null +++ b/.planning/codebase/STACK.md @@ -0,0 +1,108 @@ +--- +last_mapped_commit: 3c83c5b8aa313077e0ce43239a7b818b281b32ea +last_mapped_at: 2026-10-09 +--- +# Technology Stack + +**Analysis Date:** 2026-10-09 + +## Languages + +**Primary:** +- **Nextflow (DSL2)** — pipeline orchestration. `main.nf` (entry point, `nextflow.enable.dsl = 2`), `workflows/deepcsa.nf` (main `DEEPCSA` workflow, ~800 lines), `subworkflows/local/**` (22 subworkflows), `modules/local/**` (~60 local modules). Manifest requires `nextflowVersion = '!>=25.04.2'` (`nextflow.config`). +- **Python 3** — analysis scripts in `bin/*.py` (~90 scripts). CLI via `click`, data via `pandas`/`polars`, plotting via `matplotlib`/`seaborn`, stats via `scipy`/`statsmodels`. Container images run Python 3.10.x (e.g. conda spec `python=3.10.17` in `modules/local/createpanels/captured/main.nf`); `modules/local/samplesheet_check.nf` pins Python 3.8.3. `pyproject.toml` targets py37–py310 for lint config only. + +**Secondary:** +- **R** — 4 scripts: `bin/dNdS_run.R` (dNdScv), `bin/mutgenomes_expected_mutrisk.R`, `bin/mutrate_genome_trinuc_corrected.R`, `bin/signatures_msighdp_run.R` (mSigHdp). Bioconductor-heavy (see R libraries below). +- **Groovy** — implicit in Nextflow; used directly in `subworkflows/nf-core/utils_nfcore_pipeline/main.nf` (email/notification logic, `java.io.File`, Groovy template engine). +- **Bash/Shell** — process `script:` blocks throughout modules; `bin/createcustombed.sh`. +- **YAML** — regression configs (`assets/regressions/configs/*.yml`), parsed by `bin/regressions_editconfig.py` (PyYAML). + +## Runtime + +**Environment:** +- Nextflow >= 25.04.2 on JVM (Java 17+ implied by Nextflow 25.x) +- Nextflow plugin: `nf-schema@2.3.0` (declared in `plugins` block of `nextflow.config`) — used for params validation (`paramsSummaryMap` in `workflows/deepcsa.nf`) and schema-based input checks +- Execution engines (profiles in `nextflow.config`): Docker, Singularity, Podman, Shifter, Charliecloud, Apptainer, Conda, Mamba. Default container registry set to `quay.io` but all module `container` directives point at `docker.io` / `biocontainers` explicitly. + +**Package Manager:** +- **No application-level package manager / lockfile.** Dependencies are pinned per-process via `container` directives in module files and inline `conda` directives (e.g. `modules/local/computedepths/main.nf` → `environment.yml` with `bioconda::samtools=1.18`). +- `pyproject.toml` is **lint config only** (Black line-length 120, isort black profile) — not an installable package. +- nf-test plugins: `nft-utils@0.0.3` loaded in `nf-test.config`. + +## Frameworks + +**Core:** +- **nf-core pipeline template conventions** — boilerplate subworkflows `subworkflows/nf-core/utils_nextflow_pipeline`, `utils_nfcore_pipeline`, `utils_nfvalidation_plugin`; `PIPELINE_INITIALISATION` in `subworkflows/local/utils_nfcore_deepcsa/main.nf` (banner, params summary, completion email/notifications). +- **nf-core modules (pinned via `modules.json`, branch master, git_sha `911696ea0b62df80e900ef244d7867d177971f73`):** + - `modules/nf-core/custom/dumpsoftwareversions` — software version capture for MultiQC + - `modules/nf-core/ensemblvep` (`vep`, `veppanel`) — variant annotation + - `modules/nf-core/multiqc` — aggregated QC report + - `modules/nf-core/tabix` (`bgziptabix`, `bgziptabixquery`) — indexed TSV querying, reused under many aliases (QUERYDEPTHS, QUERYPANEL, DEPTHS.*CONS) +- **nf-core subworkflows:** `subworkflows/nf-core/vcf_annotate_ensemblvep`, `vcf_annotate_ensemblvep_panel` + +**Testing:** +- **nf-test** — `nf-test.config` (testsDir `tests`, workDir `.nf-test` or `$DEEPCSA_TEST_WORKDIR`, profile `test,singularity`, ignores `modules/nf-core/**` and `subworkflows/nf-core/**`). Main test: `tests/deepcsa.nf.test` with snapshot file `tests/deepcsa.nf.test.snap`. Per `docs/test_data.md`, tests must run on a SLURM cluster with Singularity. + +**Build/Dev:** +- No build step. Linting: Black + isort via `pyproject.toml` (applies to `bin/check_samplesheet.py` per file comment). +- Process resource/error policy centralized in `conf/base.config` (labels `cpu_low`/`cpu_medium`/`cpu_high`/`mem_low`, `error_ignore`, `error_retry`; exponential-backoff retry on exit codes 104, 130–145). + +## Key Dependencies + +**Critical (Python, imported across `bin/`):** +- `click` — CLI framework for nearly all `bin/*.py` scripts +- `pandas` — tabular mutation/depth/panel data (`bin/utils.py`, `bin/read_utils.py`, `bin/filter_cohort.py`, …) +- `polars` — high-performance panel processing (`bin/create_consensus_panel.py`, `bin/create_panel_versions.py`; conda spec `conda-forge::polars=1.30.0`) +- `matplotlib`, `seaborn` — all PDF plotting (`bin/utils_plot.py`, `bin/plot_*.py`) +- `scipy` — statistical tests (`bin/compute_hotspots_selection.py`, `bin/signatures_musical.py`) +- `statsmodels` — covariate regression models (`bin/omega_covariates_definitions.py`) +- `PyYAML` — regression config editing (`bin/regressions_editconfig.py`) +- `pybedtools` — BED interval ops (in `createpanels` containers) + +**Critical (R):** +- `dndscv` — dN/dS selection model (`bin/dNdS_run.R`) +- `BSgenome.Hsapiens.UCSC.hg38`, `GenomicRanges`, `IRanges`, `Biostrings` — genome sequence/ranges (`bin/mutgenomes_expected_mutrisk.R`) +- `tidyverse`, `data.table`, `ggplot2`, `jsonlite`, `optparse`, `Hmisc`, `R.utils`, `abind` +- `ICAMS`, `mSigHdp` — signature extraction (`bin/signatures_msighdp_run.R`) + +**Infrastructure (bioinformatics tools, containerized):** +- Ensembl VEP 111 (default; also supports 102/108 via `params.vep_cache_version` conditional containers in `modules/nf-core/ensemblvep/*/main.nf`) +- samtools 1.18 (`modules/local/computedepths`) +- tabix 1.11 (`modules/nf-core/tabix`) +- MultiQC 1.20 (`modules/nf-core/multiqc`) +- bbglab method containers — see INTEGRATIONS.md + +## Configuration + +**Environment:** +- All runtime config via Nextflow params — no `.env` files. Schema-validated by `nextflow_schema.json` (nf-schema) and `assets/schema_input.json` (samplesheet columns). +- Hardened env block in `nextflow.config`: `PYTHONNOUSERSITE=1`, `R_PROFILE_USER=/.Rprofile`, `R_ENVIRON_USER=/.Renviron`, `JULIA_DEPOT_PATH=/usr/local/share/julia`, `BGDATA_OFFLINE=TRUE`, `HOME=/tmp` — prevents host Python/R libraries leaking into containers. +- Site-specific reference paths in `conf/general_files_IRB.config` (IRB cluster: `/data/bbg/datasets/...` for VEP cache, FASTA, COSMIC, CADD, dNdScv, Oncodrive3D, NanoSeq masks, GFF3) + `singularity.cacheDir`/`libraryDir`. +- Test env var: `DEEPCSA_TEST_WORKDIR` (overrides nf-test work dir, `nf-test.config`). + +**Build/Run config files:** +- `nextflow.config` — params defaults, profiles, plugins, env, timeline/report/trace/dag reporting, manifest (name `bbglab/deepCSA`, version `1.0.1.dev`) +- `conf/base.config` — global resource limits (`max_memory 950.GB`, `max_cpus 196`, `max_time 30.d`), CPU/memory tiers, retry strategy, per-process `withName` overrides +- `conf/modules.config` — `ext.args` per module, `publishDir` routing (`${params.outdir}/processing_files/...`), label→container mapping (`deepcsa_core` → `docker.io/bbglab/deepcsa-core:0.1.0`, `bbgregressions` → `docker.io/rblancomi/bbgregressions:dev`); includes `conf/tools/*.config` and `conf/results_outputs.config` +- `conf/tools/` — tool-specific params: `panels.config`, `omega.config`, `mutdensity.config`, `oncodrive3d.config`, `oncodrivefml.config`, `hdp_sig_extraction.config`, `regressions.config` +- `conf/modes/` — analysis-mode profiles: `basic.config`, `clonal_structure.config`, `get_signatures.config` +- Other profiles: `test`, `test_real`, `test_regressions`, `mice`, `exome`, `irbcluster`, `local`, `debug`, `gitpod`, plus one per container engine +- `nextflow_schema.json` — full parameter schema (input/output, genome, grouping, hotspot, etc.) +- `tower.yml` — Seqera Platform report display config + +## Platform Requirements + +**Development:** +- Nextflow >= 25.04.2 + JVM +- A container engine (Docker locally; Singularity on cluster) or Conda/Mamba +- For tests: SLURM cluster with Singularity, queues `bbg_cpu_zen4,irb_cpu_zen4` (`tests/nextflow.config`); local execution not supported per `docs/test_data.md` + +**Production:** +- HPC via SLURM (IRB cluster profile `irbcluster` / `tests/nextflow.config`); any Nextflow executor works in principle (no executor hard-coded outside test/local configs) +- Shared filesystem for reference data (`/data/bbg/datasets/...`) and Singularity image cache (`/data/bbg/datasets/pipelines/nextflow_containers`) +- Resource ceiling: 196 CPUs, 950 GB RAM, 30 days per process (defaults in `nextflow.config`) + +--- + +*Stack analysis: 2026-10-09* diff --git a/.planning/codebase/STRUCTURE.md b/.planning/codebase/STRUCTURE.md new file mode 100644 index 00000000..d2ed6995 --- /dev/null +++ b/.planning/codebase/STRUCTURE.md @@ -0,0 +1,234 @@ +--- +last_mapped_commit: 3c83c5b8aa313077e0ce43239a7b818b281b32ea +last_mapped_at: 2026-10-09 +--- +# Codebase Structure + +**Analysis Date:** 2026-10-09 + +## Directory Layout + +``` +deepCSA/ +├── main.nf # Entry point (BBGTOOLS wrapper workflow) +├── nextflow.config # Params defaults, profiles, registries +├── nextflow_schema.json # nf-schema param validation +├── modules.json # Pinned nf-core module/subworkflow versions +├── nf-test.config # nf-test runner config +├── tower.yml # Seqera Platform report display config +├── pyproject.toml # Black/isort config for bin/ scripts +├── workflows/ # Pipeline workflow layer +│ └── deepcsa.nf # THE pipeline (workflow DEEPCSA, ~800 lines) +├── subworkflows/ +│ ├── local/ # 23 local subworkflows (analysis stages) +│ │ ├── depthanalysis/ # Depth computation + filtering +│ │ ├── createpanels/ # Consensus panel creation +│ │ ├── mutationpreprocessing/ # VEP annotation, MAF building +│ │ ├── mutationdensity/ mutationprofile/ mutability/ +│ │ ├── omega/ oncodrivefml/ oncodrive3d/ oncodriveclustl/ dnds/ +│ │ ├── signatures/ signatures_hdp/ mutatedcells/ indels/ +│ │ ├── regressions/ enrichpanels/ adjmutdensity/ +│ │ ├── plotdepths/ plotting_qc/ plottingsummary/ +│ │ └── utils_nfcore_deepcsa/ # Init: banner, versions, methods text +│ └── nf-core/ # Vendored nf-core subworkflows (5) +├── modules/ +│ ├── local/ # 97 local process modules (48 top-level dirs) +│ │ ├── bbgtools/ # bbglab tool wrappers: omega/, oncodrive3d/, +│ │ │ # oncodrivefml/, oncodriveclustl/, +│ │ │ # bbgregressions/, sitecomparison/ +│ │ ├── createpanels/ # captured/, consensus/, compare/, custombedfile/ +│ │ ├── dnds/ signatures/ downsample/ plot/ ... +│ │ └── /main.nf # One process per main.nf +│ └── nf-core/ # Vendored nf-core modules (multiqc, dumpsoftwareversions) +├── bin/ # ~90 analysis scripts (auto on PATH in tasks) +│ ├── *.py # Python CLIs (click + pandas) +│ ├── *.R # R scripts (dNdS_run.R, mutrate_genome_trinuc_corrected.R) +│ ├── utils*.py # Shared helpers (utils.py, utils_plot.py, read_utils.py...) +│ ├── test/ # pytest unit tests (test_*.py) +│ └── saturation_mutagenesis/ # Saturation kinetics scripts + notebooks +├── conf/ # Nextflow config includes +│ ├── base.config # Resource limits, labels, error strategy +│ ├── modules.config # Per-process ext.args + publishDir (78 withName blocks) +│ ├── results_outputs.config # Curated final publish paths +│ ├── test.config, test_real.config, exome.config, mice.config, local.config +│ ├── modes/ # Analysis presets: basic, clonal_structure, get_signatures +│ └── tools/ # Per-tool configs: omega, oncodrive3d, oncodrivefml, +│ # mutdensity, panels, regressions, hdp_sig_extraction +├── assets/ # Static reference data & templates +│ ├── schema_input.json # Samplesheet JSON schema +│ ├── multiqc_config.yml, email_template.*, sendmail_template.txt, slackreport.json +│ ├── omega_consequences_groupings.json, chromosome_bands/, trinucleotide_counts/ +│ ├── omega-covariates/, regressions/, build_datasets/, assess_panel/ +│ ├── example_inputs/ # Example samplesheets/features tables +│ └── useful_scripts/ # Standalone helper scripts/notebooks (not run by pipeline) +├── docs/ # usage.md, output.md, metrics.md, tools.md, input_scenarios.md... +├── tests/ # nf-test pipeline tests +│ ├── deepcsa.nf.test # Main pipeline test +│ ├── deepcsa.nf.test.snap # Snapshot file +│ ├── nextflow.config # Test-specific config +│ └── test_data/ # Test fixtures (partially gitignored) +├── test_data/ # Module-level test fixtures +├── scratchhhh/ # Developer scratch (gitignored) +└── .planning/ # GSD planning docs (this directory) +``` + +## Directory Purposes + +**`workflows/`:** +- Purpose: Top-level pipeline workflows +- Contains: `deepcsa.nf` — the single `DEEPCSA` workflow with all channel wiring +- Key files: `workflows/deepcsa.nf` + +**`subworkflows/local/`:** +- Purpose: Reusable analysis stages composed of modules +- Contains: One `main.nf` per stage with `include` + `workflow X { take/main/emit }` +- Key files: `subworkflows/local/omega/main.nf`, `subworkflows/local/depthanalysis/main.nf`, `subworkflows/local/mutationpreprocessing/main.nf` + +**`subworkflows/nf-core/`:** +- Purpose: Vendored nf-core subworkflows, pinned via `modules.json` +- Contains: `utils_nextflow_pipeline`, `utils_nfcore_pipeline`, `utils_nfvalidation_plugin`, `vcf_annotate_ensemblvep`, `vcf_annotate_ensemblvep_panel` +- Key files: `subworkflows/nf-core/utils_nfcore_pipeline/main.nf` + +**`modules/local/`:** +- Purpose: Local process definitions (one `process` per `main.nf`) +- Contains: 97 `main.nf` files; `bbgtools/` sub-tree wraps bbglab tools (omega, oncodrive*, regressions) +- Key files: `modules/local/bbgtools/omega/estimator/main.nf`, `modules/local/createpanels/consensus/main.nf`, `modules/local/table2groups/main.nf` + +**`modules/nf-core/`:** +- Purpose: Vendored nf-core modules +- Contains: `multiqc`, `custom/dumpsoftwareversions` +- Key files: `modules/nf-core/multiqc/main.nf` + +**`bin/`:** +- Purpose: Executable analysis scripts; Nextflow adds this dir to task PATH automatically +- Contains: Python CLIs (click, pandas, pybedtools, polars), R scripts, shared `utils*.py` helpers, `test/` pytest suite, `saturation_mutagenesis/` extras +- Key files: `bin/utils.py`, `bin/read_utils.py`, `bin/mut_density_simple.py`, `bin/create_consensus_panel.py`, `bin/dNdS_run.R` + +**`conf/`:** +- Purpose: Nextflow configuration includes +- Contains: `base.config` (resources/labels/retry), `modules.config` (per-process `ext.args` + publishDir), `results_outputs.config` (final outputs), `tools/*.config`, `modes/*.config`, environment configs (`test.config`, `exome.config`, `mice.config`, `local.config`) +- Key files: `conf/base.config`, `conf/modules.config` + +**`assets/`:** +- Purpose: Static data, schemas, templates shipped with the pipeline +- Contains: samplesheet schema, MultiQC config, notification templates, reference datasets (trinucleotide counts, chromosome bands, omega covariates, regressions configs), example inputs +- Key files: `assets/schema_input.json`, `assets/omega_consequences_groupings.json`, `assets/placeholder_no_file.tsv` (used as empty-channel placeholder in `workflows/deepcsa.nf`) + +**`tests/` + `test_data/`:** +- Purpose: nf-test pipeline tests and fixtures +- Contains: `tests/deepcsa.nf.test` (+ snapshot), `tests/nextflow.config`, module fixtures in `test_data/modules/` +- Key files: `tests/deepcsa.nf.test`, `nf-test.config` + +**`docs/`:** +- Purpose: User/developer documentation +- Contains: `usage.md`, `output.md`, `metrics.md`, `tools.md`, `input_scenarios.md`, `file_formatting.md`, `issue_resolution.md`, `test_data.md`, `images/` + +## Key File Locations + +**Entry Points:** +- `main.nf`: Pipeline entry; runs `PIPELINE_INITIALISATION` then `BBGTOOLS` → `DEEPCSA` +- `workflows/deepcsa.nf`: All pipeline logic; the file to read to understand data flow + +**Configuration:** +- `nextflow.config`: Params defaults, container registries, profiles (docker/singularity/conda/test/debug) +- `nextflow_schema.json`: Param definitions/validation (nf-schema) +- `conf/base.config`: Resource labels (`cpu_low/medium/high`, `mem_low`), error strategy +- `conf/modules.config`: Per-process `ext.args`, `ext.prefix`, publishDir overrides +- `conf/results_outputs.config`: Curated user-facing output paths +- `conf/modes/*.config`: Preset param bundles (`basic`, `clonal_structure`, `get_signatures`) +- `conf/tools/*.config`: Per-tool param bundles +- `assets/schema_input.json`: Samplesheet schema + +**Core Logic:** +- `bin/*.py`: All computational logic (mutation density, profiles, panels, omega QC, plotting) +- `bin/utils.py`, `bin/utils_filter.py`, `bin/utils_impacts.py`, `bin/utils_context.py`, `bin/utils_plot.py`, `bin/read_utils.py`: Shared Python helpers — import these, don't duplicate +- `bin/*.R`: dNdScv and mutation-rate R analyses + +**Testing:** +- `tests/deepcsa.nf.test`: End-to-end pipeline test (nf-test + snapshots) +- `bin/test/test_*.py`: Python unit tests (pytest) +- `nf-test.config`: Test runner config (profile `test,singularity`, ignores `modules/nf-core/**`) + +## Naming Conventions + +**Files:** +- Nextflow modules/subworkflows: always `main.nf` inside a lowercase snake_case tool directory: `modules/local/bbgtools/omega/estimator/main.nf` +- Python scripts: `snake_case.py` matching the tool/analysis name: `create_consensus_panel.py`, `mut_density_adjusted.py` +- R scripts: `snake_case.R` with tool prefix: `dNdS_run.R`, `signatures_msighdp_run.R` +- Configs: lowercase with `.config`: `base.config`, `results_outputs.config` + +**Directories:** +- Module dirs: lowercase, no separators or underscores preferred: `createpanels`, `filterdepths`, `mutatedcells` (some legacy underscores: `omega_covariates`, `select_mutdensity`) +- bbglab tool wrappers grouped under `modules/local/bbgtools///` + +**Code symbols:** +- Nextflow processes: `UPPER_SNAKE_CASE` verbs/phrases: `CREATECONSENSUSPANELS`, `OMEGA_ESTIMATOR`, `COMPUTEDEPTHS` +- Workflows: `UPPER_SNAKE_CASE`: `DEEPCSA`, `OMEGA_ANALYSIS`, `DEPTH_ANALYSIS` +- Aliases when reusing: descriptive suffixes: `MUTDENSITYALL`, `MUTDENSITYPROT`, `OMEGAMULTI`, `PREPROCESSINGGLOBALLOC` + +## Where to Add New Code + +**New analysis stage (multi-process):** +- Subworkflow: `subworkflows/local//main.nf` — follow `subworkflows/local/depthanalysis/main.nf` structure (`take:`/`main:`/`emit:`) +- Wire it in `workflows/deepcsa.nf` with an alias include and a `params.`-gated call + +**New single process:** +- Module: `modules/local///main.nf` (use `bbgtools//` for bbglab tools) +- Declare `conda "..."` + `container 'docker://...'`, `tag "$meta.id"`, meta-tuple inputs, named `emit:` outputs, `versions.yml` with `topic: versions`, and a `stub:` block +- Add per-process `ext.args`/publishDir in `conf/modules.config`; final outputs in `conf/results_outputs.config` + +**New analysis script:** +- Python: `bin/.py` — use click for CLI, pandas for tables, import shared helpers from `bin/utils.py` / `bin/read_utils.py` (same dir, so plain `import` works) +- R: `bin/.R` +- The script becomes callable by name in any process script block (Nextflow puts `bin/` on PATH) + +**New params:** +- Add default in `nextflow.config` params block +- Add schema entry in `nextflow_schema.json` +- If tool-specific, also set in `conf/tools/.config` + +**New tests:** +- Pipeline-level: `tests/deepcsa.nf.test` (nf-test); module fixtures under `test_data/modules/` +- Python unit tests: `bin/test/test_.py` (pytest) + +**New static reference data:** +- `assets//` (e.g. `assets/trinucleotide_counts/`, `assets/omega-covariates/`) + +**New utilities:** +- Shared helpers: `bin/utils*.py` (extend existing files rather than creating new ones when the topic matches) + +## Special Directories + +**`bin/`:** +- Purpose: Executables auto-mounted on task PATH by Nextflow +- Generated: No +- Committed: Yes (except `bin/__pycache__/`, `bin/saturation_mutagenesis/` notebooks partially gitignored) + +**`modules/nf-core/`, `subworkflows/nf-core/`:** +- Purpose: Vendored nf-core components pinned by `modules.json` (nf-core tools install/update) +- Generated: Managed by `nf-core modules install` — do not hand-edit +- Committed: Yes + +**`assets/`:** +- Purpose: Static pipeline data/templates +- Generated: No (except `assets/HDP_files*` which is gitignored) +- Committed: Yes (with gitignored exceptions: `assets/useful_scripts/*.ipynb`, `assets/HDP_files*`) + +**`scratchhhh/`:** +- Purpose: Developer scratch outputs (SigProfiler runs, plots, test notes) +- Generated: Yes (ad hoc) +- Committed: No (gitignored) + +**`work/`, `results/`, `.nf-test/`, `testing/`:** +- Purpose: Nextflow/nf-test runtime outputs +- Generated: Yes +- Committed: No (gitignored) + +**`.planning/`:** +- Purpose: GSD planning documents (this analysis) +- Generated: By GSD commands +- Committed: Per team convention + +--- + +*Structure analysis: 2026-10-09* diff --git a/.planning/codebase/TESTING.md b/.planning/codebase/TESTING.md new file mode 100644 index 00000000..c26a23c0 --- /dev/null +++ b/.planning/codebase/TESTING.md @@ -0,0 +1,422 @@ +--- +last_mapped_commit: e7b4ed7374de1197b1c7bcde526d062730c9bd43 +last_mapped_at: 2026-10-09 +--- +# Testing Patterns + +**Analysis Date:** 2026-10-09 + +## Test Framework + +This repo has **three test layers**: + +1. **Pipeline-level integration tests** — nf-test `nextflow_pipeline` (whole `main.nf` on SLURM) +2. **Module-level process tests** — nf-test `nextflow_process` (single local module; currently only `expand_regions`) +3. **Script-level unit tests** — Python `unittest` (stdlib, no pytest) (legacy, these type of tests should not be added, nf-test tests are much preferred in most of the cases) + +**Runner (pipeline):** +- nf-test >= 0.9.2 with plugin `nft-utils@0.0.3` (loaded in `nf-test.config`) +- Config: `nf-test.config` (root) + `tests/nextflow.config` (cluster/executor config) +- Suite: `tests/deepcsa.nf.test` +- Snapshot store: `tests/deepcsa.nf.test.snap` + +**Runner (Python):** +- `unittest` via `python -m unittest` or `python -m pytest` (pytest-compatible; no pytest config exists) +- Tests live in `bin/test/`: `test_check_samplesheet.py`, `test_check_contamination.py`, `test_mask_matrix.py`, `test_plot_selectionsideplots.py`, `test_utils_filter.py` + +**Run Commands:** + +```bash + +# Pipeline integration tests (REQUIRES SLURM cluster + Singularity — do not run locally) + +nf-test test tests/deepcsa.nf.test # whole suite +nf-test test tests/deepcsa.nf.test --tag normal # single test by tag +nf-test test tests/deepcsa.nf.test --tag omega + +# Module-level process tests (single module, fast — NOT auto-discovered) + +nf-test test modules/local/expand_regions/tests/main.nf.test + +# Python unit tests (run anywhere) + +python -m unittest discover -s bin/test -v +python -m unittest bin/test/test_utils_filter.py # single file +python -m pytest bin/test/ -v # pytest alternative +``` + +**Work directory:** `DEEPCSA_TEST_WORKDIR` env var, default `.nf-test/` in repo root (`nf-test.config`). Test outputs live under `/tests//`. + +## Test File Organization + +**Location:** +- Pipeline tests: `tests/` (separate from code, committed) +- Python unit tests: `bin/test/` (separate `test/` subdir next to the scripts they test) + +**Naming:** +- nf-test: one suite file `tests/deepcsa.nf.test` containing all pipeline tests +- Python: `test_.py` mirroring the module name (`test_utils_filter.py` → `bin/utils_filter.py`) + +**Structure:** + +``` +tests/ +├── deepcsa.nf.test # nf-test suite (5 active tests + commented templates) +├── deepcsa.nf.test.snap # MD5 snapshots of pipeline outputs +├── nextflow.config # SLURM executor, Singularity, cluster paths +├── test_data/ # committed local inputs (input.csv, input_maf.csv, +│ # input_no_bam.csv, test_mutations.maf) +└── README.md # how to run, snapshot policy, debugging guide +modules/local/expand_regions/tests/ +├── main.nf.test # module-level nf-test (nextflow_process) +└── main.nf.test.snap # MD5 snapshots of process outputs +test_data/ # repo-root fixtures for module tests +├── dummy_file.tsv +└── modules/ # PPM1D BED/TSV fixtures for expand_regions +bin/test/ +└── test_*.py # unittest files +``` + +## Pipeline Tests (nf-test) + +**Suite Organization** (`tests/deepcsa.nf.test`): + +```groovy +nextflow_pipeline { + name "Test DEEPCSA Pipeline" + script "main.nf" + config "./nextflow.config" + + test("TEST 1. Basic functionality - MAF-based processing") { + tag "input_maf" + + when { + params { + input = "${projectDir}/tests/test_data/input_maf.csv" + outdir = "$outputDir" // nf-test built-in: unique temp dir per test + input_maf = 'https://raw.githubusercontent.com/bbglab/DeepClone_protocol/main/...' + use_custom_depths = true + profileall = true + signatures = false + } + } + + then { + assert workflow.success : "Pipeline should complete without errors" + assert path("${outputDir}/mutational_profile").exists() + // Negative assertions: disabled features must NOT produce output + assert !path("${outputDir}/mutdensity").exists() + assert !path("${outputDir}/selection/omega").exists() + // Snapshot of stable output + assert snapshot(path("${outputDir}/mutational_profile/all_samples.all.profile.tsv")).match() + } + } +} +``` + +**Active tests (5):** + +| Test | Tag | Purpose | +|------|-----|---------| +| TEST 1 | `input_maf` | MAF-based run, minimal features, snapshot of profile TSV | +| TEST 1b | `input_vcf_with_depths` | VCF + custom depths run | +| TEST 2 | `omega` | Omega analysis enabled; schema + structural + snapshot assertions | +| TEST 3 | `input_maf_nodepths_validation` | Pipeline **fails** when `--input_maf` without `--use_custom_depths` | +| TEST 4 | `input_csvnobam_nodepths_validation` | Fails: no BAMs, no depths | +| TEST 5 | `input_csv_no_bam_no_depthsfile_validation` | Fails: depths enabled but no file | + +**Assertion patterns:** +- Success: `assert workflow.success : "message"` +- Failure tests: `assert workflow.failed : "message"` (validation tests assert the pipeline exits non-zero) +- Directory existence/non-existence: `path("${outputDir}/x").exists()` +- Snapshots: `snapshot(path(...)).match()` — MD5 of file content stored in `tests/deepcsa.nf.test.snap` +- Content schema checks: read file with `path(...).readLines()`, split on `\t`, assert header columns (`gene`, `sample`, `dnds`, `pvalue_adj`), row column-count consistency, and expected sample membership (TEST 2) +- Row-count assertions for non-deterministic float output: `assert dataLines.size() == 252` + +**Handling non-deterministic outputs** (TEST 2 pattern — reuse this for new float-heavy outputs): + +```groovy +// Sort by key columns to avoid floating-point order differences +def sortedRounded = filteredRows.sort { line -> /* key columns */ } + .collect { line -> + def cols = line.split('\t', -1) + (4..7).each { i -> snapshotCols[i] = String.format("%.2f", snapshotCols[i] as Double) } + snapshotCols.join('\t') + } +assert snapshot(sortedContent.join('\n')).match("omega_results") +``` + +**Test config** (`tests/nextflow.config`): +- `executor = 'slurm'`, `errorStrategy = 'retry'`, `maxRetries = 2` +- Singularity with `cacheDir`/`libraryDir` on cluster storage +- `validation.ignoreParams = ['input_maf', 'custom_depths_table']` — skips nf-schema file-existence checks for remote HTTP test inputs +- Large remote test data fetched at runtime from `bbglab/DeepClone_protocol` GitHub repo (no local download step) + +**Snapshot policy** (from `tests/README.md`): +- Regenerate only from the cluster, never locally: `nf-test test tests/deepcsa.nf.test --update-snapshot` (optionally with `--tag `) +- **Mandatory after changing default pipeline parameters** +- Review new hashes in `tests/deepcsa.nf.test.snap` before committing + +**Debugging a failed pipeline test** (from `tests/README.md`): + +```bash +cd /tests//work// +cat .command.out .command.err .command.sh +bash .command.run # reproduce exact environment +cat /tests//meta/nextflow.log +``` + +## Module Tests (nf-test `nextflow_process`) + +**Current state:** exactly **one** module-level test exists — `modules/local/expand_regions/tests/main.nf.test` for the `EXPAND_REGIONS` process (`modules/local/expand_regions/main.nf`). The other ~45 local modules in `modules/local/` have no module tests. + +**Test file convention** (follows the nf-core module test layout): + +``` +modules/local// +├── main.nf +├── meta.yml +└── tests/ + ├── main.nf.test # nextflow_process block + └── main.nf.test.snap # snapshots of process output channels +``` + +**Structure** (`modules/local/expand_regions/tests/main.nf.test`): + +```groovy +nextflow_process { + name "Test EXPAND_REGIONS process" + script "modules/local/expand_regions/main.nf" // path relative to repo root + process "EXPAND_REGIONS" // process name inside main.nf + + test("Testing a run without autoexons and autodomains, it should fail") { + when { + params { + autoexons = false + autodomains = false + subgenic_bedfile = false + } + process { + """ + // Groovy heredoc: set each input channel of the process + input[0] = tuple( + [ id:'test', single_end:false ], + file("${projectDir}/test_data/modules/consensus.exons_splice_sites.PPM1D.tsv") + ) + input[1] = file("${projectDir}/test_data/modules/PPM1D_domains.bed4.bed") + input[2] = file("${projectDir}/test_data/modules/PPM1D_exons.bed4.bed") + input[3] = file("${projectDir}/test_data/dummy_file.tsv") + """ + } + } + + then { + assert !process.success // negative test: process must fail + } + } + + test("Should run with autoexons and autodomains") { + when { + params { autoexons = true; autodomains = true; subgenic_bedfile = false } + process { """ input[0] = ...; input[1] = ... """ } + } + + then { + assert process.success + assert snapshot(process.out).match() // snapshot ALL output channels + } + } +} +``` + +**Key differences from pipeline tests:** + +| Aspect | `nextflow_pipeline` (tests/deepcsa.nf.test) | `nextflow_process` (module tests) | +|---|---|---| +| Scope | whole `main.nf`, all workflows | one process from one module file | +| Inputs | `params {}` (samplesheet, flags) | `process {}` heredoc assigning `input[N]` channels | +| Assertions | `workflow.success/failed`, published dirs | `process.success`, `process.out` channels | +| Snapshots | published output files | `snapshot(process.out).match()` — MD5 per channel element | +| Fixtures | `tests/test_data/` + remote HTTP URLs | repo-root `test_data/` (e.g. `test_data/modules/`) | +| Snapshot file | `tests/deepcsa.nf.test.snap` | `modules/local//tests/main.nf.test.snap` | +| Runtime | SLURM + Singularity, minutes–hours | still uses `tests/nextflow.config` (SLURM), but seconds–minutes | + +**Snapshot format** (`main.nf.test.snap`): JSON keyed by test name → `content` (per-channel MD5 hashes, e.g. `"panel_increased": [["...tsv:md5,6ad5da..."]]`) → `meta` (nf-test/Nextflow versions) → `timestamp`. The `meta` block records tool versions, so snapshots are regenerated when nf-test/Nextflow versions change. + +**⚠️ Discovery caveat:** `nf-test.config` sets `testsDir "tests"` and `ignore 'modules/nf-core/**/*', 'subworkflows/nf-core/**/*'`. This means: +- Running bare `nf-test test` discovers only `tests/deepcsa.nf.test` — module tests under `modules/local/**/tests/` are **not** auto-discovered and must be run by explicit path. +- The `ignore` patterns exclude nf-core modules/subworkflows from testing (they have their own CI upstream), but local modules are *not* ignored — they're just outside `testsDir`. +- Module tests still load `tests/nextflow.config` (via `configFile` in `nf-test.config`), so they submit to SLURM and use Singularity like pipeline tests. They are not runnable on a laptop. + +**Stub-run pattern:** the expand_regions test file contains a commented-out `stub = true` test template (`config { stub = true }` inside `then {}`) — the standard nf-core pattern for testing module stub blocks without executing the real script. Revive it when `main.nf` gains a `stub:` block. + +**Adding a module test (checklist):** +1. Create `modules/local//tests/main.nf.test` with a `nextflow_process` block; `script` path is relative to repo root, `process` is the process name in `main.nf` +2. Add small fixtures under repo-root `test_data/modules/` (keep them tiny — they're committed) +3. In `when { process { """...""" } }`, assign every declared input channel (`input[0]`, `input[1]`, ...) using `file("${projectDir}/...")` +4. Assert `process.success` (or `!process.success` for expected-failure cases) and `snapshot(process.out).match()` +5. Generate the snapshot: `nf-test test modules/local//tests/main.nf.test --update-snapshot` (from the cluster) +6. Commit both `main.nf.test` and `main.nf.test.snap` + +## Python Unit Tests (unittest) + +**Suite Organization** (`bin/test/test_utils_filter.py`): + +```python +#!/usr/bin/env python3 +"""Module docstring listing what is covered.""" + +import sys +import tempfile +import unittest +from pathlib import Path +import pandas as pd + +# Add the bin directory to the path to import sibling modules + +sys.path.insert(0, str(Path(__file__).parent.parent)) +from utils_filter import filter_maf, somatic_mask + +THRESHOLD = 0.3 + +# --------------------------------------------------------------------------- + +# Helpers + +# --------------------------------------------------------------------------- + +def _make_df(rows): + """Build a minimal MAF DataFrame.""" + return pd.DataFrame(...) + +class TestSomaticMask(unittest.TestCase): + """Tests for somatic_mask(maf_df, threshold).""" + + def test_all_below_threshold_is_somatic(self): + """All three VAF columns strictly below threshold → somatic True.""" + ... +``` + +**Patterns:** +- **Path bootstrap (required):** every test file starts with `sys.path.insert(0, str(Path(__file__).parent.parent))` because `bin/` scripts import each other as top-level modules +- **Class per function under test:** `TestSomaticMask`, `TestFilterMaf` — one `TestCase` class per function, named `Test` +- **Docstring = assertion explanation:** each test method's docstring states the expected behavior in plain language +- **setUp/tearDown with temp dirs:** `tempfile.mkdtemp()` in `setUp`, `shutil.rmtree` + `os.chdir` restore in `tearDown` (`bin/test/test_mask_matrix.py`) +- **Fixture-builder helpers:** module-level `_make_df()` / `_make_maf()` functions and instance methods like `create_mock_bed_file(sample_name, positions)` that write small synthetic files +- **Section comment banners** (`# ---- basic behaviour ----`) to group related tests +- **SystemExit handling:** for CLI scripts that call `sys.exit`, tests wrap calls in `try/except SystemExit` and `self.fail(...)` on unexpected exits (`bin/test/test_check_samplesheet.py`) + +**What is tested (113 test methods/classes across 5 files):** +- `check_samplesheet.py` — sample-name validation (security: shell-injection and flag-injection prevention) +- `check_contamination.py` — `compute_shared_variants` on minimal synthetic DataFrames +- `create_mask_matrix.py` + `merge_annotation_depths.apply_mask_matrix` — position × sample mask matrices from synthetic BED files +- `utils_filter.py` — somatic/germline masks, `filter_maf` branches, criteria parsing, BED extraction +- `plot_selectionsideplots.py` — plotting helpers + +## Mocking + +**Framework:** none — `unittest.mock` is **not used** anywhere in `bin/test/`. + +**Patterns:** +- Instead of mocking, tests build **real minimal in-memory fixtures** (small pandas DataFrames) or **real tiny files** in temp dirs (synthetic BED/CSV files) +- The word "mock" appears only in helper names like `create_mock_bed_file` — these create genuine files, not mock objects + +**What to Mock:** nothing currently; follow the existing approach — construct minimal real inputs. + +**What NOT to Mock:** pandas DataFrames and small text files are cheap to create for real; do not introduce `MagicMock` for them. + +## Fixtures and Factories + +**Test Data:** + +```python + +# In-memory DataFrame factory (bin/test/test_utils_filter.py) + +def _make_df(rows: list[tuple[float, float, float]]) -> pd.DataFrame: + vafs, vd_vafs, vaf_ams = zip(*rows) + return pd.DataFrame({"VAF": list(vafs), "vd_VAF": list(vd_vafs), "VAF_AM": list(vaf_ams)}) + +# On-disk fixture writer (bin/test/test_mask_matrix.py) + +def create_mock_bed_file(self, sample_name, positions): + bed_file = f"{sample_name}.flagged-pos.bed" + with open(bed_file, 'w') as f: + for chrom, start, end, filter_val in positions: + f.write(f"{chrom}\t{start}\t{end}\t{filter_val}\n") +``` + +**Location:** +- Python: fixtures are generated inline in test files (no `fixtures/` directory) +- Pipeline: committed small inputs in `tests/test_data/` (`input.csv`, `input_maf.csv`, `input_no_bam.csv`, `test_mutations.maf`); large cohort data fetched at runtime from `https://raw.githubusercontent.com/bbglab/DeepClone_protocol/main/test_datasets/deepCSA/...` + +## Coverage + +**Requirements:** None enforced. No coverage tooling configured (no `pytest.ini`, `tox.ini`, `codecov.yml`, no CI workflows — `.github/` does not exist). + +**View Coverage:** + +```bash +python -m pytest bin/test/ --cov=bin --cov-report=term-missing # if pytest-cov installed +``` + +## Test Types + +**Unit Tests:** +- Python `unittest` in `bin/test/` — pure-function tests on synthetic DataFrames/files, no network, no cluster + +**Integration Tests:** +- nf-test pipeline runs in `tests/deepcsa.nf.test` — full `main.nf` execution on SLURM with Singularity containers, validating published output trees, file schemas, row counts, and MD5 snapshots + +**E2E Tests:** +- The nf-test suite *is* the E2E layer (entire pipeline end-to-end). No browser/UI tests (not applicable). + +## Common Patterns + +**Async Testing:** Not applicable (no async code). + +**Error Testing:** + +```python + +# Python: CLI scripts exit non-zero; tests catch SystemExit (bin/test/test_check_samplesheet.py) + +try: + check_samplesheet(input_file, output_file) + self.fail("Invalid sample name was accepted") # inverted for invalid-input cases +except SystemExit: + pass # expected +``` + +```groovy +// nf-test: assert the whole pipeline fails on invalid params (tests/deepcsa.nf.test, TEST 3) +then { + assert workflow.failed : "Pipeline should fail when --input_maf is set without --use_custom_depths" +} +``` + +**Adding a new pipeline test (checklist):** +1. Add a `test("TEST N. ")` block to `tests/deepcsa.nf.test` with a unique `tag` +2. Set params in `when { params { ... } }`; use `$outputDir` for `outdir` +3. If params point to remote HTTP files, add them to `validation.ignoreParams` in `tests/nextflow.config` +4. Assert success/failure, directory presence/absence, then snapshot only deterministic outputs (round floats, sort rows first) +5. Run from the cluster; update snapshots with `--update-snapshot` if outputs changed intentionally + +**Adding a new Python unit test (checklist):** +1. Create `bin/test/test_.py` with the `sys.path.insert` bootstrap +2. One `TestCase` class per function under test; docstring each test with the expected behavior +3. Build inputs with `_make_*` helper factories or temp-dir files; clean up in `tearDown` +4. Run: `python -m unittest discover -s bin/test -v` + +## Gaps to Be Aware Of + +- **No CI:** `.github/` does not exist — neither test layer runs automatically on push/PR +- **~71 of 80 `bin/` scripts have no unit tests** — only 5 modules are covered +- **~45 of 46 local Nextflow modules have no module tests** — only `expand_regions` has a `nextflow_process` test; module tests are also not auto-discovered by `nf-test.config` (`testsDir "tests"`) +- **No pytest config / coverage gate** — coverage is unmeasured +- **nf-test requires SLURM + Singularity** — cannot run the integration suite on a laptop or in plain CI without a cluster runner +- **Commented-out test templates** at the bottom of `tests/deepcsa.nf.test` (MAF + precomputed depths integration test, trace-based process-count assertions) are ready-made patterns to revive when test assets improve + +--- + +*Testing analysis: 2026-10-09* diff --git a/CHANGELOG.md b/CHANGELOG.md index ee299ee6..17be1688 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,15 @@ Initial release of bbglab/deepCSA, created with the [nf-core](https://nf-co.re/) ### `Added` +- New output `panel_exons_protein_intervals.tsv` in the `DNA2PROTEINMAPPING` step (published to `regions/annotations/`), reporting for each exon of the panel transcripts the interval of protein coordinates it covers, using the same structure as the domains info file (`Ens_Transcr_ID`, `Begin`, `End`, `NAME`, `GENE`, `DOMAIN_ID`). +- Saturation kinetics curves (`COMPUTE_SATURATION_KINETICS` step, published to `plots/saturation_kinetics/`). For each group and each (resolution, impact) combination — genomic/residue × protein_affecting/nonsense/truncating/missense/synonymous — the step produces: + - `{group}.curves/{sites}_{impact}_empirical.pdf` — empirical discovery index curves (proportion of mutated sites vs sequencing depth) for all genes, obtained by downsampling the observed mutations with Bernoulli replicates. + - `{group}.curves/{sites}_{impact}_theoretical_empirical.pdf` — the same empirical curves overlaid with the theoretical neutral saturation curve derived from the per-site relative mutability and the synonymous mutation rate. + - `{group}.curves/{sites}_{impact}_slopes.pdf` — per-gene comparison of the rate of change (Δ proportion / Δ log10 depth) of the empirical curve against the theoretical neutral curve over identical depth intervals. + - `{group}_mutations_{sites}_rates.{impact}.tsv` — per-gene/per-site unique-mutation probabilities at each subsampling depth. + - `{group}_slopes_{sites}.{impact}.tsv` — per-gene interval slopes (empirical, theoretical and their ratio) together with the depth bounds and the proportion of covered positions at the middle of each interval, for cross-run comparison. + - Requires both `--omega` and one of the mutability-driven analyses (`--oncodrivefml`, `--oncodriveclustl` or `--oncodrive3d`) to be enabled, since it consumes the omega preprocessing mutability table and the relative mutability per site. + ### `Fixed` ### `Dependencies` diff --git a/assets/assess_panel/assess_panel.py b/assets/assess_panel/assess_panel.py index eb6aae5e..620f7fc9 100644 --- a/assets/assess_panel/assess_panel.py +++ b/assets/assess_panel/assess_panel.py @@ -9,7 +9,7 @@ import click -sys.path.append("../../bin") +sys.path.append("../../bin/") from utils_context import triplet_context_iterator diff --git a/bin/add_subgenicregions.py b/bin/add_subgenicregions.py index 655014f1..eccaead2 100755 --- a/bin/add_subgenicregions.py +++ b/bin/add_subgenicregions.py @@ -68,7 +68,7 @@ def main(panel_file, autoexons, autodomains, custom, subgenic_regions_complement # If you reach ths point, ind_start is defined and valid # Search for END position starting from ind_start - search_end = chr_data.iloc[ind_start:,:] + search_end = chr_data.iloc[ind_start:] end_matches = np.where(search_end["POS"] == row["END"])[0] end_found = len(end_matches) > 0 @@ -117,7 +117,7 @@ def main(panel_file, autoexons, autodomains, custom, subgenic_regions_complement upd_end = ind_end # Extract subgenic data and modify gene names - subgenic_data = chr_data.iloc[upd_start: upd_end + 1, :].copy() + subgenic_data = chr_data.iloc[upd_start: upd_end + 1].copy() subgenic_data["GENE"] = region_name new_data = pd.concat((new_data, subgenic_data)) diff --git a/bin/check_samplesheet.py b/bin/check_samplesheet.py index ada3e578..8af3eb05 100755 --- a/bin/check_samplesheet.py +++ b/bin/check_samplesheet.py @@ -211,7 +211,8 @@ def check_samplesheet(file_in, file_out, bam_required=False): with file_in.open(newline="") as in_handle: reader = csv.DictReader(in_handle, dialect=sniff_format(in_handle)) # Validate the existence of the expected header columns. - if not required_columns.issubset(reader.fieldnames): + fieldnames = reader.fieldnames or [] + if not required_columns.issubset(fieldnames): req_cols = ", ".join(required_columns) logger.critical(f"The sample sheet **must** contain these column headers: {req_cols}.") sys.exit(1) @@ -224,7 +225,7 @@ def check_samplesheet(file_in, file_out, bam_required=False): logger.critical(f"{str(error)} On line {i + 2}.") sys.exit(1) checker.validate_unique_samples() - header = list(reader.fieldnames) + header = list(fieldnames) # See https://docs.python.org/3.9/library/csv.html#id3 to read up on `newline=""`. with file_out.open(mode="w", newline="") as out_handle: writer = csv.DictWriter(out_handle, header, delimiter=",") diff --git a/bin/create_mask_matrix.py b/bin/create_mask_matrix.py index ba66a870..5a91592b 100755 --- a/bin/create_mask_matrix.py +++ b/bin/create_mask_matrix.py @@ -63,7 +63,7 @@ def add_bed_positions(bed_df: pd.DataFrame, entries_added = 0 for row in bed_df.itertuples(): - for pos in range(row.START, row.END + 1): + for pos in range(int(row.START), int(row.END) + 1): key = (row.CHROM, pos, sample_name) if key not in masked_positions: mask_data.append({ diff --git a/bin/mutations_custom_processing.py b/bin/mutations_custom_processing.py index 24a8643a..108ad36b 100755 --- a/bin/mutations_custom_processing.py +++ b/bin/mutations_custom_processing.py @@ -14,8 +14,32 @@ def customize_annotations(mutation_summary_file, custom_regions_file, simple = True ): """ - # TODO explain what this function does - + Override the gene and consequence annotations of the mutations falling in + custom regions. + + The custom regions file provides, per pyrimidine mutation ID (MUT_ID_pyr), + a replacement gene and consequence. Mutations matching a custom region get + their SYMBOL and Consequence (and the canonical_ counterparts) replaced by + the custom values; the VEP-derived positional annotations (Feature, + Protein_position, Amino_acids, ...) are blanked since they no longer apply, + and the derived consequence columns (Consequence_single, + Consequence_broader, Protein_affecting and their canonical_ versions) are + recomputed from the new consequence. Mutations not present in the custom + regions file keep their original annotation untouched. If the custom + regions file is empty or no mutation matches it, the input is copied to + the output unchanged. + + Parameters + ---------- + mutation_summary_file : str + Path to the annotated mutation summary table (TSV). + custom_regions_file : str + Path to the custom regions table (TSV) with at least the columns + MUT_ID_pyr, GENE and IMPACT. + customized_mutations_output : str + Path where the customized mutation summary is written (TSV). + simple : bool, optional + Unused; kept for backwards compatibility. """ # simple = ['CHROM', 'POS', 'REF', 'ALT', 'MUT_ID' , 'GENE', 'IMPACT' , 'CONTEXT_MUT', 'CONTEXT'] # rich = ['CHROM', 'POS', 'REF', 'ALT', 'MUT_ID', 'STRAND', 'GENE', 'IMPACT', 'Feature', 'Protein_position', 'Amino_acids', 'CONTEXT_MUT', 'CONTEXT'] diff --git a/bin/panels_computedna2protein.py b/bin/panels_computedna2protein.py index 55e01af0..09859ece 100755 --- a/bin/panels_computedna2protein.py +++ b/bin/panels_computedna2protein.py @@ -544,6 +544,45 @@ def get_dna2prot_depth(gene_n_transcript_info: pd.DataFrame, depth_file: str, co return dna_prot_df, exons_coord_df_final +def get_exons_protein_intervals(exons_depth: pd.DataFrame) -> pd.DataFrame: + """ + Derive the interval of protein coordinates covered by each exon of the panel transcripts, + using the same structure as the domains info file: + 'Ens_Transcr_ID', 'Begin', 'End', 'NAME', 'GENE', 'DOMAIN_ID' + + Parameters + ------------ + exons_depth : pandas.DataFrame + DataFrame containing the per-position DNA to protein mapping with depth information, + as produced by `get_dna2prot_depth`. It must contain the 'GENE', 'PROT_POS', + 'EXON_ID' and 'TRANSCRIPT_ID' columns. + + Returns + ------------ + pandas.DataFrame + A DataFrame with one row per exon containing coding positions, with the transcript ID, + the interval of protein coordinates it covers ('Begin', 'End', 1-based inclusive), + the gene-prefixed exon name ('NAME', as used in the panel files), the gene, + and the bare exon identifier ('DOMAIN_ID'). + """ + cds_positions = exons_depth[exons_depth["PROT_POS"].notna() & exons_depth["EXON_ID"].notna()] + exons_protein = cds_positions.groupby("EXON_ID").agg(Ens_Transcr_ID=("TRANSCRIPT_ID", "first"), + GENE=("GENE", "first"), + Begin=("PROT_POS", "min"), + End=("PROT_POS", "max") + ).reset_index() + exons_protein[["Begin", "End"]] = exons_protein[["Begin", "End"]].astype(int) + + # Mirror the domains file structure: 'NAME' is the gene-prefixed exon id + # (as used in the panel files) and 'DOMAIN_ID' is the bare exon identifier. + exons_protein["NAME"] = exons_protein["EXON_ID"] + exons_protein["DOMAIN_ID"] = exons_protein["EXON_ID"].str.split("--").str[-1] + exons_protein = exons_protein[["Ens_Transcr_ID", "Begin", "End", "NAME", "GENE", "DOMAIN_ID"]] + exons_protein = exons_protein.sort_values(by = ["GENE", "Begin", "NAME"]).reset_index(drop = True) + + return exons_protein + + # Plots # ---------------------------------------------------------- def plot_coverage_per_gene(depths_df: pd.DataFrame) -> None: @@ -673,6 +712,9 @@ def main(mutations_file, consensus_file, depths_file, ensembl_species, ensembl_g exons_coordinates_bed_like = exons_coord_id[['Chr', 'Start', 'End', 'ID']] exons_coordinates_bed_like.to_csv("panel_exons.bed4.bed", header = False, index = False, sep = '\t') + exons_protein_intervals = get_exons_protein_intervals(exons_depth) + exons_protein_intervals.to_csv("panel_exons_protein_intervals.tsv", header = True, index = False, sep = '\t') + plot_coverage_per_gene(exons_depth) LOG.info("All done!") diff --git a/bin/saturation_kinetics_curves.py b/bin/saturation_kinetics_curves.py new file mode 100755 index 00000000..81dc2320 --- /dev/null +++ b/bin/saturation_kinetics_curves.py @@ -0,0 +1,520 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + + +import os + +import click + +import numpy as np +import pandas as pd + +import matplotlib +matplotlib.use("Agg") # headless execution inside containers +import matplotlib.pyplot as plt +from matplotlib.backends.backend_pdf import PdfPages + +import tensorflow as tf +import tensorflow_probability as tfp +tfd = tfp.distributions + + +CONSEQUENCES_CATEGORIES = { + 'protein_affecting' : { 'nonsense', 'missense', 'essential_splice', 'protein_altering_variant'}, + 'nonsense': {'nonsense'}, + 'truncating': {'nonsense', 'essential_splice'}, + 'missense' : { 'missense' }, + 'synonymous' : { 'synonymous' } +} + +def prob_min_uniform_sample_below_cut_vec(N, n, cut): + """ + Probability that the minimum of a sampling of n elements from [N] is lower or equal than cut + Vectorized version for speed up. + Computed for arrays of (N, n, cut) triples. + """ + N = np.asarray(N, dtype=int) + n = np.asarray(n, dtype=int) + cut = np.asarray(cut, dtype=int) + + out = np.zeros_like(N, dtype=float) + out[cut >= N] = 1.0 + valid = (N > 0) & (n > 0) & (cut > 0) & (cut < N) + + for j in np.where(valid)[0]: + k = np.arange(int(n[j])) + log_terms = np.log(N[j] - cut[j] - k) - np.log(N[j] - k) + out[j] = 1 - np.exp(log_terms.sum()) + + return out + +def load_mutations(somatic_mutations_file): + somatic_mutations = pd.read_csv(somatic_mutations_file, sep='\t', low_memory=False) + mutations = somatic_mutations[ + ~(somatic_mutations['FILTER'].str.contains("not_in_exons")) + & (somatic_mutations['TYPE'] == 'SNV') + ] + + mutations_lite = mutations[['CHROM', 'POS', 'REF', 'ALT', 'SAMPLE_ID', 'ALT_DEPTH', 'ALT_DEPTH_AM']] + mutations_lite = mutations_lite.groupby(by=['CHROM', 'POS', 'REF', 'ALT'] + ).agg({'ALT_DEPTH_AM': 'sum', 'ALT_DEPTH': 'sum'}).reset_index() + return mutations_lite + + +def get_aachange_format(r): + + if r['Protein_position'] == '-': + return '-' + elif len(r['Amino_acids']) == 1: + return r['Amino_acids'] + r['Protein_position'] + r['Amino_acids'] + else: + return r['Amino_acids'].split('/')[0] + r['Protein_position'] + r['Amino_acids'].split('/')[1] + + +def collect_vep(vep_filename): + vep_panel = pd.read_csv(vep_filename, sep='\t') + dg = vep_panel[['CHROM', 'POS', 'REF', 'ALT', 'Protein_position', 'Amino_acids', 'GENE']].copy() + dg['POS'] = dg['POS'].astype(int) + dg['AACHANGE'] = dg.apply(get_aachange_format, axis=1) + return dg + + +def load_panel(consensus_panel_file, depths_file, samples_group, vep_annotations): + + df_panel = pd.read_csv(consensus_panel_file, sep='\t') + df_depth = pd.read_csv(depths_file, sep='\t') + + df_panel = pd.merge(df_panel, df_depth[['CHROM', 'POS', samples_group]], on=['CHROM', 'POS'], how='left') + df_panel.rename(columns={samples_group: 'DEPTH'}, inplace=True) + df_panel = pd.merge(df_panel, vep_annotations[['CHROM', 'POS', 'REF', 'ALT', 'AACHANGE', 'GENE']], + left_on=['CHROM', 'POS', 'REF', 'ALT', 'GENE'], + right_on=['CHROM', 'POS', 'REF', 'ALT', 'GENE'], + how='left') + + df_panel = df_panel[df_panel['AACHANGE'] != '-'] + df_panel.dropna(inplace=True) + df_panel['RESIDUE'] = df_panel['AACHANGE'].apply(lambda s: s[:-1]) + return df_panel + + + + +def empirical_discovery_index_curve(gene, mutations_dict, df_panel_dict, + subsampling_rates, + replicates=100, sites='genomic'): + + df = mutations_dict[sites] + df = df[df['GENE'] == gene] + + dg = df_panel_dict[sites] + dg = dg[dg['GENE'] == gene] + + size = dg.shape[0] + mean_depth = df['DEPTH'].mean() + + x, mean, err_low, err_high = [], [], [], [] + + for i, p in enumerate(subsampling_rates): + dist_bernoulli = tfd.Bernoulli(probs=df[f'UNIQUE_RATE_{i}'].values) + unique_mutations = np.sum(dist_bernoulli.sample(sample_shape=(replicates,)), axis=1) + y = list(unique_mutations / size) + mean += [np.mean(y)] + err_low += [np.percentile(y, 2.5)] + err_high += [np.percentile(y, 97.5)] + x += [mean_depth * p] + mean += [df.shape[0] / size] + err_low += [df.shape[0] / size] + err_high += [df.shape[0] / size] + x += [mean_depth] + + return x, mean, err_low, err_high + + +def compute_interval_slopes(x, y): + """ + Rate of change of y across each consecutive interval of x, in log-space: + + slope_i = (y[i+1] - y[i]) / (log10(x[i+1]) - log10(x[i])) + + Returns (midpoints, slopes), where midpoints is the geometric mean of the + two depth bounds of each interval. Intervals with non-increasing or + non-finite depths get a NaN slope. + """ + x = np.asarray(x, dtype=float) + y = np.asarray(y, dtype=float) + d_log = np.diff(np.log10(x)) + d_y = np.diff(y) + slopes = np.full(d_y.shape, np.nan) + valid = np.isfinite(d_log) & (d_log > 0) & np.isfinite(d_y) + slopes[valid] = d_y[valid] / d_log[valid] + midpoints = interval_midpoints(x) + return midpoints, slopes + + +def interval_midpoints(x): + """Geometric mean of the two depth bounds of each consecutive interval.""" + x = np.asarray(x, dtype=float) + return np.sqrt(x[:-1] * x[1:]) + + +def interval_mid_values(y): + """Average of the two endpoint values of each consecutive interval.""" + y = np.asarray(y, dtype=float) + return 0.5 * (y[:-1] + y[1:]) + + +def plot_empirical_discovery(gene, mutations_dict, df_panel_dict, subsampling_rates, sites='genomic', + impact="protein_affecting", + pdf=None): + click.echo("Plotting empirical discovery") + x, mean, err_low, err_high = empirical_discovery_index_curve(gene, mutations_dict, df_panel_dict, subsampling_rates, sites = sites) + fig, ax1 = plt.subplots(figsize=(2,2)) + + ax1.scatter(x, mean, s=50) + for i, m in enumerate(x): + ax1.vlines(m, err_low[i], err_high[i]) + + ax1.set_xscale('log') + if sites == 'residue': + ax1.set_ylabel('proportion of\nmutated residues') + elif sites == 'genomic': + ax1.set_ylabel('proportion of\nmutated nucleotides') + + ax1.set_xlabel('depth per residue') + ax1.spines['top'].set_visible(False) + ax1.spines['right'].set_visible(False) + ax1.set_ylim(0, max(err_high) * 1.1) + plt.title(f"{gene} ({impact}, {sites})") + if pdf is not None: + pdf.savefig(fig, bbox_inches='tight', dpi=300) + plt.close(fig) + + +def plot_slope_comparison(gene, x_empirical, y_empirical, y_theoretical, + slopes_empirical, slopes_theoretical, + sites='genomic', impact="protein_affecting", pdf=None): + """ + Compare the rate of change (log-space slope) of the empirical curve against + the theoretical neutral curve, computed over the same intervals. + + Two pages are added to the PDF: the slope against depth (log-scale) and the + slope against the proportion of positions covered at the middle of each + interval. + """ + click.echo(f"Plotting slope comparison for {gene}") + + midpoints = interval_midpoints(x_empirical) + proportions_empirical = interval_mid_values(y_empirical) + proportions_theoretical = interval_mid_values(y_theoretical) + + if sites == 'residue': + proportion_label = 'proportion of\nmutated residues' + else: + proportion_label = 'proportion of\nmutated nucleotides' + + panels = [ + ('depth per residue', midpoints, midpoints, True), + (proportion_label, proportions_empirical, proportions_theoretical, False), + ] + + for xlabel, x_emp, x_theo, logx in panels: + fig, ax1 = plt.subplots(figsize=(3, 2.5)) + + ax1.scatter(x_emp, slopes_empirical, color='brown', s=25, label='empirical') + ax1.plot(x_theo, slopes_theoretical, color='grey', lw=2, alpha=0.5, label='neutral theoretical') + + if logx: + ax1.set_xscale('log') + ax1.set_xlabel(xlabel) + ax1.set_ylabel('rate of change\n(Δ proportion / Δ log10 depth)') + ax1.spines['top'].set_visible(False) + ax1.spines['right'].set_visible(False) + ax1.legend(loc='best', fontsize=6) + + plt.title(f"{gene} ({impact}, {sites})") + if pdf is not None: + pdf.savefig(fig, bbox_inches='tight', dpi=300) + plt.close(fig) + + + +# theoretical neutral vs empirical discovery curves +# Create a PDF to save the plots +def main_empirical(sample, + mutations_dict, panel_df, + df_panel_dict, + omega_mutability_file, relative_mutability_file, + subsampling_rates, + sites='genomic', impact = "protein_affecting", logscale=False, genes_list = None, + empirical_pdf=None, combined_pdf=None, slope_pdf=None): + + # retrieve relative mutability + mutability_raw = pd.read_csv(relative_mutability_file, sep='\t', + header=None, names=['CHROM', 'POS', 'REF', 'ALT', 'MUTABILITY']) + + mutability_raw = pd.merge(mutability_raw, panel_df, on=['CHROM', 'POS', 'REF', 'ALT'], how='left') + mutability_raw = mutability_raw.rename(columns={sample: "MUTABILITY"}) + + if genes_list is None: + genes_list = sorted(panel_df['GENE'].unique()) + + mutations_lite = mutations_dict[sites] + mutations_lite = mutations_lite.assign( + VAF=mutations_lite['ALT_DEPTH'] / mutations_lite['DEPTH']) + + all_slope_records = [] + + for gene in genes_list: + try: + plot_empirical_discovery(gene, mutations_dict, df_panel_dict, subsampling_rates, sites=sites, + impact=impact, pdf=empirical_pdf) + + synonymous_mutation_rate = pd.read_csv(omega_mutability_file, sep='\t') + synonymous_mutation_rate = synonymous_mutation_rate[synonymous_mutation_rate['GENE'] == gene] + mutability_gene = mutability_raw[mutability_raw['GENE'] == gene] + mutability_gene = pd.merge(mutability_gene, synonymous_mutation_rate[['CONTEXT_MUT', sample]], on=['CONTEXT_MUT'], how='left') + mutability_gene.rename(columns={sample: 'MUTRATE'}, inplace=True) + + # discard positions in non-CDS regions, probably splicing and intronic + mutability_gene = mutability_gene[(mutability_gene['AACHANGE'] != '-') & (~mutability_gene['AACHANGE'].isnull())] + + # keep only sites of interest + mutability_gene = mutability_gene[mutability_gene["IMPACT"].isin(CONSEQUENCES_CATEGORIES[impact])] + + mutability_gene['RESIDUE'] = mutability_gene['AACHANGE'].apply(lambda s: s[:-1]) + + if sites == 'residue': + mutability_gene = mutability_gene.groupby(['GENE', 'RESIDUE']).agg({'MUTRATE': 'sum', 'DEPTH': 'mean'}).reset_index() + elif sites == 'genomic': + mutability_gene = mutability_gene.groupby(['GENE', 'POS']).agg({'MUTRATE': 'sum', 'DEPTH': 'mean'}).reset_index() + + mutability_gene['MUTABILITY'] = mutability_gene.apply(lambda s: (s['MUTRATE'] / s['DEPTH']), axis=1) + + if sites == 'residue': + mutability_gene = pd.merge(mutability_gene, mutations_lite[['GENE', 'RESIDUE', 'VAF']], on=['GENE', 'RESIDUE'], how='left') + elif sites == 'genomic': + mutability_gene = pd.merge(mutability_gene, mutations_lite[['GENE', 'POS', 'VAF']], on=['GENE', 'POS'], how='left') + + + # neutral rate + mutability_gene['RATE_NEUTRAL'] = mutability_gene['MUTABILITY'] + + # compute saturation theoretical + y_unique_neutral = [] + if sites == 'genomic': + x_theoretical = np.logspace(3, 8, num=100) + elif sites == 'residue': + x_theoretical = np.logspace(3, 7, num=100) + + for depth in x_theoretical: + unique_neutral = np.sum(1 - np.exp(-mutability_gene['RATE_NEUTRAL'].to_numpy(dtype=float) * depth)) / mutability_gene.shape[0] + y_unique_neutral.append(unique_neutral) + + # compute empirical discovery index curve + x_empirical, mean, err_low, err_high = empirical_discovery_index_curve(gene, mutations_dict, df_panel_dict, + subsampling_rates, sites=sites) + + # rate of change per interval: empirical vs theoretical + # (theoretical curve interpolated at the empirical depths so both + # slopes are computed over identical intervals) + y_theoretical_at_empirical = np.interp(np.log10(x_empirical), np.log10(x_theoretical), y_unique_neutral) + midpoints, slopes_empirical = compute_interval_slopes(x_empirical, mean) + _, slopes_theoretical = compute_interval_slopes(x_empirical, y_theoretical_at_empirical) + mid_proportions = interval_mid_values(mean) + + slope_records = [ + { + 'GENE': gene, 'SITES': sites, 'IMPACT': impact, + 'DEPTH_LOW': x_empirical[i], 'DEPTH_HIGH': x_empirical[i + 1], + 'DEPTH_MID': midpoints[i], + 'PROPORTION_COVERED': mid_proportions[i], + 'SLOPE_EMPIRICAL': slopes_empirical[i], + 'SLOPE_THEORETICAL': slopes_theoretical[i], + 'SLOPE_RATIO': (slopes_empirical[i] / slopes_theoretical[i] + if np.isfinite(slopes_theoretical[i]) and slopes_theoretical[i] != 0 + else np.nan), + } + for i in range(len(x_empirical) - 1) + ] + all_slope_records.extend(slope_records) + + # plot + fig, ax1 = plt.subplots(figsize=(2,2)) + ax1.set_xscale('log') + if logscale: + ax1.set_yscale('log') + + # empirical + ax1.scatter(x_empirical[:-1], mean[:-1], label='downsampling', color='brown', s=5) + ax1.scatter(x_empirical[-1], mean[-1], label='observed', color='white', edgecolors='brown', alpha=1, s=100) + for i, m in enumerate(x_empirical[:-2]): + plt.vlines(m, err_low[i], err_high[i], color='brown', lw=1) + + # theoretical + + ax1.plot(x_theoretical, y_unique_neutral, color='grey', lw=2, label='neutral theoretical', alpha=0.5) # neutral + if sites == 'residue': + ax1.set_ylabel('proportion of\nmutated residues') + elif sites == 'genomic': + ax1.set_ylabel('proportion of\nmutated nucleotides') + + ax1.set_xlabel('depth per residue') + + ax1.spines['top'].set_visible(False) + ax1.spines['right'].set_visible(False) + + ax1.set_ylim(0,1) + ax1.set_xlim(x_theoretical[0], x_theoretical[-1]) + + plt.title(f"{gene} ({impact}, {sites})") + + if combined_pdf is not None: + combined_pdf.savefig(fig, bbox_inches='tight', dpi=300) + + plt.close(fig) + + # slope comparison plot + plot_slope_comparison(gene, x_empirical, mean, y_theoretical_at_empirical, + slopes_empirical, slopes_theoretical, + sites=sites, impact=impact, pdf=slope_pdf) + + except Exception as e: + click.echo(f"Error occurred while processing {gene}: {e}") + continue + + # save the slope table for all genes + if all_slope_records: + df_slopes = pd.DataFrame(all_slope_records) + slope_table_file = f"{sample}_slopes_{sites}.{impact}.tsv" + df_slopes.to_csv(slope_table_file, sep="\t", index=False) + click.echo(f"Slope table saved to {slope_table_file}") + return df_slopes + + return None + +def compute_mutation_rates(mutations, name, impact, subsampling_rates_list, residue=False): + """Compute mutation rates for a given set of mutations.""" + + grp_by = ['GENE', 'RESIDUE'] if residue else ['GENE', 'POS'] + mutations_lite = mutations.groupby( + grp_by + ).agg({'ALT_DEPTH': 'sum', 'DEPTH': 'mean'}).reset_index() + + depth = mutations_lite["DEPTH"].to_numpy() + alt_depth = mutations_lite["ALT_DEPTH"].to_numpy() + + click.echo(f"Computing mutation rates for {name}") + + for i, p in enumerate(subsampling_rates_list): + n = np.floor(p * depth).astype(int) + mutations_lite[f"UNIQUE_RATE_{i}"] = prob_min_uniform_sample_below_cut_vec(depth, n, alt_depth) + + out_file = f"{name}_mutations_{'residue' if residue else 'genomic'}_rates.{impact}.tsv" + + click.echo(f"Saving mutation rates to {out_file}") + mutations_lite.to_csv(out_file, sep="\t", index=False) + + return mutations_lite + + + + +@click.command() +@click.option("--somatic-mutations-file", required=True, type=click.Path(exists=True),help="Path to the somatic mutations file") +@click.option("--vep-file", required=True, type=click.Path(exists=True),help="Path to the VEP annotation file") +@click.option("--consensus-panel-file", required=True, type=click.Path(exists=True),help="Path to the exons consensus panel file") +@click.option("--depths-file", required=True, type=click.Path(exists=True),help="Path to the depths file") +@click.option("--omega-mutability-file", required=True, type=click.Path(exists=True),help="Path to the omega mutability file") +@click.option("--relative-mutability-file", required=True, type=click.Path(exists=True),help="Path to the relative mutability file") +@click.option("--resolution", type=click.Choice(['genomic', 'residue', 'genomic,residue']), show_default=True, default='genomic,residue', help="either genomic or residue based sites") +@click.option("--group-name", type=str, default="all_samples", show_default=True, help="Name of the group/sample to be used in the code") +def cli(somatic_mutations_file, vep_file, consensus_panel_file, + omega_mutability_file, relative_mutability_file, + depths_file, resolution, group_name, + ): + subsampling_rates = np.logspace(-2, np.log10(0.9), num=20) + + click.echo(f"Analyzing {group_name}") + + mutations = load_mutations(somatic_mutations_file) + click.echo("Mutations loaded") + + vep = collect_vep(vep_file) + click.echo("VEP data collected") + + df_panel_orig = load_panel(consensus_panel_file, depths_file, group_name, vep) + click.echo(f"Panel loaded with {df_panel_orig.shape[0]} sites") + + curves_folder = f'{group_name}.curves' + os.makedirs(curves_folder, exist_ok=True) + + # df_panel represents the total number of mutable sites, + # either genomic or residue sites + for impact in CONSEQUENCES_CATEGORIES.keys(): + df_panel = df_panel_orig[df_panel_orig["IMPACT"].isin(CONSEQUENCES_CATEGORIES[impact])] + + click.echo(f"Panel filtered with {df_panel.shape[0]} sites for {impact} mutations") + mutations_lite = pd.merge(mutations, + df_panel[['CHROM', 'POS', 'REF', 'ALT', 'DEPTH', 'GENE', 'AACHANGE', 'RESIDUE']], + on=['CHROM', 'POS', 'REF', 'ALT'], + how='left') + mutations_lite.dropna(axis=0, inplace=True) # drop nan sites + + if 'residue' in resolution: + mutations_residue = compute_mutation_rates(mutations_lite, group_name, impact, subsampling_rates, residue=True) + if 'genomic' in resolution: + mutations_genomic = compute_mutation_rates(mutations_lite, group_name, impact, subsampling_rates, residue=False) + + click.echo("Mutation rates computed") + + # load mutations + mutations_dict = { + 'genomic': mutations_genomic if 'genomic' in resolution else None, + 'residue': mutations_residue if 'residue' in resolution else None + } + + df_panel_genomic = df_panel.groupby(['POS', 'GENE']).agg({'DEPTH': 'mean'}).reset_index() + df_panel_residue = df_panel.groupby(['RESIDUE', 'GENE']).agg({'DEPTH': 'mean'}).reset_index() + df_panel_dict = { + 'genomic': df_panel_genomic, + 'residue': df_panel_residue + } + + + if 'residue' in resolution: + click.echo("Plotting empirical discovery for residue sites") + with PdfPages(f'{curves_folder}/residue_{impact}_empirical.pdf') as empirical_pdf, \ + PdfPages(f'{curves_folder}/residue_{impact}_theoretical_empirical.pdf') as combined_pdf, \ + PdfPages(f'{curves_folder}/residue_{impact}_slopes.pdf') as slope_pdf: + main_empirical(group_name, mutations_dict, df_panel, df_panel_dict, + omega_mutability_file, relative_mutability_file, + subsampling_rates, + sites='residue', + impact = impact, + logscale=False, + empirical_pdf=empirical_pdf, + combined_pdf=combined_pdf, + slope_pdf=slope_pdf, + ) + + if 'genomic' in resolution: + click.echo("Plotting empirical discovery for genomic sites") + with PdfPages(f'{curves_folder}/genomic_{impact}_empirical.pdf') as empirical_pdf, \ + PdfPages(f'{curves_folder}/genomic_{impact}_theoretical_empirical.pdf') as combined_pdf, \ + PdfPages(f'{curves_folder}/genomic_{impact}_slopes.pdf') as slope_pdf: + main_empirical(group_name, mutations_dict, df_panel, df_panel_dict, + omega_mutability_file, relative_mutability_file, + subsampling_rates, + sites='genomic', + impact = impact, + logscale=False, + empirical_pdf=empirical_pdf, + combined_pdf=combined_pdf, + slope_pdf=slope_pdf, + ) + + # no per-gene files are written anymore; PDFs are closed by their context managers + + +if __name__ == '__main__': + cli() \ No newline at end of file diff --git a/bin/summarize_germline_snps.py b/bin/summarize_germline_snps.py new file mode 100644 index 00000000..2ee273b2 --- /dev/null +++ b/bin/summarize_germline_snps.py @@ -0,0 +1,452 @@ +#!/usr/bin/env python + +""" +Summarize germline SNPs from the clean and somatic mutation tables. + +The germline mutations of each sample are obtained by subtracting the somatic +calls from the clean calls (an anti-join on SAMPLE_ID and MUT_ID). Each germline +variant is then assigned to a VAF-based genotype bin, and a PCA is run on the +samples x variants genotype matrix restricted to the heterozygous-like bins. + +Outputs (written to the current working directory): + {output_prefix}.germline.mutations.tsv Germline mutations with their genotype bin and pathogenic flag. + {output_prefix}.pathogenic_snps.tsv Candidate pathogenic germline SNPs with their annotation. + {output_prefix}.pathogenic_snps_summary.tsv Number of candidate pathogenic germline SNPs per sample. + {output_prefix}.ancestry_inference.tsv Mean gnomAD population AF profile per sample and most likely ethnic group. + {output_prefix}.germline_snps_summary.pdf VAF bins, pathogenic SNPs per sample, ancestry heatmap, PCA elbow and PC1 vs PC2 scatter. +""" + +import click +import matplotlib + +matplotlib.use("Agg") + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import seaborn as sns +from matplotlib.backends.backend_pdf import PdfPages +from sklearn.decomposition import PCA +from sklearn.preprocessing import StandardScaler + +from utils_filter import germline_mask + +# VAF cuts separating the genotype bins: (0, 0.25], (0.25, 0.75], (0.75, 1] +GENOTYPE_BINS = [0, 0.25, 0.75, 1] +# Genotype bins kept for the PCA (heterozygous-like VAF ranges) +GENOTYPE_BINS_FOR_PCA = [1, 2] +NUMBER_OF_COMPONENTS = 6 +# Minimum VAF for a germline variant to be considered a SNP at all +MIN_SNP_VAF = 0.25 +# gnomAD allele frequency columns used to assess rarity; all the available ones +# must be below the threshold for a variant to be considered rare +GNOMAD_AF_COLUMNS = ["gnomADg_AF", "gnomADe_AF"] +# Population-specific gnomAD allele frequency columns used to infer the most +# likely ethnic group of each sample; the genome columns are preferred and the +# exome ones are used as fallback +GNOMAD_GENOME_POP_COLUMNS = { + "AFR": "gnomADg_AFR_AF", + "AMI": "gnomADg_AMI_AF", + "AMR": "gnomADg_AMR_AF", + "ASJ": "gnomADg_ASJ_AF", + "EAS": "gnomADg_EAS_AF", + "FIN": "gnomADg_FIN_AF", + "MID": "gnomADg_MID_AF", + "NFE": "gnomADg_NFE_AF", + "OTH": "gnomADg_OTH_AF", + "SAS": "gnomADg_SAS_AF", +} +GNOMAD_EXOME_POP_COLUMNS = { + "AFR": "gnomADe_AFR_AF", + "AMR": "gnomADe_AMR_AF", + "ASJ": "gnomADe_ASJ_AF", + "EAS": "gnomADe_EAS_AF", + "FIN": "gnomADe_FIN_AF", + "NFE": "gnomADe_NFE_AF", + "OTH": "gnomADe_OTH_AF", + "SAS": "gnomADe_SAS_AF", +} +# Value of the canonical_Protein_affecting column for protein-affecting variants +# (nonsense, missense and essential splice consequences) +PROTEIN_AFFECTING_VALUE = "protein_affecting" +# Columns reported in the per-sample pathogenic SNPs table +PATHOGENIC_SNP_COLUMNS = [ + "SAMPLE_ID", + "MUT_ID", + "canonical_SYMBOL", + "canonical_Consequence_broader", + "canonical_Amino_acids", + "VAF", + "gnomADg_AF", + "gnomADe_AF", +] + + +def subtract_somatic_mutations(clean_mutations, somatic_mutations): + """ + Remove the somatic mutations from the clean mutations to keep the germline ones. + + Parameters + ---------- + clean_mutations : pd.DataFrame + Clean mutation table (all samples), with at least SAMPLE_ID, MUT_ID and VAF. + somatic_mutations : pd.DataFrame + Somatic mutation table (all samples), with at least SAMPLE_ID, MUT_ID and VAF. + + Returns + ------- + pd.DataFrame + Rows of the clean table whose (SAMPLE_ID, MUT_ID) pair is not present in + the somatic table. + """ + merged = clean_mutations.merge( + somatic_mutations[["SAMPLE_ID", "MUT_ID", "VAF"]], + on=["SAMPLE_ID", "MUT_ID"], + suffixes=("", "_somatic"), + how="outer", + ) + germline_mutations = merged[~(merged["VAF_somatic"].notnull())].reset_index(drop=True) + return germline_mutations.drop(columns=["VAF_somatic"]) + + +def assign_genotype_bins(germline_mutations): + """ + Bin the VAF of each germline mutation into genotype classes. + + Bin 0: VAF <= 0, bin 1: (0, 0.25], bin 2: (0.25, 0.75], bin 3: (0.75, 1], + bin 4: VAF > 1. + """ + germline_mutations = germline_mutations.copy() + germline_mutations["genotype_bin"] = np.digitize( + germline_mutations["VAF"], right=True, bins=GENOTYPE_BINS + ) + return germline_mutations + + +def build_pca_matrix(germline_mutations): + """ + Build the samples x variants genotype matrix used for the PCA. + + Only the mutations in the heterozygous-like genotype bins are kept; samples + without a mutation get a 0. + + Returns + ------- + pd.DataFrame + Matrix with one row per sample and one column per germline mutation. + """ + subset = germline_mutations[germline_mutations["genotype_bin"].isin(GENOTYPE_BINS_FOR_PCA)] + pca_data = subset.pivot_table( + index="SAMPLE_ID", columns="MUT_ID", values="genotype_bin", aggfunc="first" + ) + return pca_data.fillna(0) + + +def run_pca(pca_data, number_of_components=NUMBER_OF_COMPONENTS): + """ + Scale the genotype matrix and run the PCA on it. + + Parameters + ---------- + pca_data : pd.DataFrame + Samples x variants genotype matrix. + number_of_components : int, optional + Requested number of principal components. It is capped by the number of + samples and variants available. + + Returns + ------- + pca_df : pd.DataFrame + PC coordinates per sample. + pca : sklearn.decomposition.PCA + Fitted PCA object. + """ + n_components = min(number_of_components, pca_data.shape[0], pca_data.shape[1]) + X_scaled = StandardScaler().fit_transform(pca_data) + pca = PCA(n_components=n_components) + pca_result = pca.fit_transform(X_scaled) + pca_df = pd.DataFrame( + pca_result, + columns=[f"PC{i+1}" for i in range(n_components)], + index=pca_data.index, + ) + return pca_df, pca + + +def plot_vaf_bins(germline_mutations, pdf): + """Plot the VAF distribution of the germline mutations per genotype bin.""" + fig, ax = plt.subplots(figsize=(6, 6)) + sns.boxplot(data=germline_mutations, x="genotype_bin", y="VAF", showfliers=False, ax=ax) + sns.stripplot(data=germline_mutations, x="genotype_bin", y="VAF", ax=ax) + ax.set_xlabel("Genotype bin") + ax.set_ylabel("VAF") + ax.set_title("Germline mutations VAF per genotype bin") + plt.tight_layout() + pdf.savefig(fig) + plt.close(fig) + + +def plot_elbow(pca, pdf): + """Plot the variance explained by each principal component.""" + fig, ax = plt.subplots(figsize=(6, 6)) + ax.plot(pca.explained_variance_ratio_, marker="o", linestyle="-") + ax.set_ylim(0, pca.explained_variance_ratio_.max() + 0.05) + ax.set_xlabel("PC") + ax.set_ylabel("Percent of variance explained") + ax.set_title("PCA elbow plot") + plt.tight_layout() + pdf.savefig(fig) + plt.close(fig) + + +def plot_pca_scatter(pca_df, pca, pdf): + """Scatter plot of PC1 vs PC2 with the samples annotated.""" + component_x, component_y = "PC1", "PC2" + fig, ax = plt.subplots(figsize=(6, 6)) + ax.scatter( + pca_df[component_x], + pca_df[component_y], + s=80, + alpha=0.85, + edgecolors="k", + linewidths=0.3, + ) + for sample in pca_df.index: + ax.annotate( + sample, + (pca_df.loc[sample, component_x], pca_df.loc[sample, component_y]), + fontsize=6, + ha="left", + va="bottom", + ) + ax.set_xlabel(f"{component_x} ({pca.explained_variance_ratio_[0] * 100:.1f}% variance explained)") + ax.set_ylabel(f"{component_y} ({pca.explained_variance_ratio_[1] * 100:.1f}% variance explained)") + ax.set_title("PCA of SNPs within Panel Genes") + plt.tight_layout() + pdf.savefig(fig) + plt.close(fig) + + +def flag_pathogenic_snps(germline_mutations, gnomad_af_threshold): + """ + Flag candidate pathogenic germline SNPs based on the available annotation. + + A germline SNP is considered a candidate pathogenic variant when it complies + with the germline criteria (VAF, vd_VAF and VAF_AM above the germline + threshold), has a VAF greater than 0.25 (otherwise it is not a SNP in any + sense), is protein-affecting (nonsense, missense or essential splice + according to the canonical_Protein_affecting column) and rare in the + population (all the available gnomAD allele frequencies below the + threshold). + + Parameters + ---------- + germline_mutations : pd.DataFrame + Germline mutation table with the VEP annotation columns. + gnomad_af_threshold : float + Maximum gnomAD allele frequency for a variant to be considered rare. + + Returns + ------- + pd.DataFrame + The input table with an additional IS_PATHOGENIC boolean column. + """ + germline_mutations = germline_mutations.copy() + is_protein_affecting = ( + germline_mutations["canonical_Protein_affecting"] == PROTEIN_AFFECTING_VALUE + ) + + available_af_columns = [ + col for col in GNOMAD_AF_COLUMNS if col in germline_mutations.columns + ] + if available_af_columns: + af_values = germline_mutations[available_af_columns].apply(pd.to_numeric, errors="coerce") + is_rare = (af_values.fillna(0) <= gnomad_af_threshold).all(axis="columns") + else: + print("No gnomAD allele frequency columns found; the rarity criterion is skipped") + is_rare = pd.Series(True, index=germline_mutations.index) + + # The germline criteria and the minimum SNP VAF are hard requirements + is_snp = germline_mask(germline_mutations, 0) & (germline_mutations["VAF"] > MIN_SNP_VAF) + germline_mutations["IS_PATHOGENIC"] = is_protein_affecting & is_rare & is_snp + return germline_mutations + + +def summarize_pathogenic_snps(germline_mutations): + """ + Build the per-sample and per-variant tables of candidate pathogenic germline SNPs. + + Returns + ------- + pathogenic_snps : pd.DataFrame + One row per candidate pathogenic germline SNP with its main annotation. + pathogenic_snps_summary : pd.DataFrame + One row per sample with the number of candidate pathogenic germline SNPs. + """ + pathogenic_snps = germline_mutations[germline_mutations["IS_PATHOGENIC"]].reset_index(drop=True) + reported_columns = [col for col in PATHOGENIC_SNP_COLUMNS if col in pathogenic_snps.columns] + pathogenic_snps = pathogenic_snps[reported_columns] + + pathogenic_snps_summary = ( + pathogenic_snps.groupby("SAMPLE_ID", as_index=False) + .size() + .rename(columns={"size": "N_PATHOGENIC_SNPS"}) + ) + return pathogenic_snps, pathogenic_snps_summary + + +def plot_pathogenic_snps(pathogenic_snps_summary, pdf): + """Bar plot with the number of candidate pathogenic germline SNPs per sample.""" + fig, ax = plt.subplots(figsize=(max(6, 0.3 * len(pathogenic_snps_summary)), 6)) + ax.bar(pathogenic_snps_summary["SAMPLE_ID"], pathogenic_snps_summary["N_PATHOGENIC_SNPS"]) + ax.set_xlabel("Sample") + ax.set_ylabel("Number of candidate pathogenic germline SNPs") + ax.set_title("Candidate pathogenic germline SNPs per sample") + ax.set_xticks(range(len(pathogenic_snps_summary))) + ax.set_xticklabels(pathogenic_snps_summary["SAMPLE_ID"], rotation=90, fontsize=6) + plt.tight_layout() + pdf.savefig(fig) + plt.close(fig) + + +def infer_sample_ancestry(germline_mutations): + """ + Infer the most likely ethnic group of each sample from the gnomAD + population-specific allele frequencies. + + For each sample, the mean population-specific gnomAD allele frequency of its + confident germline SNPs (VAF greater than 0.25) is computed per population; + the population with the highest mean allele frequency is reported as the + most likely ethnic group. The rationale is that a sample carrying variants + that are common in a given population is more likely to come from that + population. The OTH (other) population is reported but excluded from the + most likely group assignment since it is not an actual ethnic group. + + Returns + ------- + ancestry_profiles : pd.DataFrame + One row per sample with the mean population-specific gnomAD allele + frequency per population and the MOST_LIKELY_POPULATION column. + """ + population_columns = { + population: column + for population, column in GNOMAD_GENOME_POP_COLUMNS.items() + if column in germline_mutations.columns + } + if not population_columns: + population_columns = { + population: column + for population, column in GNOMAD_EXOME_POP_COLUMNS.items() + if column in germline_mutations.columns + } + if not population_columns: + print("No gnomAD population-specific allele frequency columns found; skipping the ancestry inference") + return pd.DataFrame() + + confident_snps = germline_mutations[germline_mutations["VAF"] > MIN_SNP_VAF] + if confident_snps.empty: + print("No confident germline SNPs (VAF > 0.25) found; skipping the ancestry inference") + return pd.DataFrame() + + af_profiles = confident_snps[list(population_columns.values())].apply(pd.to_numeric, errors="coerce") + ancestry_profiles = ( + af_profiles.assign(SAMPLE_ID=confident_snps["SAMPLE_ID"]) + .groupby("SAMPLE_ID", as_index=False) + .mean() + .rename(columns={column: population for population, column in population_columns.items()}) + ) + + assignable_populations = [pop for pop in population_columns if pop != "OTH"] + ancestry_profiles["MOST_LIKELY_POPULATION"] = ancestry_profiles[assignable_populations].idxmax(axis="columns") + return ancestry_profiles + + +def plot_ancestry_heatmap(ancestry_profiles, pdf): + """Heatmap of the mean population-specific gnomAD AF per sample.""" + populations = [col for col in ancestry_profiles.columns if col not in ("SAMPLE_ID", "MOST_LIKELY_POPULATION")] + heatmap_data = ancestry_profiles.set_index("SAMPLE_ID")[populations] + fig, ax = plt.subplots(figsize=(1.0 + 0.5 * len(populations), 1.0 + 0.4 * len(heatmap_data))) + sns.heatmap(heatmap_data, annot=True, fmt=".4f", cmap="viridis", ax=ax, cbar_kws={"label": "Mean gnomAD population AF"}) + ax.set_xlabel("gnomAD population") + ax.set_ylabel("Sample") + ax.set_title("Ancestry inference from gnomAD population AFs") + plt.tight_layout() + pdf.savefig(fig) + plt.close(fig) + + +@click.command() +@click.option( + "--clean_maf", + type=click.Path(exists=True), + required=True, + help="Path to the clean mutations file (all samples).", +) +@click.option( + "--somatic_maf", + type=click.Path(exists=True), + required=True, + help="Path to the somatic mutations file (all samples).", +) +@click.option( + "--output_prefix", + default="", + show_default=True, + help="Prefix for the output files.", +) +@click.option( + "--gnomad-af-threshold", + default=0.001, + show_default=True, + type=float, + help="Maximum gnomAD allele frequency for a germline SNP to be considered rare.", +) +def main(clean_maf, somatic_maf, output_prefix, gnomad_af_threshold): + """Derive the germline mutations and summarize them with VAF bins and a PCA.""" + clean_mutations = pd.read_table(clean_maf) + somatic_mutations = pd.read_table(somatic_maf) + print(f"Clean mutations: {clean_mutations.shape}") + print(f"Somatic mutations: {somatic_mutations.shape}") + + germline_mutations = assign_genotype_bins( + subtract_somatic_mutations(clean_mutations, somatic_mutations) + ) + print(f"Germline mutations after subtracting the somatic calls: {germline_mutations.shape}") + + germline_mutations = flag_pathogenic_snps(germline_mutations, gnomad_af_threshold) + germline_mutations.to_csv(f"{output_prefix}.germline.mutations.tsv", sep="\t", index=False) + + pathogenic_snps, pathogenic_snps_summary = summarize_pathogenic_snps(germline_mutations) + pathogenic_snps.to_csv(f"{output_prefix}.pathogenic_snps.tsv", sep="\t", index=False) + pathogenic_snps_summary.to_csv(f"{output_prefix}.pathogenic_snps_summary.tsv", sep="\t", index=False) + print(f"Candidate pathogenic germline SNPs: {pathogenic_snps.shape[0]}") + + ancestry_profiles = infer_sample_ancestry(germline_mutations) + if not ancestry_profiles.empty: + ancestry_profiles.to_csv(f"{output_prefix}.ancestry_inference.tsv", sep="\t", index=False) + print(f"Most likely population per sample:\n{ancestry_profiles[['SAMPLE_ID', 'MOST_LIKELY_POPULATION']].to_string(index=False)}") + + with PdfPages(f"{output_prefix}.germline_snps_summary.pdf") as pdf: + if germline_mutations.empty: + print("No germline mutations left after subtracting the somatic calls; skipping the plots") + return + plot_vaf_bins(germline_mutations, pdf) + if not pathogenic_snps_summary.empty: + plot_pathogenic_snps(pathogenic_snps_summary, pdf) + if not ancestry_profiles.empty: + plot_ancestry_heatmap(ancestry_profiles, pdf) + + pca_data = build_pca_matrix(germline_mutations) + if pca_data.empty: + print("No germline mutations in the genotype bins used for the PCA; skipping the PCA") + return + pca_df, pca = run_pca(pca_data) + plot_elbow(pca, pdf) + if pca_df.shape[1] >= 2: + plot_pca_scatter(pca_df, pca, pdf) + else: + print("Only one principal component available; skipping the PC1 vs PC2 scatter") + + +if __name__ == "__main__": + main() diff --git a/conf/results_outputs.config b/conf/results_outputs.config index a821bd13..6703ef65 100644 --- a/conf/results_outputs.config +++ b/conf/results_outputs.config @@ -4,7 +4,6 @@ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ -// TODO revise that all the outputs of these steps are properly copied to the output directory process { withName: ANNOTATEDEPTHS{ @@ -89,6 +88,14 @@ process { pattern: '**{tsv,pdf,png}', ] } + withName: SATURATIONKINETICS{ + publishDir = [ + path: { "${params.outdir}/plots/saturation_kinetics" }, + mode: params.publish_dir_mode, + pattern: '**{tsv,pdf}', + ] + } + withName: PLOTINTERINDIVIDUALVARIABILITY{ publishDir = [ path: { "${params.outdir}/plots/interindividual_variability" }, @@ -407,4 +414,12 @@ process { pattern: '**{tsv,pdf}', ] } + + withName: 'GERMLINE_MUTATIONS' { + publishDir = [ + path: { "${params.outdir}/germline" }, + mode: params.publish_dir_mode, + pattern: '**{tsv,pdf}', + ] + } } diff --git a/conf/tools/hdp_sig_extraction.config b/conf/tools/hdp_sig_extraction.config index 4fe758dc..f94d90c0 100644 --- a/conf/tools/hdp_sig_extraction.config +++ b/conf/tools/hdp_sig_extraction.config @@ -11,7 +11,7 @@ */ params.norm_file = "NA" -params.prior_file = "NA" // TODO this could be a file and if provided we could use it +params.prior_file = "NA" // this could be a file and if provided we could use it params.n_mut_cutoff = 50 params.sig_activity_threshold = 0 params.cohort_threshold = 0 diff --git a/docs/output.md b/docs/output.md index 26c85b78..6c906641 100644 --- a/docs/output.md +++ b/docs/output.md @@ -119,6 +119,7 @@ The directory tree below shows the maximum diversity of outputs the pipeline can │ │ └── chimerax │ ├── gene_subgenic_selection │ ├── saturation_proportions +│ ├── saturation_kinetics │ └── interindividual_variability ├── qc │ ├── trinucleotide_proportions @@ -323,6 +324,7 @@ Optional (subgenic / domain expansion): - Plot basic statistics on numbers and distribution of mutations in genes. - Plot selection results (omega, OncodriveFML, Oncodrive3D, gene/subgenic saturation, interindividual variability). +- Plot saturation kinetics curves: empirical discovery index curves (proportion of mutated sites vs sequencing depth, obtained by downsampling the observed mutations) compared against the theoretical neutral saturation curve expected from the per-site relative mutability and the synonymous mutation rate. Plotting scope can be controlled with `plot_only_allsamples`: when `true`, only cohort-level plots are generated; when `false`, plots are also produced for each defined subgroup. @@ -334,8 +336,21 @@ Plotting scope can be controlled with `plot_only_allsamples`: when `true`, only - `plots/selection/{omega,omegagloballoc,oncodrive3d}/` - `plots/gene_subgenic_selection/` - `plots/saturation_proportions/` +- `plots/saturation_kinetics/` - `plots/interindividual_variability/` +### Saturation kinetics curves + +The `COMPUTE_SATURATION_KINETICS` step quantifies how the fraction of uniquely mutated sites in a gene approaches saturation as sequencing depth increases, and whether the observed approach is faster or slower than expected under neutrality. It requires both `--omega` and one of the mutability-driven analyses (`--oncodrivefml`, `--oncodriveclustl` or `--oncodrive3d`), since it consumes the omega preprocessing mutability table and the relative mutability per site. + +For each group and each (resolution, impact) combination — genomic/residue × protein_affecting/nonsense/truncating/missense/synonymous — it produces: + +- `{group}.curves/{sites}_{impact}_empirical.pdf` — empirical discovery index curves for all genes: proportion of mutated sites as a function of depth, with 95% intervals across downsampling replicates. +- `{group}.curves/{sites}_{impact}_theoretical_empirical.pdf` — the same empirical curves overlaid with the theoretical neutral saturation curve. +- `{group}.curves/{sites}_{impact}_slopes.pdf` — per-gene comparison of the rate of change (Δ proportion / Δ log10 depth) of the empirical curve against the theoretical neutral curve, computed over identical depth intervals (the theoretical curve is interpolated at the empirical depths in log-space). +- `{group}_mutations_{sites}_rates.{impact}.tsv` — per-gene/per-site unique-mutation probabilities at each subsampling depth. +- `{group}_slopes_{sites}.{impact}.tsv` — per-gene interval slopes (`SLOPE_EMPIRICAL`, `SLOPE_THEORETICAL`, `SLOPE_RATIO`) with the depth bounds of each interval and the proportion of positions covered at its midpoint (`PROPORTION_COVERED`), intended for cross-run comparison. + ### Examples ![NeedlePlots](images/needle_plots.png) diff --git a/modules/local/bbgtools/oncodrivefml/main.nf b/modules/local/bbgtools/oncodrivefml/main.nf index f882ec88..10858b60 100644 --- a/modules/local/bbgtools/oncodrivefml/main.nf +++ b/modules/local/bbgtools/oncodrivefml/main.nf @@ -22,7 +22,6 @@ process ONCODRIVEFML { def args = task.ext.args ?: "" // "-s ${params.seed}" def prefix = task.ext.prefix ?: "" prefix = "${meta.id}${prefix}" - // TODO: See if we can provide the entire json as an input parameter """ cat > oncodrivefml_v2.mutability.conf << EOF [genome] diff --git a/modules/local/createmaskmatrix/main.nf b/modules/local/createmaskmatrix/main.nf index 2a796369..9909090e 100644 --- a/modules/local/createmaskmatrix/main.nf +++ b/modules/local/createmaskmatrix/main.nf @@ -3,8 +3,7 @@ process CREATE_MASK_MATRIX { label 'cpu_low' label 'mem_low' - - container "docker.io/bbglab/deepcsa-core:0.1.0" + label 'deepcsa_core' input: path(bed_files) // List of sample-specific flagged position BED files diff --git a/modules/local/dna2protein/main.nf b/modules/local/dna2protein/main.nf index f49d4fea..b89820c2 100644 --- a/modules/local/dna2protein/main.nf +++ b/modules/local/dna2protein/main.nf @@ -14,6 +14,7 @@ process DNA_2_PROTEIN_MAPPING { output: tuple val(meta2), path("depths_per_position_exon_gene.tsv") , emit: depths_exons_positions tuple val(meta2), path("panel_exons.bed4.bed") , emit: panel_exons_bed + tuple val(meta2), path("panel_exons_protein_intervals.tsv") , emit: panel_exons_protein_intervals tuple val(meta2), path("*.pdf") , emit: covered_genes_pdf tuple val(meta2), path("coverage*.tsv") , emit: covered_genes_tsv path "versions.yml" , topic: versions diff --git a/modules/local/germlinemuts/main.nf b/modules/local/germlinemuts/main.nf new file mode 100644 index 00000000..e595ded0 --- /dev/null +++ b/modules/local/germlinemuts/main.nf @@ -0,0 +1,49 @@ +process GERMLINE_MUTATIONS { + + tag "$meta.id" + label 'deepcsa_core' + + input: + tuple val(meta), path(clean_maf) + tuple val(meta2), path(somatic_maf) + + output: + tuple val(meta), path("*.germline.mutations.tsv"), emit: germline_mutations + tuple val(meta), path("*pathogenic_snps*.tsv") , emit: pathogenic_snps + tuple val(meta), path("*.ancestry_inference.tsv"), emit: ancestry + tuple val(meta), path("*.pdf") , emit: plots + path "versions.yml" , topic: versions + + script: + def prefix = task.ext.prefix ?: "" + prefix = "${meta.id}${prefix}" + def gnomad_af_threshold = task.ext.gnomad_af_threshold ? "--gnomad-af-threshold ${task.ext.gnomad_af_threshold}" : "" + """ + summarize_germline_snps.py \\ + --clean_maf ${clean_maf} \\ + --somatic_maf ${somatic_maf} \\ + --output_prefix ${prefix} \\ + ${gnomad_af_threshold} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python --version | sed 's/Python //g') + END_VERSIONS + """ + + stub: + def prefix = task.ext.prefix ?: "" + prefix = "${meta.id}${prefix}" + """ + touch ${prefix}.germline.mutations.tsv + touch ${prefix}.pathogenic_snps.tsv + touch ${prefix}.pathogenic_snps_summary.tsv + touch ${prefix}.ancestry_inference.tsv + touch ${prefix}.germline_snps_summary.pdf + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python --version | sed 's/Python //g') + END_VERSIONS + """ +} diff --git a/modules/local/plot/saturation/main.nf b/modules/local/plot/saturation/gene_distribution/main.nf similarity index 100% rename from modules/local/plot/saturation/main.nf rename to modules/local/plot/saturation/gene_distribution/main.nf diff --git a/modules/local/plot/saturation/meta.yml b/modules/local/plot/saturation/meta.yml deleted file mode 100644 index e69de29b..00000000 diff --git a/modules/local/saturation_kinetics/compute/main.nf b/modules/local/saturation_kinetics/compute/main.nf new file mode 100644 index 00000000..2542b9c5 --- /dev/null +++ b/modules/local/saturation_kinetics/compute/main.nf @@ -0,0 +1,50 @@ +process COMPUTE_SATURATION_KINETICS { + + tag "$meta.id" + label 'cpu_low' + + container 'docker.io/ferriolcalvet/saturation:v0.1.0' + + input: + tuple val(meta) , path(mutations), path(depths), path(omega_mutability), path(relative_mutability), path(relative_mutability_index) + tuple val(meta1), path(captured_panel_rich) + tuple val(meta2), path(expanded_panel, stageAs: "expanded_panel.tsv") + + + output: + tuple val(meta), path("**.tsv"), emit: table + tuple val(meta), path("**.pdf"), emit: plots + path "versions.yml" , topic: versions + + + script: + def prefix = task.ext.prefix ?: "" + prefix = "${meta.id}${prefix}" + """ + saturation_kinetics_curves.py \\ + --somatic-mutations-file ${mutations} \\ + --vep-file ${captured_panel_rich} \\ + --consensus-panel-file ${expanded_panel} \\ + --depths-file ${depths} \\ + --omega-mutability-file ${omega_mutability} \\ + --relative-mutability-file ${relative_mutability} \\ + --group-name ${meta.id} + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python --version | sed 's/Python //g') + END_VERSIONS + """ + + stub: + def prefix = task.ext.prefix ?: "" + prefix = "${meta.id}${prefix}" + """ + touch ${prefix}.tsv + + cat <<-END_VERSIONS > versions.yml + "${task.process}": + python: \$(python --version | sed 's/Python //g') + END_VERSIONS + """ +} diff --git a/modules/local/sortpanel/main.nf b/modules/local/sortpanel/main.nf index 87af6817..47853e43 100644 --- a/modules/local/sortpanel/main.nf +++ b/modules/local/sortpanel/main.nf @@ -3,8 +3,7 @@ process SORT_MERGED_PANEL { tag "${meta.id}" label 'mem_low' - - container "docker.io/bbglab/deepcsa-core:0.0.2-alpha" + label 'deepcsa_core' input: tuple val(meta), path(panel) diff --git a/subworkflows/local/enrichpanels/main.nf b/subworkflows/local/enrichpanels/main.nf index de56bcf1..9221d297 100644 --- a/subworkflows/local/enrichpanels/main.nf +++ b/subworkflows/local/enrichpanels/main.nf @@ -80,4 +80,5 @@ workflow ENRICHPANELS { dna2protein_mapping_depth_exons = DNA2PROTEINMAPPING.out.depths_exons_positions.first() dna2protein_mapping_panel_exons = DNA2PROTEINMAPPING.out.panel_exons_bed.first() + dna2protein_mapping_panel_exons_protein = DNA2PROTEINMAPPING.out.panel_exons_protein_intervals.first() } diff --git a/subworkflows/local/mutationpreprocessing/main.nf b/subworkflows/local/mutationpreprocessing/main.nf index 789903e2..e7daddaa 100644 --- a/subworkflows/local/mutationpreprocessing/main.nf +++ b/subworkflows/local/mutationpreprocessing/main.nf @@ -21,6 +21,7 @@ include { PLOT_MUTATIONS as PLOTSOMATICMAF } from '../../../m include { PLOT_NEEDLES as PLOTNEEDLES } from '../../../modules/local/plot/needles/main' include { DOWNSAMPLE_MUTATIONS as DOWNSAMPLEMUTS } from '../../../modules/local/downsample/mutations/main' include { COMPUTE_CONTAMINATION as CONTAMINATION } from '../../../modules/local/contamination/main' +include { GERMLINE_MUTATIONS as GERMLINEMUTATIONS } from '../../../modules/local/germlinemuts/main' workflow MUTATION_PREPROCESSING { @@ -126,6 +127,9 @@ workflow MUTATION_PREPROCESSING { // Clean mutations based on artifact filtering decisions CLEANMUTATIONS(all_clean_mutations) + channel.of([["id": "all_samples"]]) + .join(CLEANMUTATIONS.out.mutations).first() + .set{clean_muts_all_samples} // Keep only somatic mutations SOMATICMUTATIONS(CLEANMUTATIONS.out.mutations) @@ -158,6 +162,8 @@ workflow MUTATION_PREPROCESSING { PLOTNEEDLES(muts_for_plotting, sequence_information_df) + GERMLINEMUTATIONS(clean_muts_all_samples, muts_all_samples) + // Compile a BED file with all the mutations that are discarded due to: // Other sample SNP // All sites with this filter should be remove from the background. @@ -180,5 +186,6 @@ workflow MUTATION_PREPROCESSING { mutations_all_samples = muts_all_samples all_raw_vep_annotation = SUMANNOTATION.out.tab_all bedfile_clean = bedfile_updated + clean_maf_all_samples = clean_muts_all_samples } diff --git a/subworkflows/local/omega/main.nf b/subworkflows/local/omega/main.nf index a59aa2be..d97cbb4f 100644 --- a/subworkflows/local/omega/main.nf +++ b/subworkflows/local/omega/main.nf @@ -224,6 +224,7 @@ workflow OMEGA_ANALYSIS{ results_global = global_loc_results expanded_panel = expanded_panel site_comparison = site_comparison_results_flattened + preprocessing_mutab = PREPROCESSING.out.mutabs_n_mutations_tsv.map{ it -> [it[0], it[1]]} all_compiled = all_results all_globalloc_compiled = all_gloc_results diff --git a/subworkflows/local/plottingsummary/main.nf b/subworkflows/local/plottingsummary/main.nf index 70c3f3be..134f2fa2 100644 --- a/subworkflows/local/plottingsummary/main.nf +++ b/subworkflows/local/plottingsummary/main.nf @@ -1,11 +1,12 @@ + include { PLOT_SELECTION_METRICS as PLOTSELECTION } from '../../../modules/local/plot/selection_metrics/main' -include { PLOT_SATURATION as PLOTSATURATION } from '../../../modules/local/plot/saturation/main' +include { PLOT_SATURATION as PLOTSATURATION } from '../../../modules/local/plot/saturation/gene_distribution/main' include { PLOT_SATURATION_PROPORTIONS as PLOTSATURATIONPROPORTIONS } from '../../../modules/local/plot/saturation/proportions/main' include { PLOT_INTERINDIVIDUAL_VARIABILITY as PLOTINTERINDIVIDUALVARIABILITY } from '../../../modules/local/plot/interindividual_variability/main' - +include { COMPUTE_SATURATION_KINETICS as SATURATIONKINETICS } from '../../../modules/local/saturation_kinetics/compute/main' workflow PLOTTING_SUMMARY { @@ -26,6 +27,10 @@ workflow PLOTTING_SUMMARY { domain_df exons_depths_df groups_channel + depths_indv + relative_mutability + omega_mutabilities + all_clean_mutations main: @@ -68,17 +73,21 @@ workflow PLOTTING_SUMMARY { .set{ groups_results_sites } PLOTSELECTION(groups_results, seqinfo_df) - // needles with consequence type - // plot selection at cohort/group level, all the different methods available - // plot selection per domain at cohort level + PLOTSATURATION(groups_results_sites, all_samples_depth, panel, seqinfo_df, pdb_tool_df, domain_df, exons_depths_df) PLOTSATURATIONPROPORTIONS(groups_mutations, panel, full_panel_rich, expanded_panel) - // plot gene + site selection - // omega selection per domain in gene - // ? plot saturation kinetics curves + + + // plot saturation kinetics curves + groups_mutations + .join(depths_indv) + .join(omega_mutabilities) + .join(relative_mutability) + .set{ groups_mutations_depths_n_mutability } + SATURATIONKINETICS(groups_mutations_depths_n_mutability, full_panel_rich, expanded_panel) PLOTINTERINDIVIDUALVARIABILITY(samples, all_groups, panel, all_mutdensities, all_mutdensities_adjusted) diff --git a/workflows/deepcsa.nf b/workflows/deepcsa.nf index ec5babee..a281302b 100644 --- a/workflows/deepcsa.nf +++ b/workflows/deepcsa.nf @@ -165,6 +165,8 @@ workflow DEEPCSA { all_adjusted_mutdensities_file = channel.value(file("${projectDir}/assets/placeholder_no_file.tsv", checkIfExists: true)) all_compiled_stabilities = channel.empty() dndscv_table = channel.empty() + relative_mutabilities = channel.empty() + omega_mutabilities = channel.empty() // if the user wants to use custom gene groups, import the gene groups table // otherwise I am using the input csv as a dummy value channel @@ -417,6 +419,7 @@ workflow DEEPCSA { MUTPROFILEALL.out.profile, CREATEPANELS.out.exons_consensus_panel ) + relative_mutabilities = MUTABILITYALL.out.mutability if (params.profilenonprot){ MUTABILITYNONPROT(mutations_in_exons, annotated_depths, @@ -502,6 +505,7 @@ workflow DEEPCSA { CREATEPANELS.out.exons_consensus_panel, group_keys_ch ) + omega_mutabilities = OMEGA.out.preprocessing_mutab positive_selection_results = positive_selection_results.join(OMEGA.out.results, remainder: true) all_compiled_omegas = OMEGA.out.all_compiled if (params.omega_mutabilities){ @@ -694,7 +698,11 @@ workflow DEEPCSA { seqinfo_df, CREATEPANELS.out.domains_in_panel, ENRICHPANELS.out.dna2protein_mapping_depth_exons, - group_keys_ch + group_keys_ch, + annotated_depths, + relative_mutabilities, + omega_mutabilities, + MUT_PREPROCESSING.out.clean_maf_all_samples ) }