big refactor

[matthewpeterkort]

Changelog: Battenberg Refactor v3.0.2

This document outlines the significant mathematical, technical, and feature changes introduced in the refactored Battenberg pipeline compared to the original implementation.

1. Architectural Changes: "Pure R" Pipeline

• Removal of External Executable Orchestration:
  • Change: The pipeline no longer internally manages or invokes external binary executables such as alleleCounter, impute2, or java -jar beagle.jar.
  • Impact: Responsibilities such as Allele Counting and Phasing/Imputation have been factored out. The R package now strictly acts as a consumer of standard genomic file formats (BAM outputs, VCFs) generated by upstream workflow managers (e.g., Nextflow, Snakemake). This makes the package significantly lighter, more portable, and easier to containerize, as it no longer requires a complex "fat" environment with legacy binaries.
• VCF Consumption Model:
  • Change: Instead of wrapping the execution of Beagle 5, the pipeline now includes a native VCF parser (convert_beagle_to_impute) that digests the output of modern phasing tools.
  • Rationale: Decouples the statistical copy number calling from the specific version or implementation of the phasing tool.

2. Mathematical & Algorithmic Changes

Geometric Centroid Selection

• Change: When multiple valid copy number solutions (optima) are found during the grid search (rho/psi space), the algorithm now calculates the geometric centroid of all valid solutions and selects the specific optimum closest to this center.
• Rationale: The original implementation often defaulted to the solution with the absolute highest goodness-of-fit. Centroid selection ensures a more robust "central" parameter set is chosen when the solution space is flat, reducing outlier artifacts.

Recalculated Ploidy (Psi_t)

• Change: In run_clonal_ASCAT, the ploidy (psi) is partially recalculated using only high-confidence clonal segments after the initial grid search.
• Rationale: Ensures the final ploidy estimate is not skewed by subclonal noise.

Optimized Integer Copy Number Fitting

• Change: The C++ implementation (ascat_distance.cpp) explicitly optimizes the integer combination of Major/Minor alleles (checking Floor/Ceiling combinations) for every grid point to minimize the BAF squared error.
• Rationale: Provides a rigorously optimized "best fit" for integer copy numbers at every hypothetical grid point.

LOH/Deletion Constraint Logic

• Change: Applied a heuristic in calculate_solution_fast that valid solutions typically require at least some Loss-Of-Heterozygosity (LOH) or deletions (CN=0).
• Rationale: Filters out high-ploidy artifact solutions that lack biological deletion events.

Winsorization in Segmentation

• Change: Added copynumber::winsorize step prior to PCF segmentation (segmentation.R).
• Rationale: Prevents single-point outliers from distorting segment means.

3. Technical & Performance Optimizations

C++ Acceleration (Rcpp)

• Change: Core bottlenecks—distance calculations (calculate_ascat_dist_matrix_cpp) and segmentation (pcf_core.cpp)—ported to C++.
• Impact: 10x-100x speedup in grid search and segmentation.

High-Performance IO & Vectorization

• Change:
  • vroom: Used for instant reading/merging of large BAF files (haplotype.R).
  • collapse: Replaced base R stats with collapse::fsum/fmean.
  • Vectorization: Grid search logic fully vectorized to remove nested R loops.
• Impact: Massive reduction in I/O overhead and compute time.

Workflow Checkpoints & Resume

• Change: Added preprocessed_data_dir and phasing_results_dir arguments to battenberg().
• Impact: Allows skipping expensive upstream steps (allele counting, imputation) when re-running segmentation with new parameters.

Dynamic Thread Budgeting

• Change: Explicit calculation of total threads (chromosomes_in_parallel × threads_per_chromosome) and dynamic setting of OMP_NUM_THREADS and MKL_NUM_THREADS.
• Impact: Prevents system lockups due to thread oversubscription.

Memory Management

• Change:
  • Chunked Back-Transformation: Processes LRR/BAF vectors in chunks to prevent OOM.
  • Smart Downsampling: Plotting functions downsample data to ~500k points to prevent graphics device hangs.

4. Stability & Robustness

Critical Probe Misalignment Fix

• Change: Explicit logic in run_ascat_enhanced to subset vectors using names (lrr[names(baf)]) rather than position.
• Impact: Fixes data corruption caused by upstream filtering offsets.

Failsafe Grid Search

• Change: Automatic fallback to Full Grid Search or Top-N Search if the optimized local minima search yields no results.

5. New Features

Configurable Grid Search

• New Arguments:
  • n_neighbors_search: Limit search to top N closest points.
  • psi_step / rho_step: Custom grid resolution.
  • local_min_window_size: Adjustable local minima window.

Structured Logging

• Change: Full migration to logger package for timestamped, leveled logs (INFO, DEBUG, ERROR).

Early Termination

• Change: early_termination = TRUE flag to stop grid search once a high-quality solution is found.

6. Code Structure & Dependencies

Modularization

• Change: Split monolithic scripts (clonal_ascat.R) into focused modules (run_clonal_ascat.R, clonal_ascat_calc.R, etc.).

New Dependencies

• Rcpp, RcppRoll (Acceleration)
• collapse, vroom, tictoc, fs (IO/Perf)
• logger, cli, optparse (Interface)
• S4Vectors, IRanges (Genomics)

License

• Change: Updated to AGPL-3 since that is what was listed in the LICENSE file.

7. Container & CLI Support

• Change: Added Dockerfile, Makefile, and cli/optparse support.
• Impact: Facilitates robust command-line usage and reproducible containerized deployment.
• Registry: A Singularity-compatible Docker image is available at quay.io/ohsu-comp-bio/battenberg.