Add Logistic PCA (LPCA) analysis

.github/workflows/build-containers.yml

@@ -73,3 +73,8 @@ jobs:
       # of detected changes.
       always_build: true
       context: ./
+  build-and-remove-lpca:
+    uses: "./.github/workflows/build-and-remove-template.yml"
+    with:
+      path: docker-wrappers/LPCA
+      container: reedcompbio/lpca

Snakefile

@@ -91,6 +91,16 @@ def make_final_input(wildcards):
         final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}ensemble-pathway.txt',out_dir=out_dir,sep=SEP,dataset=dataset_labels,algorithm_params=algorithms_with_params))
         final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}jaccard-matrix.txt',out_dir=out_dir,sep=SEP,dataset=dataset_labels,algorithm_params=algorithms_with_params))
         final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}jaccard-heatmap.png',out_dir=out_dir,sep=SEP,dataset=dataset_labels,algorithm_params=algorithms_with_params))
+
+    if _config.config.analysis_include_lpca:
+        final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}{algorithm}-lpca-scores.csv',out_dir=out_dir, sep=SEP, dataset=dataset_labels, algorithm=algorithms_mult_param_combos))
+        final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}{algorithm}-lpca.png',out_dir=out_dir, sep=SEP,dataset=dataset_labels,algorithm=algorithms_mult_param_combos))
+        final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}{algorithm}-lpca-coordinates.txt',out_dir=out_dir, sep=SEP, dataset=dataset_labels,algorithm=algorithms_mult_param_combos))
+        final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}{algorithm}-lpca-binary-matrix.csv',out_dir=out_dir, sep=SEP, dataset=dataset_labels,algorithm=algorithms_mult_param_combos))
+        final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}lpca.png',out_dir=out_dir, sep=SEP, dataset=dataset_labels))
+        final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}lpca-scores.csv',out_dir=out_dir, sep=SEP, dataset=dataset_labels))
+        final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}lpca-coordinates.txt',out_dir=out_dir, sep=SEP, dataset=dataset_labels))
+        final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}lpca-binary-matrix.csv',out_dir=out_dir, sep=SEP, dataset=dataset_labels))
 
     if _config.config.analysis_include_ml_aggregate_algo:
         final_input.extend(expand('{out_dir}{sep}{dataset}-ml{sep}{algorithm}-pca.png',out_dir=out_dir,sep=SEP,dataset=dataset_labels,algorithm=algorithms_mult_param_combos))
@@ -356,6 +366,7 @@ rule ml_analysis:
         ml.hac_horizontal(summary_df, output.hac_image_horizontal, output.hac_clusters_horizontal, **hac_params)
         ml.pca(summary_df, output.pca_image, output.pca_variance, output.pca_coordinates, **pca_params)
 
+
 # Calculated Jaccard similarity between output pathways for each dataset
 rule jaccard_similarity:
     input:
@@ -403,6 +414,32 @@ rule ml_analysis_aggregate_algo:
         ml.hac_horizontal(summary_df, output.hac_image_horizontal, output.hac_clusters_horizontal, **hac_params)
         ml.pca(summary_df, output.pca_image, output.pca_variance, output.pca_coordinates, **pca_params)
 
+rule lpca_analysis_all:
+    input:
+        pathways = expand('{out_dir}{sep}{{dataset}}-{algorithm_params}{sep}pathway.txt', out_dir=out_dir, sep=SEP, algorithm_params=algorithms_with_params)
+    output:
+        lpca_scores = SEP.join([out_dir, '{dataset}-ml', 'lpca-scores.csv']),
+        lpca_png = SEP.join([out_dir, '{dataset}-ml', 'lpca.png']),
+        lpca_coord = SEP.join([out_dir, '{dataset}-ml', 'lpca-coordinates.txt']),
+        lpca_matrix = SEP.join([out_dir, '{dataset}-ml', 'lpca-binary-matrix.csv'])
+    run:
+        from spras.analysis import lpca
+        summary_df = ml.summarize_networks(input.pathways)
+        lpca.run_lpca(
+            summary_df,
+            output.lpca_scores,
+            output.lpca_matrix,
+            k=_config.config.lpca_params.k,
+            m=_config.config.lpca_params.m,
+            cv=_config.config.lpca_params.cv,
+            container_settings=container_settings
+        )
+        lpca.plot_lpca(
+            output.lpca_scores,
+            output.lpca_png,
+            output.lpca_coord
+        )
+
 # Ensemble the output pathways for each dataset per algorithm
 rule ensemble_per_algo:
     input:

config/config.yaml

@@ -272,3 +272,15 @@ analysis:
     # adds evaluation per algorithm per dataset-goldstandard pair
     # evaluation per algorithm will not run unless ml include and ml aggregate_per_algorithm are set to true
     aggregate_per_algorithm: true
+  lpca:
+    # if true, runs LPCA in addition to the existing PCA (both analyses will run in parallel).
+    # Running both helps compare the two methods on the same data.
+    include: false
+    # number of principal components to compute
+    k: 2
+    # fixed value of the logisticPCA tuning parameter m, used when cv is false.
+    # default of 6 was selected based on cross-validation experiments across multiple
+    # SPRAS algorithms (see https://github.com/Jeebjean/lpca-spras for details)
+    m: 6
+    # if true, choose m by cross-validation; if false, use the fixed m above
+    cv: false

docker-wrappers/LPCA/Dockerfile

@@ -0,0 +1,31 @@
+# Logistic PCA (logisticPCA) wrapper for SPRAS.
+# Invoked by spras/analysis/lpca.py, which supplies the full command
+# (Rscript /app/run_lpca.R ... or /app/run_cv.R ...), so no ENTRYPOINT is set.
+
+# Pinned R version for reproducibility. Bump deliberately, not to :latest.
+FROM rocker/r-base:4.4.2
+
+LABEL org.opencontainers.image.source="https://github.com/Reed-CompBio/spras"
+LABEL org.opencontainers.image.description="Logistic PCA (logisticPCA) wrapper for SPRAS"
+
+# System libraries needed to compile ggplot2 (a hard Import of logisticPCA)
+# and its dependency stack from source on Debian.
+RUN apt-get update && apt-get install -y --no-install-recommends \
+        libcurl4-openssl-dev \
+        libssl-dev \
+        libxml2-dev \
+        libfontconfig1-dev \
+        libfreetype6-dev \
+        libpng-dev \
+        libtiff5-dev \
+        libjpeg-dev \
+    && rm -rf /var/lib/apt/lists/*
+
+# logisticPCA is still on CRAN (last published 2016) and pulls in ggplot2.
+RUN Rscript -e "install.packages(c('logisticPCA', 'rARPACK'), repos='https://cran.r-project.org')" \
+    && Rscript -e "library(logisticPCA); library(rARPACK)"
+
+COPY run_lpca.R /app/run_lpca.R
+COPY run_cv.R /app/run_cv.R
+
+WORKDIR /app

docker-wrappers/LPCA/README.md

@@ -0,0 +1,57 @@
+# LPCA (Logistic PCA) wrapper
+
+Docker image: https://hub.docker.com/r/reedcompbio/lpca
+
+This wrapper runs [logisticPCA](https://github.com/andland/logisticPCA)
+
+This wrapper runs [logisticPCA](https://github.com/andland/logisticPCA)
+([Landgraf & Lee, 2020](https://doi.org/10.1016/j.jmva.2020.104668)) as a SPRAS analysis step. It reduces the binary
+edge-by-run matrix built from a set of pathway reconstruction outputs to a small
+number of components and reports the proportion of deviance explained.
+
+The analysis is driven by the `analysis.lpca` config block and the
+`lpca_analysis` Snakemake rule, and is implemented in `spras/analysis/lpca.py`.
+
+## Configuration
+
+    analysis:
+      lpca:
+        include: false   # run the LPCA analysis per algorithm
+        k: 2             # number of principal components
+        m: 6             # fixed logisticPCA tuning parameter, used when cv is false
+        cv: false        # true: choose m by cross-validation; false: use the fixed m
+
+LPCA only runs for algorithms with multiple parameter combinations, so that the
+binary matrix has more than one column. It also needs a reasonable number of
+observations to be meaningful; very small inputs (such as the bundled example
+datasets) produce degenerate results, which is why it is disabled by default.
+
+## Scripts
+
+The image contains two R scripts under `/app`:
+
+- `run_lpca.R <input> <output> <k> <m>`: runs logisticPCA with a fixed `m` and
+  writes the scores CSV plus a sibling `<output basename>_deviance.txt`.
+- `run_cv.R <input> <output> <k>`: cross-validates `m` over 1..20 and writes a
+  CSV with a `best_m` column (plus a `_curve.csv` with the full CV curve). Only
+  used when `cv: true`.
+
+Both read a CSV whose first column holds row labels and whose remaining columns
+are binary (0/1) features, and coerce missing values to 0.
+
+## Dependency note
+
+`logisticPCA` declares `ggplot2` as a hard `Imports` dependency, so building the
+image compiles the ggplot2 stack. The Dockerfile installs the required Debian
+system libraries for that. To build from a source CRAN mirror the `repos`
+argument already points at `https://cran.r-project.org`.
+
+## Building and publishing the image
+
+For the SPRAS default registry to resolve the image, it must be published as
+`docker.io/reedcompbio/lpca:v1`, which requires access to the `reedcompbio`
+Docker Hub organization:
+
+    docker build -t reedcompbio/lpca:v1 docker-wrappers/lpca/
+    docker push reedcompbio/lpca:v1
+

docker-wrappers/LPCA/run_cv.R

@@ -0,0 +1,36 @@
+# run_cv.R
+# Finds the optimal m for a given k using cross-validation
+
+args = commandArgs(trailingOnly = TRUE)
+input_file = args[1]
+output_file = args[2]
+k = as.integer(args[3])
+
+set.seed(42)
+
+# Load data
+library(logisticPCA)
+data = read.csv(input_file, row.names = NULL)
+data = data[, -1]
+data_matrix = as.matrix(data)
+data_matrix[is.na(data_matrix)] = 0
+
+# Cross-validation over m, fixed k
+cv_result = cv.lpca(data_matrix, ks = k, ms = 1:20)
+best_m = which.min(cv_result)
+
+cat("Cross-validation done for k =", k, "\n")
+cat("Best m:", best_m, "\n")
+
+# Save best m
+write.csv(data.frame(k = k, best_m = best_m), output_file, row.names = FALSE)
+
+# Save full CV curve (all m values and their reconstruction error)
+cv_curve_file = sub("\\.csv$", "_curve.csv", output_file)
+cv_df = data.frame(
+  m = 1:20,
+  reconstruction_error = as.numeric(cv_result),
+  is_best = (1:20) == best_m
+)
+write.csv(cv_df, cv_curve_file, row.names = FALSE)
+cat("CV curve saved to", cv_curve_file, "\n")

docker-wrappers/LPCA/run_lpca.R

@@ -0,0 +1,30 @@
+# Read command line arguments
+args = commandArgs(trailingOnly = TRUE)
+input_file = args[1]
+output_file = args[2]
+k = as.integer(args[3])
+m = as.numeric(args[4])
+
+# Load data
+library(logisticPCA)
+data = read.csv(input_file, row.names = NULL)
+row_labels = data[, 1]
+data = data[, -1]
+data_matrix = as.matrix(data)
+data_matrix[is.na(data_matrix)] = 0
+
+# Run LPCA
+model = logisticPCA(data_matrix, k = k, m = m, partial_decomp = TRUE)
+
+# Save scores
+scores = model$PCs
+rownames(scores) = row_labels
+write.csv(scores, output_file, row.names = TRUE)
+
+# Save deviance explained
+deviance_file = sub("\\.csv$", "_deviance.txt", output_file)
+writeLines(as.character(model$prop_deviance_expl), deviance_file)
+
+cat("LPCA done! Scores saved to", output_file, "\n")
+cat("Score dimensions:", nrow(scores), "x", ncol(scores), "\n")
+cat("Proportion of deviance explained:", model$prop_deviance_expl, "\n")

spras/analysis/lpca.py

@@ -0,0 +1,141 @@
+# Logistic PCA analysis for SPRAS
+# Runs LPCA on the binary edge x algorithm matrix produced by summarize_networks.
+# Configured through the analysis.lpca block (k, m, cv, transpose).
+
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import pandas as pd
+import seaborn as sns
+
+from spras.analysis.ml import create_palette
+from spras.config.container_schema import ProcessedContainerSettings
+from spras.containers import prepare_volume, run_container_and_log
+
+# Published as docker.io/reedcompbio/lpca:v1. Only the suffix is given here;
+# the registry prefix is resolved from the container settings.
+LPCA_CONTAINER_SUFFIX = 'lpca:v1'
+LPCA_WORK_DIR = '/app'
+
+
+def run_lpca(
+    dataframe: pd.DataFrame,
+    output_scores: str,
+    output_matrix: str,
+    k: int = 2,
+    m: float = 6,
+    cv: bool = False,
+    container_settings=None
+) -> None:
+    """
+    Runs Logistic PCA on the binary edge x algorithm matrix built from SPRAS
+    algorithm output files.
+
+    @param dataframe: binary dataframe of edge comparison between algorithms from summarize_networks
+    @param output_matrix: path to write the binary matrix CSV (used as Docker volume mount)
+    @param output_scores: path to write the LPCA PC scores CSV
+    @param k: number of principal components (default 2)
+    @param m: fixed logisticPCA tuning parameter, used when cv is False
+    @param cv: if True, determine m by cross-validation; if False, use the
+        fixed m directly (default False)
+    @param container_settings: configure the container runtime (Docker or Singularity)
+
+    Note: KDE-based parameter selection (used by PCA) always uses PCA scores, even when LPCA
+    is also enabled. KDE integration with LPCA may be added in a future update.
+    """
+    if not container_settings:
+        container_settings = ProcessedContainerSettings()
+
+    # Step 1: build the binary edge x algorithm matrix
+    matrix = dataframe
+    # Transpose so runs are the observations, mirroring the classic PCA analysis
+    matrix = matrix.T
+    print(f'LPCA: Matrix shape: {matrix.shape}')
+
+    # Step 2: write the binary matrix for Docker volume mount
+    output_dir = Path(output_scores).parent
+    output_dir.mkdir(parents=True, exist_ok=True)
+    Path(output_matrix).parent.mkdir(parents=True, exist_ok=True)
+    matrix_path = output_matrix
+    matrix.to_csv(matrix_path)
+
+    # Step 3: mount the matrix and the scores output
+    volumes = []
+    bind_path, mapped_matrix = prepare_volume(matrix_path, LPCA_WORK_DIR, container_settings)
+    volumes.append(bind_path)
+    bind_path, mapped_scores = prepare_volume(output_scores, LPCA_WORK_DIR, container_settings)
+    volumes.append(bind_path)
+
+    # Step 4: choose m, optionally via cross-validation
+    if cv:
+        algo_name = Path(output_scores).name.replace('-lpca-scores.csv', '')
+        cv_output_path = str(output_dir / f'{algo_name}-lpca_cv_result.csv')
+        bind_path, mapped_cv_output = prepare_volume(cv_output_path, LPCA_WORK_DIR, container_settings)
+        volumes.append(bind_path)
+
+        print(f'LPCA: Running cross-validation with k={k}...')
+        command_cv = ['Rscript', '/app/run_cv.R', mapped_matrix, mapped_cv_output, str(k)]
+        run_container_and_log('LPCA-CV', LPCA_CONTAINER_SUFFIX, command_cv, volumes,
+                              LPCA_WORK_DIR, None, container_settings)
+
+        if not Path(cv_output_path).exists():
+            raise FileNotFoundError(
+                f'LPCA: Cross-validation output not found at {cv_output_path}. '
+                'Check the LPCA Docker container logs for errors.'
+            )
+        m_used = pd.read_csv(cv_output_path)['best_m'][0]
+        print(f'LPCA: Best m found by CV: {m_used}')
+    else:
+        m_used = m
+        print(f'LPCA: Using fixed m={m_used}')
+
+    # Step 5: run LPCA with the chosen k and m
+    print(f'LPCA: Running LPCA with k={k}, m={m_used}...')
+    command_lpca = ['Rscript', '/app/run_lpca.R', mapped_matrix, mapped_scores, str(k), str(m_used)]
+    run_container_and_log('LPCA', LPCA_CONTAINER_SUFFIX, command_lpca, volumes,
+                          LPCA_WORK_DIR, None, container_settings)
+
+    print(f'LPCA: Done! Scores saved to {output_scores}')
+
+def plot_lpca(scores_file: str, output_png: str, output_coord: str, labels: bool = True) -> None:
+    """
+    Creates a scatterplot of the first two LPCA principal components.
+    @param scores_file: path to the LPCA scores CSV file
+    @param output_png: path to save the scatterplot PNG
+    @param output_coord: path to save the PC coordinates
+    @param labels: if True, adds algorithm labels to the plot
+    """
+    scores = pd.read_csv(scores_file, index_col=0)
+
+    if scores.empty:
+        print('LPCA: Scores file is empty, skipping plot.')
+        return
+
+    # Extract algorithm names from the index
+    column_names = [idx.split('-')[-3] if '-' in idx else idx for idx in scores.index]
+
+    X = scores.values
+    fig, ax = plt.subplots(figsize=(10, 8))
+
+    label_color_map = create_palette(column_names)
+    sns.scatterplot(x=X[:, 0], y=X[:, 1], hue=column_names, palette=label_color_map, s=70, ax=ax)
+
+    if labels:
+        for i, label in enumerate(scores.index):
+            ax.annotate(label, (X[i, 0], X[i, 1]), fontsize=6, alpha=0.7)
+
+    ax.set_xlabel('PC1')
+    ax.set_ylabel('PC2')
+    ax.set_title('Logistic PCA')
+
+    plt.tight_layout()
+
+    # Save PNG
+    Path(output_png).parent.mkdir(parents=True, exist_ok=True)
+    plt.savefig(output_png, dpi=200)
+    plt.close()
+
+    # Save coordinates
+    coord_df = pd.DataFrame(X, columns=['PC1', 'PC2'], index=scores.index)
+    coord_df.to_csv(output_coord)
+    print(f'LPCA: Plot saved to {output_png}')

spras/analysis/ml.py

@@ -156,7 +156,8 @@ def pca(dataframe: pd.DataFrame, output_png: str | PathLike, output_var: str | P
 
     # center binary data by subtracting the column-wise mean
     # allows PCA to focus on edge inclusion patterns across runs rather than raw output volume.
-    # TODO: replace PCA https://github.com/Reed-CompBio/spras/issues/271
+    # TODO: consider replacing PCA with LPCA for binary data https://github.com/Reed-CompBio/spras/issues/271
+    # LPCA is now available as an alternative analysis (analysis.lpca in config)
     scaler = StandardScaler(with_std=False)
     scaler.fit(X)  # compute mean inclusion rate per edge
     X_scaled = scaler.transform(X)