diff --git a/.github/workflows/basic_checks.yaml b/.github/workflows/basic_checks.yaml index 9bc5d72..40e90cb 100644 --- a/.github/workflows/basic_checks.yaml +++ b/.github/workflows/basic_checks.yaml @@ -49,7 +49,14 @@ jobs: - name: Build pkgdown run: | - PATH=$PATH:$HOME/bin/ Rscript -e 'pkgdown::build_site(".")' + PATH=$PATH:$HOME/bin/ Rscript -e 'pkgdown::build_site(".")' + + - name: Add Joint-RPCA validation report + run: | + mkdir -p docs/filterRPCAInput_validation + cp \ + vignettes/filterRPCAInput_validation/joint-rpca-mia-gemelli-comparison.html \ + docs/filterRPCAInput_validation/index.html - name: Upload pkgdown artifact if: github.event_name == 'push' @@ -143,4 +150,4 @@ jobs: push: ${{ github.event_name != 'pull_request' }} tags: | ${{ env.IMAGE }}:latest - ${{ env.IMAGE }}:${{ env.GIT_SHA }} \ No newline at end of file + ${{ env.IMAGE }}:${{ env.GIT_SHA }} diff --git a/vignettes/filterRPCAInput_validation/SETUP.md b/vignettes/filterRPCAInput_validation/SETUP.md new file mode 100644 index 0000000..d8f252d --- /dev/null +++ b/vignettes/filterRPCAInput_validation/SETUP.md @@ -0,0 +1,242 @@ +# Reproduce the mia–Gemelli Joint-RPCA validation + +This guide runs the numbered scripts and then renders one complete HTML report. +You only need to enter the locations on your computer once in each environment. + +## Reproducibility specification + +| Item | Fixed value | +|---|---| +| mia version | `1.21.6` | +| mia commit | `25fa7bb51792566131f2bd8d07b4ebd1fde54b38` | +| mia data | Revised `ibdmdb` main MGX and MTX experiments | +| Run configuration | `config.tsv` | +| Feature prevalence | `0.10` in R = `10%` in Gemelli | +| Components | `3` | +| Maximum iterations | `5` | +| Test samples | `10` | +| Random seed | `42` | + +The R runner checks the mia version and commit. The Python runner checks the +SHA-256 hashes of the inputs exported by R. Both implementations use the same +saved training/test split. + +## Required software + +- Git +- R and RStudio +- Quarto, normally included with RStudio +- Windows Subsystem for Linux (WSL) +- Miniforge or Miniconda inside WSL + +## Project layout + +Extract or copy the supplied files into one folder. This folder is called +`PROJECT_DIR` throughout this guide: + +```text +/ +├── config.tsv +├── SETUP.md +├── joint-rpca-mia-gemelli-comparison.qmd +└── scripts/ + ├── 00_check_gemelli_environment.py + ├── 01_run_mia_filter_joint_rpca.R + ├── 02_run_gemelli_filter_joint_rpca.py + └── 03_compare_mia_gemelli_corrected.R +``` + +The analysis creates `input/` and `results/` automatically. + +## Step 1: enter your Windows paths once + +Open PowerShell. Replace the three example paths below with locations on your +computer. You should only need to edit these three lines: + +```powershell +$PROJECT_DIR = "C:\path\to\filterRPCAInput_validation" +$MIA_CLONE_DIR = "C:\path\to\your\mia" +$MIA_VALIDATION_DIR = "C:\path\to\mia_validation_25fa7bb5" +``` + +- `PROJECT_DIR` is the folder containing this file and `config.tsv`. +- `MIA_CLONE_DIR` is your existing Git clone of the mia repository. +- `MIA_VALIDATION_DIR` is a new location for the clean, pinned worktree. + +Keep this PowerShell window open so the variables remain available. + +## Step 2: create a clean mia checkout + +Do not modify or discard changes in your normal mia checkout. Create a separate +worktree at the pinned commit: + +```powershell +Set-Location $MIA_CLONE_DIR +git fetch upstream +git worktree add --detach $MIA_VALIDATION_DIR 25fa7bb51792566131f2bd8d07b4ebd1fde54b38 +git -C $MIA_VALIDATION_DIR status --short +git -C $MIA_VALIDATION_DIR rev-parse HEAD +``` + +The status command should print nothing. The final command should print the +pinned commit shown above. + +## Step 3: prepare R and run mia + +Open RStudio. Enter your two relevant Windows paths once, using forward slashes +in R. Replace only the first two path strings: + +```r +project_dir <- normalizePath( + "C:/path/to/filterRPCAInput_validation", + winslash = "/" +) +mia_repo <- normalizePath( + "C:/path/to/mia_validation_25fa7bb5", + winslash = "/" +) + +Sys.setenv( + MIA_GEMELLI_PROJECT = project_dir, + MIA_REPO = mia_repo +) +``` + +Install the workflow packages and any missing dependencies of the pinned mia +checkout. This is needed only once for the selected R library: + +```r +install.packages(c("devtools", "digest", "knitr")) +devtools::install_deps(mia_repo, dependencies = TRUE) +``` + +Run mia and export the common inputs: + +```r +source(file.path( + project_dir, + "scripts", + "01_run_mia_filter_joint_rpca.R" +)) +``` + +A successful run ends with `mia run completed successfully.` and creates: + +```text +input/frequency_010/ +results/frequency_010/mia/ +``` + +## Step 4: prepare Gemelli in WSL + +From PowerShell, enter WSL: + +```powershell +wsl +``` + +Windows and WSL name the same folder differently. For example, Windows drive +`C:` is `/mnt/c` in WSL, while drive `D:` is `/mnt/d`. Set the WSL form of your +project path once, replacing only this line: + +```bash +PROJECT_DIR_WSL="/mnt/c/path/to/filterRPCAInput_validation" +cd "$PROJECT_DIR_WSL" +``` + +Enable Conda and activate the environment: + +```bash +source ~/miniforge3/etc/profile.d/conda.sh +conda activate mia-gemelli +``` + +If the environment does not exist, create it and install the pinned Gemelli +source instead: + +```bash +conda create -n mia-gemelli python=3.10 pip -y +conda activate mia-gemelli +python -m pip install \ + "gemelli @ git+https://github.com/biocore/gemelli.git@c53c9ee958948683bc4216cf5fa447e2c3ac4806" +``` + +Run the environment/input check and Gemelli: + +```bash +python scripts/00_check_gemelli_environment.py --project-dir . +python scripts/02_run_gemelli_filter_joint_rpca.py --project-dir . +``` + +Run the second command only after the first passes. A successful run creates +`results/frequency_010/gemelli/`. Then leave WSL with `exit`. + +## Step 5: compare the implementations + +Return to the same RStudio session, where `project_dir` is already defined: + +```r +source(file.path( + project_dir, + "scripts", + "03_compare_mia_gemelli_corrected.R" +)) +``` + +A successful run ends with `Comparison completed successfully.` and creates +`results/frequency_010/comparison/`. + +## Step 6: render the complete report + +Open `joint-rpca-mia-gemelli-comparison.qmd` in RStudio and select **Render**. +The report reads all validated result tables and figures; it does not repeat the +Joint-RPCA computation. + +Alternatively, in the original PowerShell window: + +```powershell +Set-Location $PROJECT_DIR +quarto render joint-rpca-mia-gemelli-comparison.qmd +``` + +The output is `joint-rpca-mia-gemelli-comparison.html`. Tables and figures are +embedded in this standalone file. + +## Correct execution order + +```text +01 mia/R + → 00 environment and input check + → 02 Gemelli/Python + → 03 comparison/R + → render the QMD +``` + +## Common mistakes + +### Mixing Windows and WSL paths + +Use a Windows path in PowerShell and RStudio. Use `/mnt//...` in +WSL. Do not use a WSL path in PowerShell or a Windows path inside WSL. + +### Opening a new terminal or R session + +PowerShell variables last only for the current PowerShell session. R variables +last only for the current R session. If you open a new session, repeat only the +short path-definition block for that environment. + +### Rendering before comparison + +The QMD stops if required tables or figures are missing. Complete scripts 01, +02, and 03 before rendering. + +### Using the wrong prevalence unit + +Do not change `0.10` to `10` in `config.tsv`. R uses the proportion `0.10`; the +Python runner converts it to Gemelli's `10%` representation. + +### Using nonzero count thresholds with IBDMDB + +The revised IBDMDB inputs are abundance-like and non-integer. This workflow +permits prevalence filtering but stops if nonzero count thresholds are used. +Validate count thresholds separately with an integer-count dataset. diff --git a/vignettes/filterRPCAInput_validation/config.tsv b/vignettes/filterRPCAInput_validation/config.tsv new file mode 100644 index 0000000..f1a299f --- /dev/null +++ b/vignettes/filterRPCAInput_validation/config.tsv @@ -0,0 +1,2 @@ +run_label data_variant min_sample_count min_feature_count min_feature_frequency n_components max_iterations n_test_samples seed expected_mia_version expected_mia_commit +frequency_010 main 0 0 0.10 3 5 10 42 1.21.6 25fa7bb51792566131f2bd8d07b4ebd1fde54b38 diff --git a/vignettes/filterRPCAInput_validation/joint-rpca-mia-gemelli-comparison.html b/vignettes/filterRPCAInput_validation/joint-rpca-mia-gemelli-comparison.html new file mode 100644 index 0000000..f665f57 --- /dev/null +++ b/vignettes/filterRPCAInput_validation/joint-rpca-mia-gemelli-comparison.html @@ -0,0 +1,5077 @@ + + + + + + + + + +Reproducible Joint-RPCA validation + + + + + + + + + + + + + + + + + + + + + + + + +
+ +
+ +
+
+

Reproducible Joint-RPCA validation

+

mia and Gemelli on the revised IBDMDB MGX–MTX data

+
+ + + +
+ + + + +
+ + + +
+ + +
+

Purpose

+

This report checks whether Joint-RPCA in mia and Gemelli produces the same numerical structure when both implementations receive:

+
    +
  • the revised IBDMDB MGX and MTX data;
    • +
  • identical exported input matrices;
    • +
  • equivalent filtering thresholds;
    • +
  • the same training/test split; and
    • +
  • the same number of components and iterations.
    • +
+

The full installation and execution instructions are in SETUP.md. The computational details remain in the numbered scripts; this report contains only the small amount of code needed to display their validated outputs. No machine-specific path is fixed in this report: it uses MIA_GEMELLI_PROJECT when defined, or the folder from which the report is rendered.

+
+
+
+ +
+
+Important +
+
+
+

Filtering and rCLR agreement must be established before differences in Joint-RPCA outputs are interpreted.

+
+
+
+
+

Reproducibility specification

+
+

Run configuration

+
+
+ + ++++++++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Configuration used for this run
run_labeldata_variantmin_sample_countmin_feature_countmin_feature_frequencyn_componentsmax_iterationsn_test_samplesseedexpected_mia_versionexpected_mia_commitgemelli_feature_frequency_percent
frequency_010main000.13510421.21.625fa7bb51792566131f2bd8d07b4ebd1fde54b3810
+
+
+

The configuration stores feature prevalence as an R-style proportion. Thus 0.10 in config.tsv is passed to Gemelli as 10%.

+
+
+

Software versions

+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Pinned mia, R, and platform provenance
itemvalue
mia_version1.21.6
mia_commit25fa7bb51792566131f2bd8d07b4ebd1fde54b38
expected_mia_version1.21.6
expected_mia_commit25fa7bb51792566131f2bd8d07b4ebd1fde54b38
data_variantmain
frequency_unitproportion_0_to_1
r_versionR version 4.5.3 (2026-03-11 ucrt)
platformx86_64-w64-mingw32
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Python and Gemelli provenance
packageversion
python3.8.20
gemelli0.0.12
numpy1.23.5
pandas2.0.3
scikit-bio0.5.9
biom-format2.1.13
gemelli-commitnot-recorded
+
+
+
+
+

Input data

+
+
+ + +++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Selected IBDMDB experiments and alternative experiments
modalityselected_variantselected_assaymain_featuresselected_featuresshared_samplesalternative_experiments
MGXmainmgx17017060kingdom;phylum;class;order;family;genus;strain;original
MTXmainmtx80080060gene_function;gene_taxon
+
+
+ + ++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Dimensions, sparsity, and detected numerical scale
modalityfeaturessampleszero_fractionnon_integer_fractiondetected_scale
MGX170600.68107840.3189216non_integer_abundance
MTX800600.30872920.6912708non_integer_abundance
+
+
+

The revised IBDMDB assays contain non-integer abundance values. Therefore this run evaluates zero/nonzero prevalence filtering. Nonzero count-depth thresholds should be validated separately with an integer-count dataset.

+
+
+
+

Main result

+
+
+ + ++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Stage-wise mia–Gemelli validation summary
stagestatusprimary_resultinterpretation
ConfigurationPASSidentical fields = 11/11All shared filtering, rank, iteration, split-size, and seed settings must match.
FilteringPASSexact retained ID sets = TRUE; minimum Jaccard = 1This must match exactly before interpreting the ordination.
rclr preprocessingPASSminimum missingness agreement = 1; maximum absolute difference = 7.5495e-15Missingness should match exactly and numeric differences should be at floating-point scale.
Sample scoresDESCRIPTIVEminimum aligned Spearman correlation = 1Higher aligned correlations indicate similar component-wise sample placement.
Sample subspaceDESCRIPTIVEProcrustes correlation = 1; relative error = 2.6281e-15High correlation and low relative error indicate similar global sample geometry.
Feature loadingsDESCRIPTIVEminimum aligned Spearman correlation = 1Higher aligned correlations indicate similar modality-specific feature patterns.
Singular valuesDESCRIPTIVEmaximum relative difference = 6.5857e-15Compare corresponding components after applying only the sample-derived sign correction.
Percent varianceDESCRIPTIVEmaximum absolute percentage-point difference = 1.4921e-13Gemelli proportions were converted to percent before comparison.
Sample distancesDESCRIPTIVESpearman correlation = 1; Pearson correlation = 1High correlations indicate similar pairwise sample geometry.
CV trajectoryDESCRIPTIVEmaximum RMSE across mean/std trajectories = 1.0402e-13This is diagnostic and may vary with package and linear-algebra versions.
+
+
+

PASS is used only for exact or tolerance-based validation gates. DESCRIPTIVE stages report agreement metrics without imposing an arbitrary correlation threshold.

+
+
+

Filtering and preprocessing

+
+

Configuration parity

+
+
+ + ++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Configuration recorded independently by mia and Gemelli
fieldmia_valuegemelli_valueidentical
run_labelfrequency_010frequency_010TRUE
data_variantmainmainTRUE
min_sample_count00TRUE
min_feature_count00TRUE
min_feature_frequency0.10.1TRUE
n_components33TRUE
max_iterations55TRUE
n_test_samples1010TRUE
seed4242TRUE
expected_mia_version1.21.61.21.6TRUE
expected_mia_commit25fa7bb51792566131f2bd8d07b4ebd1fde54b3825fa7bb51792566131f2bd8d07b4ebd1fde54b38TRUE
+
+
+
+
+

Filter summaries

+
+
+ + +++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
mia filtering summary
modalityfeatures_beforefeatures_after_first_passfeatures_aftersamples_beforesamples_after_first_passsamples_after
MGX170119119606060
MTX800785785606060
+
+
+ + +++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Gemelli filtering summary
modalityfeatures_beforefeatures_aftersamples_beforesamples_afterr_frequency_proportiongemelli_frequency_percent
MGX17011960600.110
MTX80078560600.110
+
+
+ + +++++++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Agreement of retained feature and sample identifiers
modalitymia_featuresgemelli_featuresfeature_sets_identicalfeature_order_identicalfeature_jaccardmia_samplesgemelli_samplessample_sets_identicalsample_order_identicalsample_jaccard
MGX119119TRUETRUE16060TRUETRUE1
MTX785785TRUETRUE16060TRUETRUE1
+
+
+

The retained identifier sets must agree exactly. Otherwise, the two methods did not analyse the same observations.

+
+
+

Robust CLR

+
+
+ + ++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Agreement of rCLR values and missingness masks
modalitycommon_featurescommon_samplesmissingness_agreementmax_absolute_differencemean_absolute_difference
MGX11960100
MTX78560100
+
+
+

The missingness agreement should equal one. Finite-value differences should be limited to floating-point precision.

+
+
+
+

Joint-RPCA comparison

+
+

Component alignment

+

Components are matched by the permutation that maximizes the absolute sample-score correlations. Signs are determined from sample scores and then applied consistently to feature loadings. Sample scores are additionally scaled for direct comparison of their numerical magnitude.

+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Component permutation, sign, and sample-score scale
mia_componentgemelli_componentsignsample_score_scale_factor
PC1PC111
PC2PC2-11
PC3PC311
+
+
+
+
+

Sample scores and sample subspace

+
+
+ + ++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Component-wise sample-score agreement
mia_componentgemelli_componentsignscale_factorspearmanpearson
PC1PC11111
PC2PC2-1111
PC3PC31111
+
+
+ + ++++++ + + + + + + + + + + + + + + + + +
Orthogonal Procrustes agreement of the sample subspace
common_samplescomponentsrelative_frobenius_errorprocrustes_correlation
60301
+
+
+
+
+
+
+

+
Aligned and scaled sample scores from mia and Gemelli.
+
+
+
+
+
+
+

Feature loadings

+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Feature-loading agreement by modality and component
modalitycomponentfeaturesspearmanpearson
MGXPC111911
MGXPC211911
MGXPC311911
MTXPC178511
MTXPC278511
MTXPC378511
+
+
+
+
+
+
+

+
Feature loadings after applying sample-derived component matching and signs.
+
+
+
+
+
+
+

Singular values and explained variation

+
+
+ + ++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Singular values after component matching
componentgemelli_componentmia_valuegemelli_valueabsolute_differencerelative_difference
PC1PC124.9105224.9105200
PC2PC219.2986019.2986000
PC3PC313.4864313.4864300
+
+
+ + +++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Percentage of variation explained by matched components
componentgemelli_componentmia_percentgemelli_percentabsolute_difference
PC1PC152.8179752.817970
PC2PC231.7006231.700620
PC3PC315.4814115.481410
+
+
+
+
+
+
+

+
Percentage of variation explained by matched components.
+
+
+
+
+
+
+

Pairwise sample distances

+
+
+ + +++++++ + + + + + + + + + + + + + + + + + + +
Agreement of all unique pairwise sample distances
common_samplesspearmanpearsonmean_absolute_differencemax_absolute_difference
601100
+
+
+

Distance agreement is invariant to component signs and therefore provides an important validation of the global sample geometry.

+
+
+

Cross-validation trajectory

+
+
+ + +++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Cross-validation values by iteration
iterationgemelli_mean_CVmia_mean_CVmean_CV_differencegemelli_std_CVmia_std_CVstd_CV_difference
1122.39220122.39220098.11210898.1121080
275.2041775.20417029.47112629.4711260
351.1356951.13569034.98925434.9892540
426.4879026.48790014.49526414.4952640
530.8030530.8030508.9897918.9897910
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + +
Summary agreement of cross-validation trajectories
metriciterationsspearmanpearsonroot_mean_squared_error
mean_CV5110
std_CV5110
+
+
+
+
+
+
+

+
Iteration-wise cross-validation metrics from Gemelli and mia.
+
+
+
+
+
+
+
+

Interpretation

+
    +
  • A filtering mismatch invalidates downstream algorithm attribution.
    • +
  • An rCLR mismatch indicates preprocessing disagreement rather than a Joint-RPCA difference.
    • +
  • Component-wise correlations assess latent ordering, while scale factors quantify normalization differences.
    • +
  • Procrustes and pairwise-distance comparisons assess overall sample geometry.
    • +
  • Cross-validation differences may also reflect numerical-library versions and should be reported rather than hidden by rescaling.
    • +
+
+
+

Detailed output files

+

The report presents every required validation result. Row-level diagnostic files—including aligned sample scores, aligned feature loadings, and pairwise distance values—remain in:

+
results/<run_label>/comparison/
+

They are retained separately because embedding thousands of diagnostic rows would make the report harder to inspect without adding scientific information.

+
+
+

Report environment

+
+
+
R version 4.5.3 (2026-03-11 ucrt)
+Platform: x86_64-w64-mingw32/x64
+Running under: Windows 11 x64 (build 26200)
+
+Matrix products: default
+  LAPACK version 3.12.1
+
+locale:
+[1] LC_COLLATE=English_United States.utf8 
+[2] LC_CTYPE=English_United States.utf8   
+[3] LC_MONETARY=English_United States.utf8
+[4] LC_NUMERIC=C                          
+[5] LC_TIME=English_United States.utf8    
+
+time zone: Europe/Helsinki
+tzcode source: internal
+
+attached base packages:
+[1] stats     graphics  grDevices utils     datasets  methods   base     
+
+loaded via a namespace (and not attached):
+ [1] htmlwidgets_1.6.4 compiler_4.5.3    fastmap_1.2.0     cli_3.6.6        
+ [5] tools_4.5.3       htmltools_0.5.9   otel_0.2.0        rstudioapi_0.19.0
+ [9] yaml_2.3.12       rmarkdown_2.31    knitr_1.51        jsonlite_2.0.0   
+[13] xfun_0.60         digest_0.6.39     rlang_1.3.0       png_0.1-9        
+[17] evaluate_1.0.5
+
+
+
+ +
+ + +
+ + + + + \ No newline at end of file diff --git a/vignettes/filterRPCAInput_validation/joint-rpca-mia-gemelli-comparison.qmd b/vignettes/filterRPCAInput_validation/joint-rpca-mia-gemelli-comparison.qmd new file mode 100644 index 0000000..a3a1938 --- /dev/null +++ b/vignettes/filterRPCAInput_validation/joint-rpca-mia-gemelli-comparison.qmd @@ -0,0 +1,382 @@ +--- +title: "Reproducible Joint-RPCA validation" +subtitle: "mia and Gemelli on the revised IBDMDB MGX–MTX data" +format: + html: + toc: true + toc-depth: 3 + embed-resources: true + code-fold: true + code-summary: "Show report code" + df-print: paged +engine: knitr +execute: + echo: false + warning: false + message: false + error: false +--- + +```{r} +#| label: setup +#| include: false + +project_root <- Sys.getenv("MIA_GEMELLI_PROJECT", unset = "") +if (!nzchar(project_root)) { + project_root <- getwd() +} +project_root <- normalizePath(project_root, winslash = "/", mustWork = TRUE) + +config_path <- file.path(project_root, "config.tsv") +if (!file.exists(config_path)) { + stop( + "config.tsv was not found. Follow SETUP.md and render from the ", + "project root.", + call. = FALSE + ) +} + +config <- read.delim(config_path, check.names = FALSE) +if (nrow(config) != 1L) { + stop("config.tsv must contain exactly one run.", call. = FALSE) +} + +run_label <- as.character(config$run_label[[1L]]) +input_dir <- file.path(project_root, "input", run_label) +mia_dir <- file.path(project_root, "results", run_label, "mia") +gemelli_dir <- file.path(project_root, "results", run_label, "gemelli") +comparison_dir <- file.path( + project_root, "results", run_label, "comparison" +) +figures_dir <- file.path(comparison_dir, "figures") + +read_result <- function(directory, filename) { + path <- file.path(directory, filename) + if (!file.exists(path)) { + stop( + "Missing required result: ", path, + "\nRun the workflow in SETUP.md before rendering.", + call. = FALSE + ) + } + read.delim( + path, + check.names = FALSE, + na.strings = c("NA", "NaN", "") + ) +} + +show_table <- function(x, caption, digits = 6L) { + knitr::kable( + x, + caption = caption, + digits = digits, + row.names = FALSE + ) +} + +show_result <- function(directory, filename, caption, digits = 6L) { + show_table(read_result(directory, filename), caption, digits) +} + +show_figure <- function(filename) { + path <- file.path(figures_dir, filename) + if (!file.exists(path)) { + stop("Missing required figure: ", path, call. = FALSE) + } + knitr::include_graphics(path) +} + +comparison_summary <- read_result( + comparison_dir, + "comparison_summary.tsv" +) +``` + +# Purpose + +This report checks whether Joint-RPCA in `mia` and Gemelli produces the same +numerical structure when both implementations receive: + +- the revised IBDMDB MGX and MTX data; +- identical exported input matrices; +- equivalent filtering thresholds; +- the same training/test split; and +- the same number of components and iterations. + +The full installation and execution instructions are in +[`SETUP.md`](SETUP.md). The computational details remain in the numbered +scripts; this report contains only the small amount of code needed to display +their validated outputs. No machine-specific path is fixed in this report: it +uses `MIA_GEMELLI_PROJECT` when defined, or the folder from which the report is +rendered. + +::: {.callout-important} +Filtering and rCLR agreement must be established before differences in +Joint-RPCA outputs are interpreted. +::: + +# Reproducibility specification + +## Run configuration + +```{r} +config_display <- config +config_display$gemelli_feature_frequency_percent <- + 100 * config_display$min_feature_frequency +show_table(config_display, "Configuration used for this run", digits = 8L) +``` + +The configuration stores feature prevalence as an R-style proportion. Thus +`0.10` in `config.tsv` is passed to Gemelli as `10%`. + +## Software versions + +```{r} +show_result( + mia_dir, + "mia_provenance.tsv", + "Pinned mia, R, and platform provenance" +) +show_result( + gemelli_dir, + "gemelli_versions.tsv", + "Python and Gemelli provenance" +) +``` + +## Input data + +```{r} +show_result( + input_dir, + "dataset_structure.tsv", + "Selected IBDMDB experiments and alternative experiments" +) +show_result( + input_dir, + "input_data_summary.tsv", + "Dimensions, sparsity, and detected numerical scale", + digits = 8L +) +``` + +The revised IBDMDB assays contain non-integer abundance values. Therefore this +run evaluates zero/nonzero prevalence filtering. Nonzero count-depth +thresholds should be validated separately with an integer-count dataset. + +# Main result + +```{r} +show_table( + comparison_summary, + "Stage-wise mia–Gemelli validation summary", + digits = 8L +) +``` + +`PASS` is used only for exact or tolerance-based validation gates. +`DESCRIPTIVE` stages report agreement metrics without imposing an arbitrary +correlation threshold. + +# Filtering and preprocessing + +## Configuration parity + +```{r} +show_result( + comparison_dir, + "config_comparison.tsv", + "Configuration recorded independently by mia and Gemelli" +) +``` + +## Filter summaries + +```{r} +show_result( + mia_dir, + "mia_filter_summary.tsv", + "mia filtering summary" +) +show_result( + gemelli_dir, + "gemelli_filter_summary.tsv", + "Gemelli filtering summary" +) +show_result( + comparison_dir, + "filter_id_comparison.tsv", + "Agreement of retained feature and sample identifiers", + digits = 8L +) +``` + +The retained identifier sets must agree exactly. Otherwise, the two methods +did not analyse the same observations. + +## Robust CLR + +```{r} +show_result( + comparison_dir, + "rclr_comparison.tsv", + "Agreement of rCLR values and missingness masks", + digits = 12L +) +``` + +The missingness agreement should equal one. Finite-value differences should +be limited to floating-point precision. + +# Joint-RPCA comparison + +## Component alignment + +Components are matched by the permutation that maximizes the absolute +sample-score correlations. Signs are determined from sample scores and then +applied consistently to feature loadings. Sample scores are additionally +scaled for direct comparison of their numerical magnitude. + +```{r} +show_result( + comparison_dir, + "component_alignment.tsv", + "Component permutation, sign, and sample-score scale", + digits = 8L +) +``` + +## Sample scores and sample subspace + +```{r} +show_result( + comparison_dir, + "sample_score_correlations.tsv", + "Component-wise sample-score agreement", + digits = 8L +) +show_result( + comparison_dir, + "procrustes_summary.tsv", + "Orthogonal Procrustes agreement of the sample subspace", + digits = 8L +) +``` + +```{r} +#| label: sample-score-figure +#| fig-cap: "Aligned and scaled sample scores from mia and Gemelli." + +show_figure("subject_loadings_python_vs_r.png") +``` + +## Feature loadings + +```{r} +show_result( + comparison_dir, + "feature_loading_correlations.tsv", + "Feature-loading agreement by modality and component", + digits = 8L +) +``` + +```{r} +#| label: feature-loading-figure +#| fig-cap: "Feature loadings after applying sample-derived component matching and signs." + +show_figure("feature_loadings_python_vs_r.png") +``` + +## Singular values and explained variation + +```{r} +show_result( + comparison_dir, + "singular_value_comparison.tsv", + "Singular values after component matching", + digits = 10L +) +show_result( + comparison_dir, + "variance_comparison.tsv", + "Percentage of variation explained by matched components", + digits = 8L +) +``` + +```{r} +#| label: explained-variation-figure +#| fig-cap: "Percentage of variation explained by matched components." + +show_figure("proportion_explained_python_vs_r.png") +``` + +## Pairwise sample distances + +```{r} +show_result( + comparison_dir, + "distance_comparison.tsv", + "Agreement of all unique pairwise sample distances", + digits = 10L +) +``` + +Distance agreement is invariant to component signs and therefore provides an +important validation of the global sample geometry. + +## Cross-validation trajectory + +```{r} +show_result( + comparison_dir, + "cv_iteration_values.tsv", + "Cross-validation values by iteration", + digits = 10L +) +show_result( + comparison_dir, + "cv_comparison.tsv", + "Summary agreement of cross-validation trajectories", + digits = 10L +) +``` + +```{r} +#| label: cross-validation-figure +#| fig-cap: "Iteration-wise cross-validation metrics from Gemelli and mia." + +show_figure("cv_metrics_python_vs_r.png") +``` + +# Interpretation + +- A filtering mismatch invalidates downstream algorithm attribution. +- An rCLR mismatch indicates preprocessing disagreement rather than a + Joint-RPCA difference. +- Component-wise correlations assess latent ordering, while scale factors + quantify normalization differences. +- Procrustes and pairwise-distance comparisons assess overall sample geometry. +- Cross-validation differences may also reflect numerical-library versions and + should be reported rather than hidden by rescaling. + +# Detailed output files + +The report presents every required validation result. Row-level diagnostic +files—including aligned sample scores, aligned feature loadings, and pairwise +distance values—remain in: + +```text +results//comparison/ +``` + +They are retained separately because embedding thousands of diagnostic rows +would make the report harder to inspect without adding scientific information. + +# Report environment + +```{r} +sessionInfo() +``` diff --git a/vignettes/filterRPCAInput_validation/scripts/00_check_gemelli_environment.py b/vignettes/filterRPCAInput_validation/scripts/00_check_gemelli_environment.py new file mode 100644 index 0000000..abdb0a4 --- /dev/null +++ b/vignettes/filterRPCAInput_validation/scripts/00_check_gemelli_environment.py @@ -0,0 +1,110 @@ +#!/usr/bin/env python3 +"""Minimal preflight check for the Gemelli validation environment.""" + +from __future__ import annotations + +import argparse +import csv +import hashlib +import importlib +import importlib.metadata +import inspect +import sys +from pathlib import Path + + +PACKAGES = { + "gemelli": "gemelli", + "biom-format": "biom", + "numpy": "numpy", + "pandas": "pandas", + "scikit-bio": "skbio", + "scipy": "scipy", +} + + +def checksum(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as handle: + for chunk in iter(lambda: handle.read(1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def check_inputs(project_dir: Path) -> None: + """Check the files exported by script 01 without duplicating its logic.""" + + with (project_dir / "config.tsv").open(encoding="utf-8", newline="") as file: + rows = list(csv.DictReader(file, delimiter="\t")) + if len(rows) != 1 or not rows[0].get("run_label"): + raise ValueError("config.tsv must contain one row with a run_label") + + input_dir = project_dir / "input" / rows[0]["run_label"] + manifest = input_dir / "input_manifest.tsv" + with manifest.open(encoding="utf-8", newline="") as file: + entries = list(csv.DictReader(file, delimiter="\t")) + if not entries or set(entries[0]) != {"file", "sha256"}: + raise ValueError("input_manifest.tsv must contain file and sha256 columns") + + for entry in entries: + path = input_dir / entry["file"] + if not path.is_file() or checksum(path) != entry["sha256"]: + raise ValueError(f"Missing or changed R-exported input: {path}") + print(f"OK: verified {len(entries)} R-exported inputs") + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--project-dir", type=Path) + args = parser.parse_args() + + print("Python:", sys.version.split()[0], "at", sys.executable) + modules: dict[str, object] = {} + problems: list[str] = [] + for distribution, module_name in PACKAGES.items(): + try: + modules[module_name] = importlib.import_module(module_name) + version = importlib.metadata.version(distribution) + print(f"OK: {distribution} {version}") + except (ImportError, importlib.metadata.PackageNotFoundError) as exc: + problems.append(f"{distribution}: {exc}") + + if "gemelli" in modules: + from gemelli.rpca import joint_rpca, rpca_table_processing + + required = { + joint_rpca: { + "sample_metadata", + "train_test_column", + "min_sample_count", + "min_feature_count", + "min_feature_frequency", + }, + rpca_table_processing: { + "min_sample_count", + "min_feature_count", + "min_feature_frequency", + }, + } + for function, names in required.items(): + missing = names - set(inspect.signature(function).parameters) + if missing: + problems.append( + f"{function.__name__} is missing: {', '.join(sorted(missing))}" + ) + + if args.project_dir: + try: + check_inputs(args.project_dir.expanduser().resolve()) + except (OSError, ValueError) as exc: + problems.append(str(exc)) + + if problems: + print("\nPreflight failed:", *(f"\n- {item}" for item in problems)) + return 1 + print("\nEnvironment and input check passed.") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/vignettes/filterRPCAInput_validation/scripts/01_run_mia_filter_joint_rpca.R b/vignettes/filterRPCAInput_validation/scripts/01_run_mia_filter_joint_rpca.R new file mode 100644 index 0000000..302297a --- /dev/null +++ b/vignettes/filterRPCAInput_validation/scripts/01_run_mia_filter_joint_rpca.R @@ -0,0 +1,367 @@ +#!/usr/bin/env Rscript + +# Run the pinned mia filterRPCAInput + Joint-RPCA validation and export the +# exact inputs later consumed by Gemelli. + +options(stringsAsFactors = FALSE) + +args <- commandArgs(trailingOnly = TRUE) +project_dir <- Sys.getenv("MIA_GEMELLI_PROJECT", unset = "") +if (!nzchar(project_dir)) project_dir <- if (length(args)) args[[1L]] else getwd() +project_dir <- normalizePath(project_dir, winslash = "/", mustWork = TRUE) + +mia_repo <- Sys.getenv("MIA_REPO", unset = "") +if (!nzchar(mia_repo)) { + stop("Set MIA_REPO to the clean, pinned mia checkout.", call. = FALSE) +} +mia_repo <- normalizePath(mia_repo, winslash = "/", mustWork = TRUE) + +check <- function(ok, message) { + if (!isTRUE(ok)) stop(message, call. = FALSE) +} +write_tsv <- function(x, path) { + write.table(x, path, sep = "\t", quote = FALSE, row.names = FALSE, na = "NA") +} +write_matrix <- function(x, path, id_name) { + x <- as.matrix(x) + out <- data.frame(id = rownames(x), x, check.names = FALSE) + names(out)[[1L]] <- id_name + write_tsv(out, path) +} +write_ids <- function(x, path, id_name) { + out <- data.frame(id = x, check.names = FALSE) + names(out) <- id_name + write_tsv(out, path) +} + +# Configuration --------------------------------------------------------- + +config <- read.delim(file.path(project_dir, "config.tsv"), check.names = FALSE) +required <- c( + "run_label", "data_variant", "min_sample_count", "min_feature_count", + "min_feature_frequency", "n_components", "max_iterations", + "n_test_samples", "seed", "expected_mia_version", "expected_mia_commit" +) +check(nrow(config) == 1L, "config.tsv must contain exactly one row.") +check(!length(setdiff(required, names(config))), "config.tsv is missing columns.") + +value <- function(name, type = c("character", "numeric", "integer")) { + switch(match.arg(type), + character = as.character(config[[name]][[1L]]), + numeric = as.numeric(config[[name]][[1L]]), + integer = as.integer(config[[name]][[1L]]) + ) +} +run_label <- value("run_label") +data_variant <- value("data_variant") +min_sample_count <- value("min_sample_count", "numeric") +min_feature_count <- value("min_feature_count", "numeric") +min_feature_frequency <- value("min_feature_frequency", "numeric") +n_components <- value("n_components", "integer") +max_iterations <- value("max_iterations", "integer") +n_test_samples <- value("n_test_samples", "integer") +seed <- value("seed", "integer") +expected_mia_version <- value("expected_mia_version") +expected_mia_commit <- tolower(value("expected_mia_commit")) + +whole_nonnegative <- function(x) { + length(x) == 1L && is.finite(x) && x >= 0 && x == floor(x) +} +check(nzchar(run_label), "run_label must not be empty.") +check(data_variant %in% c("main", "legacy_original_mgx"), "Unknown data_variant.") +check(whole_nonnegative(min_sample_count), "min_sample_count must be a non-negative integer.") +check(whole_nonnegative(min_feature_count), "min_feature_count must be a non-negative integer.") +check( + length(min_feature_frequency) == 1L && is.finite(min_feature_frequency) && + min_feature_frequency >= 0 && min_feature_frequency <= 1, + "min_feature_frequency must be between 0 and 1." +) +check(!is.na(n_components) && n_components >= 2L, "n_components must be at least 2.") +check(!is.na(max_iterations) && max_iterations %in% 1:100, "max_iterations must be 1--100.") +check(!is.na(n_test_samples) && n_test_samples >= 1L, "n_test_samples must be positive.") +check(!is.na(seed), "seed must be an integer.") +check(grepl("^[0-9a-f]{40}$", expected_mia_commit), "expected_mia_commit must be a full SHA.") + +input_dir <- file.path(project_dir, "input", run_label) +result_dir <- file.path(project_dir, "results", run_label, "mia") +dir.create(input_dir, recursive = TRUE, showWarnings = FALSE) +dir.create(result_dir, recursive = TRUE, showWarnings = FALSE) + +# Pinned mia checkout ---------------------------------------------------- + +packages <- c( + "devtools", "digest", "MultiAssayExperiment", "SingleCellExperiment", + "SummarizedExperiment" +) +missing <- packages[!vapply(packages, requireNamespace, logical(1L), quietly = TRUE)] +check(!length(missing), paste("Install required R packages:", paste(missing, collapse = ", "))) + +git <- function(...) { + output <- suppressWarnings(system2( + "git", c("-C", shQuote(mia_repo), ...), stdout = TRUE, stderr = TRUE + )) + check(is.null(attr(output, "status")) || attr(output, "status") == 0L, + paste(output, collapse = "\n")) + trimws(paste(output, collapse = "\n")) +} +actual_mia_commit <- tolower(git("rev-parse", "HEAD")) +check(identical(actual_mia_commit, expected_mia_commit), paste( + "mia commit mismatch. Expected", expected_mia_commit, "but found", actual_mia_commit +)) +check(!nzchar(git("status", "--porcelain", "--untracked-files=no")), + "The pinned mia checkout has tracked changes.") +devtools::load_all(mia_repo, quiet = TRUE) +actual_mia_version <- as.character(utils::packageVersion("mia")) +check(identical(actual_mia_version, expected_mia_version), paste( + "mia version mismatch. Expected", expected_mia_version, "but found", actual_mia_version +)) +check(all(vapply( + c("filterRPCAInput", "transformAssay", "getJointRPCA"), + exists, logical(1L), envir = asNamespace("mia"), inherits = FALSE +)), "The pinned mia checkout is missing a required function.") + +# Data and common export ------------------------------------------------ + +data_env <- new.env(parent = emptyenv()) +data_path <- file.path(mia_repo, "data", "ibdmdb.rda") +check(file.exists(data_path), paste("Dataset not found:", data_path)) +load(data_path, envir = data_env) +check(exists("ibdmdb", data_env, inherits = FALSE), "Could not load ibdmdb.") +ibdmdb_current <- get("ibdmdb", data_env, inherits = FALSE) +check(methods::is(ibdmdb_current, "MultiAssayExperiment"), "ibdmdb has the wrong class.") +check(all(c("MGX", "MTX") %in% names(ibdmdb_current)), "ibdmdb needs MGX and MTX.") + +mgx <- ibdmdb_current[["MGX"]] +if (data_variant == "legacy_original_mgx") { + check("original" %in% SingleCellExperiment::altExpNames(mgx), + "MGX has no 'original' alternative experiment.") + mgx <- SingleCellExperiment::altExp(mgx, "original") +} +experiments <- list(MGX = mgx, MTX = ibdmdb_current[["MTX"]]) +assay_types <- c(MGX = "mgx", MTX = "mtx") +for (name in names(experiments)) { + check(assay_types[[name]] %in% SummarizedExperiment::assayNames(experiments[[name]]), + paste("Missing", assay_types[[name]], "assay in", name)) +} + +shared_samples <- sort(Reduce(intersect, lapply(experiments, colnames))) +check(length(shared_samples) > n_test_samples, "Too few shared samples for the test split.") + +feature_maps <- Map(function(se, modality) { + se <- se[, shared_samples, drop = FALSE] + original <- rownames(se) + exported <- paste0(modality, "::", make.unique(as.character(original), sep = "_dup")) + rownames(se) <- exported + list( + experiment = se, + map = data.frame( + modality = modality, original_feature_id = original, + exported_feature_id = exported + ) + ) +}, experiments, names(experiments)) +experiments <- setNames(lapply(feature_maps, `[[`, "experiment"), names(experiments)) +feature_id_map <- do.call(rbind, lapply(feature_maps, `[[`, "map")) +check(!anyDuplicated(feature_id_map$exported_feature_id), "Exported feature IDs are not unique.") +write_tsv(feature_id_map, file.path(input_dir, "feature_id_map.tsv")) + +dataset_structure <- do.call(rbind, lapply(names(experiments), function(name) { + source <- ibdmdb_current[[name]] + data.frame( + modality = name, + selected_variant = if (name == "MGX") data_variant else "main", + selected_assay = assay_types[[name]], + main_features = nrow(source), + selected_features = nrow(experiments[[name]]), + shared_samples = ncol(experiments[[name]]), + alternative_experiments = paste(SingleCellExperiment::altExpNames(source), collapse = ";") + ) +})) +write_tsv(dataset_structure, file.path(input_dir, "dataset_structure.tsv")) + +input_matrices <- setNames(Map(function(se, assay) { + SummarizedExperiment::assay(se, assay) +}, experiments, assay_types[names(experiments)]), names(experiments)) + +data_summary <- do.call(rbind, lapply(names(input_matrices), function(name) { + x <- as.matrix(input_matrices[[name]]) + check(all(is.finite(x)) && all(x >= 0), paste(name, "contains invalid values.")) + non_integer <- abs(x - round(x)) > 1e-8 + data.frame( + modality = name, features = nrow(x), samples = ncol(x), + zero_fraction = mean(x == 0), non_integer_fraction = mean(non_integer), + detected_scale = if (any(non_integer)) "non_integer_abundance" else "integer_like" + ) +})) +check( + !any(data_summary$non_integer_fraction > 0) || + (min_sample_count == 0 && min_feature_count == 0), + "IBDMDB is non-integer: use zero count thresholds and a prevalence filter." +) +write_tsv(data_summary, file.path(input_dir, "input_data_summary.tsv")) + +input_paths <- setNames( + file.path(input_dir, paste0("input_", names(input_matrices), ".tsv")), + names(input_matrices) +) +for (name in names(input_matrices)) { + write_matrix(input_matrices[[name]], input_paths[[name]], "feature_id") +} + +# Filtering ------------------------------------------------------------- + +filter_one <- function(se, assay) { + se <- mia::filterRPCAInput( + se, assay.type = assay, min.sample.count = NULL, + min.feature.count = min_feature_count, + min.feature.frequency = min_feature_frequency + ) + mia::filterRPCAInput( + se, assay.type = assay, min.sample.count = min_sample_count, + min.feature.count = NULL, min.feature.frequency = NULL + ) +} +filter_all <- function(x) { + setNames(Map(filter_one, x, assay_types[names(x)]), names(x)) +} +keep_shared <- function(x, stage) { + shared <- sort(Reduce(intersect, lapply(x, colnames))) + check(length(shared) > 0L, paste("No shared samples remain after", stage, "filtering.")) + lapply(x, function(se) se[, shared, drop = FALSE]) +} + +first_pass <- keep_shared(filter_all(experiments), "first-pass") +filtered <- keep_shared(filter_all(first_pass), "second-pass") +shared_filtered <- colnames(filtered[[1L]]) + +filter_summary <- do.call(rbind, lapply(names(experiments), function(name) { + data.frame( + modality = name, + features_before = nrow(experiments[[name]]), + features_after_first_pass = nrow(first_pass[[name]]), + features_after = nrow(filtered[[name]]), + samples_before = ncol(experiments[[name]]), + samples_after_first_pass = ncol(first_pass[[name]]), + samples_after = ncol(filtered[[name]]) + ) +})) +write_tsv(filter_summary, file.path(result_dir, "mia_filter_summary.tsv")) +for (name in names(filtered)) { + write_ids(rownames(filtered[[name]]), + file.path(result_dir, paste0("mia_retained_features_", name, ".tsv")), + "feature_id") + write_ids(colnames(filtered[[name]]), + file.path(result_dir, paste0("mia_retained_samples_", name, ".tsv")), + "sample_id") +} + +# Shared split and checksums -------------------------------------------- + +set.seed(seed) +test_samples <- sort(sample(shared_filtered, n_test_samples)) +split <- data.frame( + sample_id = shared_samples, + train_test = ifelse(shared_samples %in% test_samples, "test", "train") +) +split_path <- file.path(input_dir, "sample_metadata.tsv") +write_tsv(split, split_path) + +manifest_paths <- c( + unname(input_paths), split_path, + file.path(input_dir, c( + "feature_id_map.tsv", "dataset_structure.tsv", "input_data_summary.tsv" + )) +) +manifest <- data.frame( + file = basename(manifest_paths), + sha256 = vapply(manifest_paths, function(path) { + digest::digest(file = path, algo = "sha256", serialize = FALSE) + }, character(1L)) +) +write_tsv(manifest, file.path(input_dir, "input_manifest.tsv")) +write_tsv(config, file.path(result_dir, "mia_config_used.tsv")) + +# rCLR and Joint-RPCA --------------------------------------------------- + +analysis <- filtered +for (name in names(analysis)) { + analysis[[name]] <- mia::transformAssay( + analysis[[name]], assay.type = assay_types[[name]], + method = "rclr", impute = FALSE, name = "rclr" + ) + write_matrix( + SummarizedExperiment::assay(analysis[[name]], "rclr"), + file.path(result_dir, paste0("mia_rclr_", name, ".tsv")), + "feature_id" + ) +} + +set.seed(seed) +result <- mia::getJointRPCA( + MultiAssayExperiment::MultiAssayExperiment(experiments = analysis), + experiments = names(analysis), + assay.types = rep("rclr", length(analysis)), + ncomponents = n_components, + max.iterations = max_iterations, + test.set = test_samples +) + +write_matrix(result, file.path(result_dir, "mia_sample_scores.tsv"), "sample_id") +write_matrix(attr(result, "rotation"), + file.path(result_dir, "mia_feature_loadings.tsv"), "feature_id") + +write_named_vector <- function(x, path, value_name) { + check(!is.null(x), paste("Joint-RPCA result is missing", value_name)) + components <- names(x) + if (is.null(components)) components <- paste0("PC", seq_along(x)) + out <- data.frame(component = components, value = as.numeric(x)) + names(out)[[2L]] <- value_name + write_tsv(out, path) +} +write_named_vector(attr(result, "percentVar"), + file.path(result_dir, "mia_percent_variance.tsv"), "percent") +write_named_vector(attr(result, "varExplained"), + file.path(result_dir, "mia_singular_values.tsv"), "value") + +cv <- as.data.frame(attr(result, "cv_error")) +check(nrow(cv) > 0L, "Joint-RPCA result has no CV trajectory.") +cv$iteration <- seq_len(nrow(cv)) +write_tsv(cv, file.path(result_dir, "mia_cv_error.tsv")) +write_matrix(as.matrix(stats::dist(result)), + file.path(result_dir, "mia_distance_all_samples.tsv"), "sample_id") + +reconstruction <- attr(result, "reconstruct_error") +if (!is.null(reconstruction)) { + write_tsv(data.frame(modality = names(analysis), error = as.numeric(reconstruction)), + file.path(result_dir, "mia_reconstruction_error.tsv")) +} + +# Provenance ------------------------------------------------------------ + +provenance <- data.frame( + item = c( + "mia_version", "mia_commit", "expected_mia_version", + "expected_mia_commit", "data_variant", "frequency_unit", + "r_version", "platform" + ), + value = c( + actual_mia_version, actual_mia_commit, expected_mia_version, + expected_mia_commit, data_variant, "proportion_0_to_1", + R.version.string, R.version$platform + ) +) +write_tsv(provenance, file.path(result_dir, "mia_provenance.tsv")) +saveRDS(list( + config = config, provenance = provenance, data_summary = data_summary, + dataset_structure = dataset_structure, filter_summary = filter_summary, + test_samples = test_samples, + filtered_mae = MultiAssayExperiment::MultiAssayExperiment(experiments = filtered), + result = result +), file.path(result_dir, "mia_result.rds")) +capture.output(sessionInfo(), file = file.path(result_dir, "mia_session_info.txt")) + +cat("\nmia run completed successfully.\n") +cat("Run label:", run_label, "\nPinned mia commit:", actual_mia_commit, "\n") +print(filter_summary, row.names = FALSE) +cat("Test samples:", paste(test_samples, collapse = ", "), "\n") +cat("Results:", result_dir, "\n") diff --git a/vignettes/filterRPCAInput_validation/scripts/02_run_gemelli_filter_joint_rpca.py b/vignettes/filterRPCAInput_validation/scripts/02_run_gemelli_filter_joint_rpca.py new file mode 100644 index 0000000..1f3f7bc --- /dev/null +++ b/vignettes/filterRPCAInput_validation/scripts/02_run_gemelli_filter_joint_rpca.py @@ -0,0 +1,258 @@ +#!/usr/bin/env python3 +"""Run Gemelli on the exact inputs exported by the mia validation.""" + +from __future__ import annotations + +import argparse +import hashlib +import importlib.metadata +import json +import sys +from pathlib import Path + +import numpy as np +import pandas as pd +from biom import Table +from gemelli.preprocessing import matrix_rclr +from gemelli.rpca import joint_rpca, rpca_table_processing + + +CONFIG_TYPES = { + "run_label": str, + "data_variant": str, + "min_sample_count": int, + "min_feature_count": int, + "min_feature_frequency": float, + "n_components": int, + "max_iterations": int, + "n_test_samples": int, + "seed": int, + "expected_mia_version": str, + "expected_mia_commit": str, +} + + +def read_config(path: Path) -> dict[str, object]: + frame = pd.read_csv(path, sep="\t") + missing = set(CONFIG_TYPES) - set(frame.columns) + if len(frame) != 1 or missing: + raise ValueError( + "config.tsv must contain one row and all required columns; missing: " + + ", ".join(sorted(missing)) + ) + row = frame.iloc[0] + config = {name: cast(row[name]) for name, cast in CONFIG_TYPES.items()} + config["expected_mia_commit"] = str(config["expected_mia_commit"]).lower() + + for name in ("min_sample_count", "min_feature_count"): + raw = float(row[name]) + if not np.isfinite(raw) or raw < 0 or not raw.is_integer(): + raise ValueError(f"{name} must be a non-negative integer") + if not 0 <= float(config["min_feature_frequency"]) <= 1: + raise ValueError("min_feature_frequency must be between 0 and 1") + if int(config["n_components"]) < 2: + raise ValueError("n_components must be at least 2") + if not 1 <= int(config["max_iterations"]) <= 100: + raise ValueError("max_iterations must be between 1 and 100") + return config + + +def sha256(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as handle: + for chunk in iter(lambda: handle.read(1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def verify_manifest(input_dir: Path) -> pd.DataFrame: + manifest = pd.read_csv(input_dir / "input_manifest.tsv", sep="\t", dtype=str) + if set(manifest.columns) != {"file", "sha256"}: + raise ValueError("input_manifest.tsv must contain file and sha256 columns") + manifest["observed_sha256"] = [ + sha256(input_dir / filename) for filename in manifest["file"] + ] + manifest["identical"] = manifest["sha256"] == manifest["observed_sha256"] + if not manifest["identical"].all(): + changed = manifest.loc[~manifest["identical"], "file"].tolist() + raise ValueError("R-exported inputs changed: " + ", ".join(changed)) + return manifest + + +def read_table(path: Path) -> tuple[pd.DataFrame, Table]: + frame = pd.read_csv(path, sep="\t", index_col="feature_id") + frame.index, frame.columns = frame.index.astype(str), frame.columns.astype(str) + if frame.empty or frame.index.has_duplicates or frame.columns.has_duplicates: + raise ValueError(f"Empty table or duplicate IDs: {path}") + frame = frame.apply(pd.to_numeric, errors="raise").astype(float) + values = frame.to_numpy() + if not np.isfinite(values).all() or (values < 0).any(): + raise ValueError(f"Values must be finite and non-negative: {path}") + table = Table( + values, + observation_ids=frame.index.tolist(), + sample_ids=frame.columns.tolist(), + ) + return frame, table + + +def read_metadata(path: Path, sample_order: list[str]) -> pd.DataFrame: + metadata = pd.read_csv(path, sep="\t", dtype=str) + if not {"sample_id", "train_test"}.issubset(metadata.columns): + raise ValueError("sample_metadata.tsv needs sample_id and train_test") + if metadata["sample_id"].duplicated().any(): + raise ValueError("sample_metadata.tsv contains duplicate sample IDs") + metadata = metadata.set_index("sample_id") + if metadata.index.tolist() != sample_order: + raise ValueError("Metadata and input tables must have the same sample order") + if set(metadata["train_test"]) != {"train", "test"}: + raise ValueError("train_test must contain both train and test") + return metadata + + +def filter_two_pass(tables: list[Table], **thresholds: float) -> list[Table]: + """Reproduce the two filtering passes inside Gemelli Joint-RPCA.""" + + def process(items: list[Table]) -> list[Table]: + return [rpca_table_processing(table.copy(), **thresholds) for table in items] + + def keep_shared(items: list[Table]) -> list[Table]: + shared = set.intersection(*(set(x.ids(axis="sample")) for x in items)) + if not shared: + raise ValueError("No shared samples remain after filtering") + order = [x for x in items[0].ids(axis="sample") if x in shared] + return [x.filter(order, axis="sample", inplace=False) for x in items] + + return keep_shared(process(keep_shared(process(tables)))) + + +def write_ids(ids: list[str], path: Path, column: str) -> None: + pd.DataFrame({column: ids}).to_csv(path, sep="\t", index=False) + + +def version(name: str) -> str: + try: + return importlib.metadata.version(name) + except importlib.metadata.PackageNotFoundError: + return "not-installed" + + +def installed_commit(name: str) -> str: + try: + text = importlib.metadata.distribution(name).read_text("direct_url.json") + return str(json.loads(text or "{}").get("vcs_info", {}).get("commit_id", "not-recorded")) + except (importlib.metadata.PackageNotFoundError, json.JSONDecodeError): + return "not-recorded" + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--project-dir", required=True, type=Path) + project = parser.parse_args().project_dir.expanduser().resolve() + config = read_config(project / "config.tsv") + run = str(config["run_label"]) + input_dir = project / "input" / run + output = project / "results" / run / "gemelli" + output.mkdir(parents=True, exist_ok=True) + + manifest = verify_manifest(input_dir) + manifest.to_csv( + output / "gemelli_input_manifest_verification.tsv", sep="\t", index=False + ) + + modalities = ["MGX", "MTX"] + inputs = [read_table(input_dir / f"input_{name}.tsv") for name in modalities] + frames, tables = map(list, zip(*inputs)) + sample_order = frames[0].columns.tolist() + if any(frame.columns.tolist() != sample_order for frame in frames[1:]): + raise ValueError("All input tables must have the same sample order") + if set(frames[0].index) & set(frames[1].index): + raise ValueError("Feature IDs must be unique across modalities") + + metadata = read_metadata(input_dir / "sample_metadata.tsv", sample_order) + if (metadata["train_test"] == "test").sum() != int(config["n_test_samples"]): + raise ValueError("The saved test-sample count differs from config.tsv") + + frequency = float(config["min_feature_frequency"]) + thresholds = { + "min_sample_count": int(config["min_sample_count"]), + "min_feature_count": int(config["min_feature_count"]), + "min_feature_frequency": 100 * frequency, + } + filtered = filter_two_pass(tables, **thresholds) + + summaries = [] + for modality, raw, table in zip(modalities, tables, filtered): + features = list(table.ids(axis="observation")) + samples = list(table.ids(axis="sample")) + summaries.append( + { + "modality": modality, + "features_before": raw.shape[0], + "features_after": table.shape[0], + "samples_before": raw.shape[1], + "samples_after": table.shape[1], + "r_frequency_proportion": frequency, + "gemelli_frequency_percent": 100 * frequency, + } + ) + write_ids(features, output / f"gemelli_retained_features_{modality}.tsv", "feature_id") + write_ids(samples, output / f"gemelli_retained_samples_{modality}.tsv", "sample_id") + + rclr = matrix_rclr(table.matrix_data.toarray().T).T + pd.DataFrame(rclr, index=features, columns=samples).rename_axis( + "feature_id" + ).to_csv(output / f"gemelli_rclr_{modality}.tsv", sep="\t", na_rep="NA") + + filter_summary = pd.DataFrame(summaries) + filter_summary.to_csv(output / "gemelli_filter_summary.tsv", sep="\t", index=False) + pd.DataFrame([config]).to_csv(output / "gemelli_config_used.tsv", sep="\t", index=False) + + np.random.seed(int(config["seed"])) + ordination, distance, cv = joint_rpca( + [table.copy() for table in tables], + sample_metadata=metadata, + train_test_column="train_test", + n_components=int(config["n_components"]), + rclr_transform_tables=True, + max_iterations=int(config["max_iterations"]), + **thresholds, + ) + + ordination.samples.rename_axis("sample_id").to_csv( + output / "gemelli_sample_scores.tsv", sep="\t" + ) + ordination.features.rename_axis("feature_id").to_csv( + output / "gemelli_feature_loadings.tsv", sep="\t" + ) + pd.DataFrame( + { + "component": ordination.proportion_explained.index, + "percent": 100 * ordination.proportion_explained.to_numpy(), + } + ).to_csv(output / "gemelli_percent_variance.tsv", sep="\t", index=False) + pd.DataFrame( + {"component": ordination.eigvals.index, "value": ordination.eigvals.to_numpy()} + ).to_csv(output / "gemelli_singular_values.tsv", sep="\t", index=False) + distance.to_data_frame().rename_axis("sample_id").to_csv( + output / "gemelli_distance_all_samples.tsv", sep="\t" + ) + cv.assign(iteration=np.arange(1, len(cv) + 1)).to_csv( + output / "gemelli_cv_error.tsv", sep="\t", index=False + ) + + packages = ["gemelli", "numpy", "pandas", "scikit-bio", "biom-format"] + pd.DataFrame( + { + "package": ["python", *packages, "gemelli-commit"], + "version": [sys.version.split()[0], *(version(x) for x in packages), installed_commit("gemelli")], + } + ).to_csv(output / "gemelli_versions.tsv", sep="\t", index=False) + + print("\nGemelli run completed successfully.") + print(filter_summary.to_string(index=False)) + print("Results:", output) + + +if __name__ == "__main__": + main() diff --git a/vignettes/filterRPCAInput_validation/scripts/03_compare_mia_gemelli_corrected.R b/vignettes/filterRPCAInput_validation/scripts/03_compare_mia_gemelli_corrected.R new file mode 100644 index 0000000..e44fe99 --- /dev/null +++ b/vignettes/filterRPCAInput_validation/scripts/03_compare_mia_gemelli_corrected.R @@ -0,0 +1,535 @@ +#!/usr/bin/env Rscript + +# Compare mia and Gemelli filtering, rCLR, and Joint-RPCA outputs. + +options(stringsAsFactors = FALSE) + +args <- commandArgs(trailingOnly = TRUE) +project_dir <- Sys.getenv("MIA_GEMELLI_PROJECT", unset = "") +if (!nzchar(project_dir)) project_dir <- if (length(args)) args[[1L]] else getwd() +project_dir <- normalizePath(project_dir, winslash = "/", mustWork = TRUE) + +config <- read.delim(file.path(project_dir, "config.tsv"), check.names = FALSE) +if (nrow(config) != 1L) stop("config.tsv must contain one row.", call. = FALSE) +run_label <- as.character(config$run_label[[1L]]) +mia_dir <- file.path(project_dir, "results", run_label, "mia") +gemelli_dir <- file.path(project_dir, "results", run_label, "gemelli") +output_dir <- file.path(project_dir, "results", run_label, "comparison") +figure_dir <- file.path(output_dir, "figures") +dir.create(figure_dir, recursive = TRUE, showWarnings = FALSE) + +read_matrix <- function(path, id_name) { + x <- read.delim(path, check.names = FALSE, na.strings = c("NA", "NaN", "")) + if (!id_name %in% names(x)) stop("Missing ", id_name, " in ", path, call. = FALSE) + ids <- as.character(x[[id_name]]) + x[[id_name]] <- NULL + x <- as.matrix(x) + storage.mode(x) <- "double" + rownames(x) <- ids + x +} +read_ids <- function(path, id_name) { + as.character(read.delim(path, check.names = FALSE)[[id_name]]) +} +write_result <- function(x, name) { + write.table( + x, file.path(output_dir, name), sep = "\t", quote = FALSE, + row.names = FALSE, na = "NA" + ) +} +cor_safe <- function(x, y, method = "spearman") { + suppressWarnings(stats::cor(x, y, method = method, use = "pairwise.complete.obs")) +} +cor_label <- function(x, y) { + keep <- is.finite(x) & is.finite(y) + if (sum(keep) < 3L) return("Spearman rho=NA, p=NA") + test <- suppressWarnings(stats::cor.test( + x[keep], y[keep], method = "spearman", exact = FALSE + )) + paste0( + "Spearman rho=", formatC(unname(test$estimate), digits = 4L), + ", p=", format.pval(test$p.value, digits = 3L, eps = 1e-16) + ) +} +align_rows <- function(x, y) { + ids <- sort(intersect(rownames(x), rownames(y))) + if (!length(ids)) stop("Compared matrices have no common IDs.", call. = FALSE) + list(x = x[ids, , drop = FALSE], y = y[ids, , drop = FALSE]) +} + +# Match components using the sample scores. Sign changes and permutations are +# mathematically equivalent ordination representations, so they are corrected +# before any component-wise comparison. +permutations <- function(x) { + if (length(x) == 1L) return(matrix(x, nrow = 1L)) + do.call(rbind, lapply(seq_along(x), function(i) { + cbind(x[[i]], permutations(x[-i])) + })) +} +align_components <- function(reference, candidate) { + k <- ncol(reference) + if (k != ncol(candidate) || k > 8L) { + stop("Component counts must match and be at most 8.", call. = FALSE) + } + correlations <- abs(stats::cor(reference, candidate, use = "pairwise.complete.obs")) + correlations[!is.finite(correlations)] <- -Inf + orders <- permutations(seq_len(k)) + score <- apply(orders, 1L, function(order) { + sum(correlations[cbind(seq_len(k), order)]) + }) + if (!any(is.finite(score))) stop("Components could not be aligned.", call. = FALSE) + + order <- orders[which.max(score), ] + candidate_names <- colnames(candidate)[order] + candidate <- candidate[, order, drop = FALSE] + signs <- vapply(seq_len(k), function(j) { + value <- cor_safe(reference[, j], candidate[, j], "pearson") + if (!is.finite(value) || value == 0) 1 else sign(value) + }, numeric(1L)) + aligned <- sweep(candidate, 2L, signs, "*") + scale <- vapply(seq_len(k), function(j) { + denominator <- sum(aligned[, j]^2) + if (!is.finite(denominator) || denominator == 0) 1 else + sum(reference[, j] * aligned[, j]) / denominator + }, numeric(1L)) + scaled <- sweep(aligned, 2L, scale, "*") + colnames(aligned) <- colnames(scaled) <- colnames(reference) + list( + aligned = aligned, scaled = scaled, order = order, signs = signs, + scale = scale, mia_names = colnames(reference), + gemelli_names = candidate_names + ) +} + +# Configuration and filtering ------------------------------------------ + +mia_config <- read.delim(file.path(mia_dir, "mia_config_used.tsv"), check.names = FALSE) +gemelli_config <- read.delim( + file.path(gemelli_dir, "gemelli_config_used.tsv"), check.names = FALSE +) +fields <- union(names(mia_config), names(gemelli_config)) +get_config <- function(x, field) { + if (field %in% names(x)) as.character(x[[field]][[1L]]) else NA_character_ +} +config_comparison <- do.call(rbind, lapply(fields, function(field) { + mia_value <- get_config(mia_config, field) + gemelli_value <- get_config(gemelli_config, field) + data.frame( + field = field, mia_value = mia_value, gemelli_value = gemelli_value, + identical = identical(mia_value, gemelli_value) + ) +})) +write_result(config_comparison, "config_comparison.tsv") + +modalities <- c("MGX", "MTX") +jaccard <- function(x, y) length(intersect(x, y)) / length(union(x, y)) +filter_id_comparison <- do.call(rbind, lapply(modalities, function(modality) { + mia_features <- read_ids( + file.path(mia_dir, paste0("mia_retained_features_", modality, ".tsv")), + "feature_id" + ) + gemelli_features <- read_ids( + file.path(gemelli_dir, paste0("gemelli_retained_features_", modality, ".tsv")), + "feature_id" + ) + mia_samples <- read_ids( + file.path(mia_dir, paste0("mia_retained_samples_", modality, ".tsv")), + "sample_id" + ) + gemelli_samples <- read_ids( + file.path(gemelli_dir, paste0("gemelli_retained_samples_", modality, ".tsv")), + "sample_id" + ) + data.frame( + modality = modality, + mia_features = length(mia_features), + gemelli_features = length(gemelli_features), + feature_sets_identical = setequal(mia_features, gemelli_features), + feature_order_identical = identical(mia_features, gemelli_features), + feature_jaccard = jaccard(mia_features, gemelli_features), + mia_samples = length(mia_samples), + gemelli_samples = length(gemelli_samples), + sample_sets_identical = setequal(mia_samples, gemelli_samples), + sample_order_identical = identical(mia_samples, gemelli_samples), + sample_jaccard = jaccard(mia_samples, gemelli_samples) + ) +})) +write_result(filter_id_comparison, "filter_id_comparison.tsv") + +# rCLR ------------------------------------------------------------------ + +rclr_comparison <- do.call(rbind, lapply(modalities, function(modality) { + pair <- align_rows( + read_matrix(file.path(mia_dir, paste0("mia_rclr_", modality, ".tsv")), "feature_id"), + read_matrix(file.path(gemelli_dir, paste0("gemelli_rclr_", modality, ".tsv")), "feature_id") + ) + samples <- sort(intersect(colnames(pair$x), colnames(pair$y))) + x <- pair$x[, samples, drop = FALSE] + y <- pair$y[, samples, drop = FALSE] + finite <- is.finite(x) & is.finite(y) + difference <- abs(x[finite] - y[finite]) + data.frame( + modality = modality, + common_features = nrow(x), + common_samples = ncol(x), + missingness_agreement = mean(is.na(x) == is.na(y)), + max_absolute_difference = if (length(difference)) max(difference) else NA_real_, + mean_absolute_difference = if (length(difference)) mean(difference) else NA_real_ + ) +})) +write_result(rclr_comparison, "rclr_comparison.tsv") + +# Sample scores and subspace -------------------------------------------- + +score_pair <- align_rows( + read_matrix(file.path(mia_dir, "mia_sample_scores.tsv"), "sample_id"), + read_matrix(file.path(gemelli_dir, "gemelli_sample_scores.tsv"), "sample_id") +) +mia_scores <- score_pair$x +gemelli_scores <- score_pair$y +alignment <- align_components(mia_scores, gemelli_scores) +gemelli_scores_aligned <- alignment$aligned +gemelli_scores_scaled <- alignment$scaled + +component_alignment <- data.frame( + mia_component = alignment$mia_names, + gemelli_component = alignment$gemelli_names, + sign = alignment$signs, + sample_score_scale_factor = alignment$scale +) +write_result(component_alignment, "component_alignment.tsv") + +sample_score_correlations <- do.call(rbind, lapply(seq_len(ncol(mia_scores)), function(j) { + data.frame( + mia_component = alignment$mia_names[[j]], + gemelli_component = alignment$gemelli_names[[j]], + sign = alignment$signs[[j]], + scale_factor = alignment$scale[[j]], + spearman = cor_safe(mia_scores[, j], gemelli_scores_aligned[, j]), + pearson = cor_safe(mia_scores[, j], gemelli_scores_aligned[, j], "pearson") + ) +})) +write_result(sample_score_correlations, "sample_score_correlations.tsv") + +aligned_sample_scores <- cbind( + data.frame(sample_id = rownames(mia_scores), check.names = FALSE), + setNames(as.data.frame(mia_scores), paste0("mia_", colnames(mia_scores))), + setNames( + as.data.frame(gemelli_scores_scaled), + paste0("gemelli_aligned_scaled_", colnames(gemelli_scores_scaled)) + ) +) +write_result(aligned_sample_scores, "aligned_sample_scores.tsv") + +mia_centered <- scale(mia_scores, center = TRUE, scale = FALSE) +gemelli_centered <- scale(gemelli_scores, center = TRUE, scale = FALSE) +decomposition <- svd(t(gemelli_centered) %*% mia_centered) +rotation <- decomposition$u %*% t(decomposition$v) +gemelli_procrustes <- gemelli_centered %*% rotation +procrustes_summary <- data.frame( + common_samples = nrow(mia_centered), + components = ncol(mia_centered), + relative_frobenius_error = norm(mia_centered - gemelli_procrustes, "F") / + norm(mia_centered, "F"), + procrustes_correlation = cor_safe( + as.numeric(mia_centered), as.numeric(gemelli_procrustes), "pearson" + ) +) +write_result(procrustes_summary, "procrustes_summary.tsv") + +# Feature loadings ------------------------------------------------------ + +loading_pair <- align_rows( + read_matrix(file.path(mia_dir, "mia_feature_loadings.tsv"), "feature_id"), + read_matrix(file.path(gemelli_dir, "gemelli_feature_loadings.tsv"), "feature_id") +) +mia_loadings <- loading_pair$x +gemelli_loadings <- loading_pair$y[, alignment$order, drop = FALSE] +gemelli_loadings_aligned <- sweep(gemelli_loadings, 2L, alignment$signs, "*") +colnames(gemelli_loadings_aligned) <- colnames(mia_loadings) + +feature_loading_correlations <- do.call(rbind, lapply(modalities, function(modality) { + keep <- startsWith(rownames(mia_loadings), paste0(modality, "::")) + do.call(rbind, lapply(seq_len(ncol(mia_loadings)), function(j) { + data.frame( + modality = modality, + component = colnames(mia_loadings)[[j]], + features = sum(keep), + spearman = cor_safe(mia_loadings[keep, j], gemelli_loadings_aligned[keep, j]), + pearson = cor_safe( + mia_loadings[keep, j], gemelli_loadings_aligned[keep, j], "pearson" + ) + ) + })) +})) +write_result(feature_loading_correlations, "feature_loading_correlations.tsv") + +aligned_feature_loadings <- cbind( + data.frame(feature_id = rownames(mia_loadings), check.names = FALSE), + setNames(as.data.frame(mia_loadings), paste0("mia_", colnames(mia_loadings))), + setNames( + as.data.frame(gemelli_loadings_aligned), + paste0("gemelli_aligned_", colnames(gemelli_loadings_aligned)) + ) +) +write_result(aligned_feature_loadings, "aligned_feature_loadings.tsv") + +# Scalars, distances, and CV -------------------------------------------- + +align_component_table <- function(mia_name, gemelli_name, value_name) { + mia <- read.delim(file.path(mia_dir, mia_name), check.names = FALSE) + gemelli <- read.delim(file.path(gemelli_dir, gemelli_name), check.names = FALSE) + mia <- mia[match(alignment$mia_names, mia$component), , drop = FALSE] + gemelli <- gemelli[match(alignment$gemelli_names, gemelli$component), , drop = FALSE] + if (anyNA(mia$component) || anyNA(gemelli$component)) { + stop("A scalar output is missing an aligned component.", call. = FALSE) + } + list(mia = mia[[value_name]], gemelli = gemelli[[value_name]]) +} + +singular <- align_component_table( + "mia_singular_values.tsv", "gemelli_singular_values.tsv", "value" +) +singular_value_comparison <- data.frame( + component = alignment$mia_names, + gemelli_component = alignment$gemelli_names, + mia_value = singular$mia, + gemelli_value = singular$gemelli, + absolute_difference = abs(singular$mia - singular$gemelli), + relative_difference = abs(singular$mia - singular$gemelli) / + pmax(abs(singular$mia), .Machine$double.eps) +) +write_result(singular_value_comparison, "singular_value_comparison.tsv") + +variance <- align_component_table( + "mia_percent_variance.tsv", "gemelli_percent_variance.tsv", "percent" +) +variance_comparison <- data.frame( + component = alignment$mia_names, + gemelli_component = alignment$gemelli_names, + mia_percent = variance$mia, + gemelli_percent = variance$gemelli, + absolute_difference = abs(variance$mia - variance$gemelli) +) +write_result(variance_comparison, "variance_comparison.tsv") + +distance_pair <- align_rows( + read_matrix(file.path(mia_dir, "mia_distance_all_samples.tsv"), "sample_id"), + read_matrix(file.path(gemelli_dir, "gemelli_distance_all_samples.tsv"), "sample_id") +) +distance_ids <- rownames(distance_pair$x) +mia_distance <- distance_pair$x[, distance_ids, drop = FALSE] +gemelli_distance <- distance_pair$y[, distance_ids, drop = FALSE] +upper <- upper.tri(mia_distance) +indices <- which(upper, arr.ind = TRUE) +distance_pairwise_values <- data.frame( + sample_1 = rownames(mia_distance)[indices[, "row"]], + sample_2 = colnames(mia_distance)[indices[, "col"]], + mia = mia_distance[upper], + gemelli = gemelli_distance[upper] +) +write_result(distance_pairwise_values, "distance_pairwise_values.tsv") +distance_comparison <- data.frame( + common_samples = length(distance_ids), + spearman = cor_safe(mia_distance[upper], gemelli_distance[upper]), + pearson = cor_safe(mia_distance[upper], gemelli_distance[upper], "pearson"), + mean_absolute_difference = mean(abs(mia_distance[upper] - gemelli_distance[upper])), + max_absolute_difference = max(abs(mia_distance[upper] - gemelli_distance[upper])) +) +write_result(distance_comparison, "distance_comparison.tsv") + +mia_cv <- read.delim(file.path(mia_dir, "mia_cv_error.tsv"), check.names = FALSE) +gemelli_cv <- read.delim(file.path(gemelli_dir, "gemelli_cv_error.tsv"), check.names = FALSE) +iterations <- intersect(mia_cv$iteration, gemelli_cv$iteration) +mia_cv <- mia_cv[match(iterations, mia_cv$iteration), , drop = FALSE] +gemelli_cv <- gemelli_cv[match(iterations, gemelli_cv$iteration), , drop = FALSE] +cv_iteration_values <- data.frame( + iteration = iterations, + gemelli_mean_CV = gemelli_cv$mean_CV, + mia_mean_CV = mia_cv$mean_CV, + mean_CV_difference = mia_cv$mean_CV - gemelli_cv$mean_CV, + gemelli_std_CV = gemelli_cv$std_CV, + mia_std_CV = mia_cv$std_CV, + std_CV_difference = mia_cv$std_CV - gemelli_cv$std_CV +) +write_result(cv_iteration_values, "cv_iteration_values.tsv") + +cv_comparison <- do.call(rbind, lapply(c("mean_CV", "std_CV"), function(metric) { + difference <- mia_cv[[metric]] - gemelli_cv[[metric]] + data.frame( + metric = metric, iterations = length(iterations), + spearman = cor_safe(mia_cv[[metric]], gemelli_cv[[metric]]), + pearson = cor_safe(mia_cv[[metric]], gemelli_cv[[metric]], "pearson"), + root_mean_squared_error = sqrt(mean(difference^2)) + ) +})) +write_result(cv_comparison, "cv_comparison.tsv") + +# Summary --------------------------------------------------------------- + +configuration_exact <- isTRUE(all(config_comparison$identical)) +filtering_exact <- isTRUE(all( + filter_id_comparison$feature_sets_identical & + filter_id_comparison$sample_sets_identical & + filter_id_comparison$feature_jaccard == 1 & + filter_id_comparison$sample_jaccard == 1 +)) +rclr_exact <- isTRUE(all( + rclr_comparison$missingness_agreement == 1 & + rclr_comparison$max_absolute_difference <= 1e-10 +)) + +comparison_summary <- data.frame( + stage = c( + "Configuration", "Filtering", "rclr preprocessing", "Sample scores", + "Sample subspace", "Feature loadings", "Singular values", + "Percent variance", "Sample distances", "CV trajectory" + ), + status = c( + if (configuration_exact) "PASS" else "CHECK", + if (filtering_exact) "PASS" else "CHECK", + if (rclr_exact) "PASS" else "CHECK", + rep("DESCRIPTIVE", 7L) + ), + primary_result = c( + sprintf("identical fields = %s/%s", sum(config_comparison$identical), nrow(config_comparison)), + sprintf( + "exact retained ID sets = %s; minimum Jaccard = %s", + filtering_exact, + signif(min(filter_id_comparison$feature_jaccard, filter_id_comparison$sample_jaccard), 5L) + ), + sprintf( + "minimum missingness agreement = %s; maximum absolute difference = %s", + signif(min(rclr_comparison$missingness_agreement), 5L), + signif(max(rclr_comparison$max_absolute_difference), 5L) + ), + paste("minimum aligned Spearman correlation =", signif(min(sample_score_correlations$spearman), 5L)), + sprintf( + "Procrustes correlation = %s; relative error = %s", + signif(procrustes_summary$procrustes_correlation, 5L), + signif(procrustes_summary$relative_frobenius_error, 5L) + ), + paste("minimum aligned Spearman correlation =", signif(min(feature_loading_correlations$spearman), 5L)), + paste("maximum relative difference =", signif(max(singular_value_comparison$relative_difference), 5L)), + paste("maximum absolute percentage-point difference =", signif(max(variance_comparison$absolute_difference), 5L)), + sprintf( + "Spearman correlation = %s; Pearson correlation = %s", + signif(distance_comparison$spearman, 5L), signif(distance_comparison$pearson, 5L) + ), + paste("maximum RMSE across mean/std trajectories =", signif(max(cv_comparison$root_mean_squared_error), 5L)) + ), + interpretation = c( + "All shared filtering, rank, iteration, split-size, and seed settings must match.", + "This must match exactly before interpreting the ordination.", + "Missingness should match exactly and numeric differences should be at floating-point scale.", + "Higher aligned correlations indicate similar component-wise sample placement.", + "High correlation and low relative error indicate similar global sample geometry.", + "Higher aligned correlations indicate similar modality-specific feature patterns.", + "Compare corresponding components after applying only the sample-derived sign correction.", + "Gemelli proportions were converted to percent before comparison.", + "High correlations indicate similar pairwise sample geometry.", + "This is diagnostic and may vary with package and linear-algebra versions." + ) +) +write_result(comparison_summary, "comparison_summary.tsv") + +# Four report figures --------------------------------------------------- + +axis_limits <- function(x, y) { + values <- c(x, y) + values <- values[is.finite(values)] + if (!length(values)) stop("A plot has no finite values.", call. = FALSE) + limits <- range(values) + if (diff(limits) == 0) { + padding <- max(abs(limits[[1L]]) * 0.05, 1e-8) + limits <- limits + c(-padding, padding) + } + limits +} +scatter <- function( + x, y, xlab, ylab, colour = "#2878B5", cex = 1, + title_prefix = "" +) { + limits <- axis_limits(x, y) + plot( + x, y, xlab = xlab, ylab = ylab, + main = paste0(title_prefix, cor_label(x, y)), + pch = 19, col = colour, cex = cex, xlim = limits, ylim = limits + ) + abline(0, 1, col = "#CC3311", lwd = 1.5, lty = 2) +} +draw_subjects <- function() { + k <- ncol(mia_scores) + par(mfrow = c(ceiling(k / min(3L, k)), min(3L, k)), + mar = c(4.5, 4.5, 3.5, 1), oma = c(0, 0, 3, 0), pty = "s") + for (j in seq_len(k)) scatter( + gemelli_scores_scaled[, j], mia_scores[, j], + paste("Gemelli / Python aligned and scaled", colnames(mia_scores)[[j]]), + paste("Mia / R", colnames(mia_scores)[[j]]) + ) + mtext("Subject/sample loadings: Gemelli (Python) versus current Mia (R)", + side = 3, outer = TRUE, font = 2, cex = 1.2) +} +draw_features <- function() { + k <- ncol(mia_loadings) + par(mfrow = c(length(modalities), k), mar = c(4.2, 4.2, 3.5, 1), + oma = c(0, 0, 3, 0), pty = "s") + for (modality in modalities) { + keep <- startsWith(rownames(mia_loadings), paste0(modality, "::")) + for (j in seq_len(k)) scatter( + gemelli_loadings_aligned[keep, j], mia_loadings[keep, j], + paste("Gemelli / Python", colnames(mia_loadings)[[j]]), + paste("Mia / R", colnames(mia_loadings)[[j]]), + if (modality == "MGX") "#AA3377" else "#EE7733", 0.55, + paste0(modality, ": ") + ) + } + mtext("Feature loadings by modality: Gemelli (Python) versus current Mia (R)", + side = 3, outer = TRUE, font = 2, cex = 1.2) +} +draw_variance <- function() { + par(mfrow = c(1, 1), mar = c(5, 5, 4, 1), oma = c(0, 0, 0, 0), pty = "m") + maximum <- max(variance_comparison[c("gemelli_percent", "mia_percent")], na.rm = TRUE) + barplot( + rbind(variance_comparison$gemelli_percent, variance_comparison$mia_percent), + beside = TRUE, names.arg = variance_comparison$component, + col = c("#EE7733", "#3366AA"), xlab = "Principal component", + ylab = "Percent variance explained", + main = "Proportion explained: Gemelli (Python) versus Mia (R)", + ylim = c(0, maximum * 1.15) + ) + legend("topright", c("Gemelli / Python", "Mia / R"), + fill = c("#EE7733", "#3366AA"), bty = "n") +} +draw_cv <- function() { + par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3.5, 1), + oma = c(0, 0, 3, 0), pty = "s") + for (metric in c("mean_CV", "std_CV")) scatter( + gemelli_cv[[metric]], mia_cv[[metric]], + paste("Gemelli / Python", metric), paste("Mia / R", metric) + ) + mtext("Cross-validation metrics by iteration", side = 3, + outer = TRUE, font = 2, cex = 1.2) +} +save_png <- function(name, width, height, draw) { + png(file.path(figure_dir, name), width = width, height = height, res = 160) + tryCatch(draw(), finally = dev.off()) +} +save_png("subject_loadings_python_vs_r.png", 1900, 700, draw_subjects) +save_png("feature_loadings_python_vs_r.png", 1900, 1250, draw_features) +save_png("proportion_explained_python_vs_r.png", 1100, 800, draw_variance) +save_png("cv_metrics_python_vs_r.png", 1500, 700, draw_cv) + +comparison_pdf <- file.path(output_dir, "mia_gemelli_comparison_plots.pdf") +pdf(comparison_pdf, width = 11, height = 8.5, onefile = TRUE) +tryCatch({ + draw_subjects() + draw_features() + draw_variance() + draw_cv() +}, finally = dev.off()) + +cat("\nComparison completed successfully.\n") +print(comparison_summary, row.names = FALSE) +cat("\nFour-plot comparison PDF:", comparison_pdf, "\n") +cat("Outputs:", output_dir, "\n")