cfWGS-MM-MRD

Analysis code for "Cell-free DNA Whole Genome Sequencing for Non-Invasive Minimal Residual Disease Detection in Multiple Myeloma".

The repository contains the R scripts used to process the clinical and cfWGS data, construct the analysis tables, fit and evaluate the MRD classifiers, and produce the main and supplementary figures and tables. The scripts retain their original analysis-based numbering so that saved objects and manuscript records can be traced back to the code that generated them.

The files are not one simple linear pipeline: several scripts branch into mutation, copy-number, fragmentomics, MRDetect, dilution-series, longitudinal, and survival analyses. The tables below describe the purpose, inputs, outputs, and manuscript contribution of each script. The docs/SCRIPT_GUIDE.md explains how the scripts connect, and docs/FIGURE_TABLE_MAP.md provides a concise panel-level crosswalk.

The analysis can be run as documented individual stages or through the repository runner. The 50-repeat patient-grouped nested cross-validation is run separately because it consists of three long, parameterized model blocks.

Raw data files are not included in this repository (see Data availability below).

Run commands from the repository root. Before running a script, stage the protected inputs listed in its header and confirm that the preceding stage has produced the required objects. Use Pipeline scripts for the source-data and figure stages and Patient-grouped repeated nested cross-validation for the analysis used in Figures 3A and 4A, Extended Data Figures 5A, 7A and 9A-B, and Supplementary Table 4.

Check the expected local inputs and inspect the planned script order without running an analysis:

Rscript prepare_local_inputs.R --check
Rscript run_pipeline.R

The input check is read-only unless --report <path> or --copy-missing is supplied.

run_pipeline.R is a read-only dry run unless --execute is supplied. Its ordered script list is config/source_pipeline.tsv. The routine model-scoring stage applies the preserved models and thresholds; it does not retrain them.

Several numbered scripts also call manuscript_output_helpers.R to place a labeled copy of a figure, table, or source-data file under the ignored final_manuscript_objects/ directory. This helper requires the artifact map in docs/manuscript_artifact_source_map.tsv.

The numbered scripts are still organized by analysis question, not by final figure number. Each script header lists the same outputs shown below. When a script runs, the files used in the manuscript are copied or saved into final_manuscript_objects/ under the final figure/table label.

Each numbered script header lists the manuscript figure/table outputs that the script creates or updates. Use docs/FIGURE_TABLE_MAP.md for the current, manually reviewed panel-to-script map.

The stage numbers describe the broad order of the analysis. Scripts within a stage may form separate branches, and optional/support scripts are identified below. Run commands from the repository root.

For the platform-calibration panels, run 1_7D before 1_8E and run 1_8C before 1_8D.

Figure 3D and Figure 4C use the same cohort assignment but different frontline eligibility rules. Figure 3D includes all 42 landmark samples with an evaluable BM-informed cfWGS call; 39 of those also have MFC or clonoSEQ data. Figure 4C first requires an evaluable blood-informed cfWGS call and at least one of MFC or clonoSEQ, leaving 41 of 46 blood-call-evaluable landmark samples.

The panel reconstructions in 5_1 should be compared with the final figures before regenerated workbooks replace the reviewed manuscript copies. The currently retained final workbooks contain 17 main-figure sheets and 54 Extended Data sheets. The generated 18/56-sheet workbooks additionally contain Fig1B, ED10A, and ED10B, which are absent from the retained final copies. Figure 1B is assembled from exported counts, and the final Extended Data Figure 10 tracks come from the separate ichorCNA plotting workflow.

The repository generates the individual panels and their source data. Final multi-panel figure assembly is performed separately and assembly scripts are not included.

Final_Scripts/
├── 1_0_*.R … 4_2_*.R      ← clinical, feature, figure, model, and outcome scripts
├── 5_1_*.R … 5_4_*.R      ← figure source-data export, workbook build, and validation
├── 6_12_*.R … 6_17_*.R    ← 50-repeat grouped nested-CV and plotting scripts
├── config/source_pipeline.tsv ← ordered non-CV runner plan
├── run_pipeline.R          ← dry-run and execution entry point
├── helpers.R               ← shared utility functions
├── manuscript_output_helpers.R ← labeled figure/table copy helper
├── docs/                    ← script guide and figure/table crosswalk
└── README.md               ← this file

Prerequisites

• R ≥ 4.2
• Key libraries include tidyverse, readxl, data.table, lubridate, ComplexHeatmap, circlize, ChromHeatMap, maftools, GenomicRanges, gtsummary, officer, flextable, gt, pROC, patchwork, rmda, exact2x2, survival, survminer, tableone, ggbreak, ggridges, ggpubr, GGally, pbapply, PRROC, VennDiagram, caret, glmnet, Matrix, DescTools, GeneCycle, RColorBrewer, fuzzyjoin, rstatix, openxlsx, scales, viridis, and writexl. See each script header for its dependencies.

Source-script workflow

1. Obtain the required inputs described in Data availability and place them at the project-relative paths listed in each script header.
2. Run from the repository root. Scripts use the project-relative paths listed in their headers.
3. Follow the stage table above and the input/output header in each script. Several branches are independent, so only run the stages needed for the intended figure or table.
4. Run the repeated nested-CV workflow below separately when reproducing its figures and Supplementary Table 4.

To inspect the retained non-CV workflow without writing files:

Rscript run_pipeline.R

To execute that workflow after all protected inputs have been staged:

Rscript run_pipeline.R --execute

After completing the grouped-CV commands below and generating their panels, rebuild and validate the figure source-data workbooks with:

Rscript run_pipeline.R --execute --only 5_3 --include-post-cv

The manuscript's internal-validation analysis assigns patients—not individual sample rows—to both outer and inner folds. Thus, all repeated samples belonging to a patient remain on the same side of every split. This workflow is separate from the saved full-training models and thresholds used to score the test cohort; it does not overwrite those test-scoring objects.

The paper uses the complete 32-model library, five outer folds repeated 50 times, five inner folds repeated five times, 2,000 patient-clustered bootstrap replicates, and base seed 20260731. These are now the script defaults. The commands below still specify them explicitly so the saved command records the complete resampling design.

Run the three model blocks from the repository root with new, unused run IDs:

Rscript 6_12_Patient_Grouped_Repeated_Nested_CV.R \
  --run-id paper_50rep_bm_v1 --model-library all \
  --models BM_Sites,BM_cVAF,BM_Raw_cVAF,BM_All_Mutation_Features,BM_Combined_Mutation_Zscores,BM_Mutation_Fragmentomics_Full,BM_Mutation_Fragmentomics_Min,Fragmentomics_BMRestricted_Full,Fragmentomics_BMRestricted_Min,Fragmentomics_BMRestricted_FS,Fragmentomics_BMRestricted_Mean_Coverage,Fragmentomics_BMRestricted_Proportion_Short,Fragmentomics_BMRestricted_Tumor_Fraction \
  --outer-repeats 50 --inner-repeats 5 --bootstrap-reps 2000 \
  --seed 20260731

Rscript 6_12_Patient_Grouped_Repeated_Nested_CV.R \
  --run-id paper_50rep_blood_v1 --model-library all \
  --models Blood_Sites,Blood_cVAF,Blood_Raw_cVAF,Blood_All_Mutation_Features,Blood_Combined_Mutation_Zscores,Blood_Mutation_Fragmentomics_Full,Blood_Mutation_Fragmentomics_Min,Fragmentomics_BloodRestricted_Full,Fragmentomics_BloodRestricted_Min,Fragmentomics_BloodRestricted_FS,Fragmentomics_BloodRestricted_Mean_Coverage,Fragmentomics_BloodRestricted_Proportion_Short,Fragmentomics_BloodRestricted_Tumor_Fraction \
  --outer-repeats 50 --inner-repeats 5 --bootstrap-reps 2000 \
  --seed 20260731

Rscript 6_12_Patient_Grouped_Repeated_Nested_CV.R \
  --run-id paper_50rep_fullfrag_v1 --model-library all \
  --models Fragmentomics_FullCohort_Full,Fragmentomics_FullCohort_Min,Fragmentomics_FullCohort_FS,Fragmentomics_FullCohort_Mean_Coverage,Fragmentomics_FullCohort_Proportion_Short,Fragmentomics_FullCohort_Tumor_Fraction \
  --outer-repeats 50 --inner-repeats 5 --bootstrap-reps 2000 \
  --seed 20260731

Each run writes all fold assignments, held-out predictions, fitted settings, thresholds, warnings, summaries, session information, and QC results under Output_tables_2025/patient_grouped_repeated_nested_cv/<run-id>/. It refuses to overwrite an existing run directory and writes RUN_COMPLETE only after its checks pass. The preserved model/validation input files listed in the 6_12 header must be present and non-writable.

After all three blocks complete, assemble the 32-model result:

Rscript 6_13_Assemble_All_Model_Grouped_CV_Results.R \
  --source-runs paper_50rep_bm_v1,paper_50rep_blood_v1,paper_50rep_fullfrag_v1 \
  --output-run-id paper_50rep_combined_v1

Generate the main ROC panels and the Extended Data operating-point panels from that completed combined result:

Rscript 6_17_Generate_Compact_All_Model_Grouped_CV_ROC_Panels.R \
  --input-run-id paper_50rep_combined_v1 \
  --output-run-id paper_50rep_main_roc_v1

Rscript 6_14_Generate_Grouped_CV_Manuscript_Replacement_Panels.R \
  --input-run-id=paper_50rep_combined_v1 \
  --output-run-id=paper_50rep_extended_data_v1

The 6_14_Generate_* script uses --name=value arguments; the preceding 6_12, 6_13, and 6_17 commands use --name value arguments.

• The grouped-CV AUC estimates internal discrimination for a previously unseen patient, conditional on the already-selected feature specification.
• Inner folds tune glmnet and generate out-of-fold predictions for threshold selection; outer-held-out patients are used only for evaluation.
• Patient-clustered bootstrap intervals preserve all repeated samples from a resampled patient.
• The model specifications were selected before this repeated validation; the grouped-CV AUC therefore evaluates those fixed specifications rather than a new feature-discovery process.
• None of these scripts overwrites the saved full-training models, thresholds, older CV objects, or test-cohort predictions.

Outputs are written to new versioned directories under Output_tables_2025/patient_grouped_repeated_nested_cv/ and Output_figures_2025/patient_grouped_repeated_nested_cv/.

Because the input files contain protected patient information they are not distributed with this repository. This includes:

• Clinical metadata spreadsheets (SPORE, M4, IMMAGINE cohorts)
• WGS MAF files, ichorCNA segmentation outputs, and Ig-caller structural variant calls
• MRDetect cumulative-VAF CSV files
• Fragmentomics nucleosome-distance and fragment-score files
• EasyM proteomic MRD data

To reproduce the analysis, request the de-identified data from the study authors.

The pipeline produces:

• CSV / XLSX tables - cleaned metadata, mutation counts, MRD feature matrices, concordance statistics, sensitivity/specificity tables (written to Output_tables_2025/)
• PDF / PNG figures - the manuscript panels generated in R (primarily under Final Tables and Figures/ and Output_figures_2025/); Figure 1B is assembled from an exported source table, and the final Extended Data Figure 10 genome-wide tracks come from a separate ichorCNA plotting workflow
• Support-only QA outputs - non-manuscript review exports are kept in explicitly named support folders such as Final Tables and Figures/swim_plot_support/ and Output_tables_2025_updated/support_only_bam_archive_diagnostics/
• RDS objects - intermediate R data objects passed between scripts

For questions about the analysis, contact Dory Abelman or open an issue on this repository.