---
title: "Analyzing human RNA-Seq data with CLAMP"
author: "Marc Subirana-Granés"
package: CLAMP
output: 
  BiocStyle::html_document:
    fig_width: 7
    fig_height: 5
vignette: >
  %\VignetteIndexEntry{Analyzing human RNA-Seq data with CLAMP}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    eval = TRUE,
    warning = FALSE,
    message = FALSE,
    fig.align = "center",
    out.width = "80%",
    dev = "png"
)

options(bitmapType = "cairo")
set.seed(1)
```

```{r logo, echo=FALSE, eval=TRUE, out.width='40%'}
knitr::include_graphics("../man/figures/clamp.png", dpi = 800)
```

# Installation

Install the released version of CLAMP from Bioconductor:

```{r install, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("CLAMP")
```

The development version can be installed from GitHub:

```{r install-dev, eval=FALSE}
BiocManager::install("chikinalab/CLAMP")
```

# Load packages

```{r load-packages}
library(CLAMP)
library(CLAMPData)
library(dplyr)
library(rsvd)
library(glmnet)
library(Matrix)
library(rhdf5)
library(data.table)
library(bigstatsr)
library(here)
library(AnnotationDbi)
library(org.Hs.eg.db)
library(DT)
library(DiagrammeR)
```

# Introduction

The CLAMP (**C**urated **L**atent-variable **A**nalysis with **M**olecular **P**riors) 
package provides a two-stage framework to extract interpretable latent
variables from high-dimensional transcriptomic data.
It combines a standard matrix decomposition (**CLAMPbase**) with
pathway-guided factor refinement (**CLAMPfull**), enabling:

1. **Dimensionality reduction** of gene expression matrices  
2. **Incorporation of prior knowledge** (e.g., Gene Ontology or MSigDB)  
3. **Adaptive regularization** of gene weights based on their agreement with pathway priors  
4. **Cross-validation** to select robust latent variables  

In CLAMPfull, pathway information is integrated through an adaptive **variance prior** that dynamically
modulates the contribution of each gene according to how well its latent signal aligns with pathway predictions.
This mechanism allows CLAMP to emphasize biologically consistent genes while maintaining flexibility
to discover novel, data-driven components.  
By combining prior-guided regularization with scalable matrix updates, CLAMPfull produces interpretable
latent variables that capture both known and emergent biological processes across large transcriptomic datasets.

## Workflow overview

```{r CLAMP-workflow, echo=FALSE, warning=FALSE, message=FALSE, fig.width=12, fig.height=9, out.width='100%'}
grViz('
 digraph CLAMP_pipeline {
   graph [
     layout=dot, rankdir=TB, splines=ortho,
     ranksep=0.6, nodesep=0.45, margin=0.03
   ]
   node [fontname=Helvetica, fontsize=16, margin="0.10,0.05", penwidth=1.2]
   edge [arrowsize=0.7, penwidth=1.1]

   // ---- Inputs & Priors (gray) ----
   DF     [shape=folder,   style=filled, fillcolor=Gray90, label="data.frame"]
   H5     [shape=folder,   style=filled, fillcolor=Gray90, label="H5 file"]
   GO [shape=cylinder, style=filled, fillcolor=Gray90,
       label="GO Biological Process"]
   MSIGDB [shape=cylinder, style=filled, fillcolor=Gray90,
           label="MSigDB Hallmark"]
   GTEx   [shape=cylinder, style=filled, fillcolor=Gray90, label="GTEx Tissues"]
   GMTFILE [shape=folder,  style=filled, fillcolor=Gray90, label="GMT file"]

   // ---- Functions (white) ----
   node [shape=rectangle, style=filled, fillcolor=White]
   cpmDF      [label="cpmCLAMP()"]
   preDF      [label="preprocessCLAMP()"]
   zDF        [label="zscoreCLAMP()"]
   selKMem    [label="select_svd_k()"]
   svdMem     [label="compute_svd()"]
   cpmFBM     [label="cpmCLAMPFBM()"]
   preFBM     [label="preprocessCLAMPFBM()"]
   zFBM       [label="zscoreCLAMPFBM()"]
   selKFBM    [label="select_svd_k()"]
   svdFBM     [label="compute_svd()"]
   numPC      [label="select_clamp_k()"]
   getGMTfunc [label="getGMT()"]
   readGMT    [label="read_gmt()"]
   gmtList    [label="gmtListToSparseMat()"]
   matchPaths [label="getMatchedPathwayMat()"]
   base       [label="CLAMPbase()"]
   full       [label="CLAMPfull()"]

   // ---- Outputs (light-blue tabs) ----
   node [shape=tab, style=filled, fillcolor=LightBlue]
   cpmOut       [label="Y_cpm"]
   filtered     [label="Y_filtered, rowStats"]
   zscored      [label="Y_zscored"]
   svdK         [label="svd_k"]
   svdRes       [label="svdres"]
   clampK       [label="clamp_k"]
   cpmFBMOut    [label="FBM_cpm"]
   filteredFBM  [label="FBM_filtered, rowStats"]
   zscoredFBM   [label="FBM_zscored"]
   svdKFBM      [label="svd_k"]
   svdResFBM    [label="svdres"]
   pathMat      [label="pathMat"]
   matchedPaths [label="matchedPathways"]
   baseResult   [label="CLAMPbase.result"]
   Zmat         [label="Z"]
   Bmat         [label="B"]
   summaryTbl   [label="summary"]
   Umat         [label="U"]

   // ---- In-memory branch ----
   DF -> cpmDF -> cpmOut -> preDF -> filtered -> zDF -> 
        zscored -> selKMem -> svdK
   svdK -> svdMem
   zscored -> svdMem -> svdRes
   svdRes -> numPC
   svdK -> numPC
   numPC -> clampK
   svdRes -> base
   clampK -> base
   zscored -> base
   base -> baseResult -> full

   // ---- On-disk branch ----
   H5 -> cpmFBM -> cpmFBMOut -> preFBM -> filteredFBM -> 
        zFBM -> zscoredFBM -> selKFBM -> svdKFBM
   svdKFBM -> svdFBM
   zscoredFBM -> svdFBM -> svdResFBM
   svdResFBM -> numPC
   svdKFBM -> numPC
   svdResFBM -> base
   zscoredFBM -> base

   // ---- Prior-knowledge branch ----
   GO -> getGMTfunc
   MSIGDB -> getGMTfunc
   GTEx -> getGMTfunc
   GMTFILE -> readGMT
   getGMTfunc -> gmtList
   readGMT -> gmtList
   gmtList -> pathMat -> matchPaths -> matchedPaths
   matchedPaths -> full

   // ---- CLAMPfull outputs ----
   full -> Zmat
   full -> Bmat
   full -> summaryTbl
   full -> Umat
 }',
    width = "100%",
    height = "900px"
)
```

# Quick Start

## Available datasets

CLAMPData ships three curated datasets used throughout this vignette. You can
inspect them with `list_clamp_data()`:

```{r list-clamp-data}
list_clamp_data()
```

We provide three examples:

1. **Data-frame example (whole blood):**  
    A small dataset loaded entirely into memory. Shows basic preprocessing,
    z-scoring, and running CLAMP without on-disk storage.  

2. **HDF5 example (Alzheimer’s brain):**  
    Demonstrates how to import expression from an HDF5 file, create a file-backed FBM object,
    and process larger datasets using the FBM interface.

3. **Table example (pancreatic islets):**  
    Illustrates reading a tab-delimited count file and comparing conditions via the B matrix.

Each example follows these steps:

1. **Load & preprocess data** (filter by mean/variance, then z-score).
2. **Compute SVD** and **infer the model dimension `k`**.
3. **Prepare pathway annotations** via `getGMT()` and construct the prior matrix.
4. **Run CLAMPbase** with the pre-computed SVD result and `k` to initialize the latent variables.
5. **Run CLAMPfull** to refine latent variables using pathway priors and variance-adaptive regularization, producing the final model and summary statistics.

## Computing the SVD and inferring k

CLAMP requires a truncated Singular Value Decomposition (SVD) of the z-scored expression matrix as input.
The choice of SVD function depends on dataset size:

- **Small to medium datasets (in-memory):** Use `rsvd::rsvd()` from the **rsvd** package. This is efficient for matrices that fit comfortably in RAM.

- **Large datasets (file-backed):** Use `bigstatsr::big_randomSVD()` for file-backed matrices (FBM). This function computes the SVD without loading the entire matrix into memory, enabling analysis of datasets too large for RAM.

After computing the SVD, infer the optimal number of latent variables
(`clamp_k`) using `num.pc()`:

# Example 1: In-memory data.frame (Whole Blood)

This example uses whole-blood RNA-seq data available from GEO under accession
[GSE130824](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE130824)
(Homo sapiens, 36 samples). It demonstrates the standard in-memory CLAMP
workflow for bulk transcriptomic data from peripheral blood.

## Load example data

In this chunk, we load the whole-blood expression matrix.

```{r load-data, message=FALSE}
data("dataWholeBlood") # expression matrix
dim(dataWholeBlood) # genes x samples
dataWholeBlood[1:6, 1:6] # genes x samples
```

## Preprocess and z-score

We first CPM-normalize the data (when needed), filter for genes with mean 
expression ≥ 0.5 and variance ≥ 0.1, and then apply z-score normalization.

```{r preprocess, message=TRUE}
#  CPM normalization
dataWholeBlood_cpm <- cpmCLAMP(dataWholeBlood)

# Filter and compute row statistics
prep_wb <- preprocessCLAMP(
    Y = dataWholeBlood_cpm,
    mean_cutoff = 0.5,
    var_cutoff = 0.1
)

# Extract filtered matrix and rowStats
wb_Y_filtered <- prep_wb$Y_filtered
wb_rowStats <- prep_wb$rowStats


# Z-score normalization
wb_Y_z <- zscoreCLAMP(
    Y_filtered = wb_Y_filtered,
    rowStats = wb_rowStats
)
```

## Compute SVD and infer k

We compute the SVD using `select_svd_k()` and `compute_svd()`, then select `clamp_k` with `select_clamp_k()`.

```{r wb-svd, message=FALSE}
# Select SVD rank and compute SVD
wb_svd_k   <- select_svd_k(wb_Y_z)
wb_svd     <- compute_svd(wb_Y_z, k = wb_svd_k)

# Select clamp_k (elbow method by default)
wb_clamp_k <- select_clamp_k(wb_svd, n_samples = ncol(wb_Y_z), svd_k = wb_svd_k)
wb_clamp_k
```

## CLAMPbase initialization

We initialize latent variables using `CLAMPbase`, providing the pre-computed SVD
and inferred `k`.

The argument `adaptive.p` defines the percentile used to determine the adaptive
sparsity threshold applied to each latent variable's gene loadings.
During alternating updates, negative entries in `Z` are treated as noise, and 
CLAMP estimates a cutoff based on the `adaptive.p` quantile of these negative
values. All genes with loadings below this cutoff are set to zero.

This produces **data-driven sparsity**, automatically filtering weak or noisy 
signals while retaining genes with the strongest positive contributions.
Lower values of `adaptive.p` (e.g., 0.01) result in stronger sparsity, while 
higher values (e.g., 0.1) retain more genes.
The default `adaptive.p = 0.05` typically yields interpretable, well-separated 
latent variables in large transcriptomic datasets.

```{r CLAMPbase, message=FALSE}
wb_baseRes <- CLAMPbase(
    Y = wb_Y_z,
    svdres = wb_svd,
    clamp_k = wb_clamp_k
)
```

## Prepare pathway priors

Next, we build a prior matrix from curated gene sets and compute the Chat 
object for `CLAMPfull`.

```{r priors-download, eval=FALSE}
# How to download pathway and cell marker libraries from Enrichr.
# Not run during vignette build to avoid network calls; pre-fetched
# .rds files are loaded in the next chunk instead.
enrichr_url <- "https://maayanlab.cloud/Enrichr/geneSetLibrary"
gmtList <- list(
    CellMarkers = getGMT(
        paste0(enrichr_url, "?mode=text&libraryName=CellMarker_2024"),
        "CellMarker_2024"
    ),
    KEGG = getGMT(
        paste0(enrichr_url, "?mode=text&libraryName=KEGG_2021_Human"),
        "KEGG_2021_Human"
    )
)
```

```{r priors, message=FALSE}
# Load pre-fetched gene set libraries bundled with the package
gmtList <- list(
    CellMarkers = readRDS(
        system.file("extdata", "CellMarker_2024.rds", package = "CLAMP")
    ),
    KEGG = readRDS(
        system.file("extdata", "KEGG_2021_Human.rds", package = "CLAMP")
    )
)

# Combine into a single sparse matrix
pathMatCell <- gmtListToSparseMat(gmtList)

# Load additional xCell reference matrix
data("xCell")

# Match pathways to the gene space of whole blood
matchedPathsWB <- getMatchedPathwayMatList(
    pathMatCell, xCell,
    new.genes = rownames(dataWholeBlood),
    min.genes = 2
)
```

**Note**: GMT files can also be loaded from local storage using `read_gmt()`.
This allows you to integrate custom or curated gene set libraries, such as 
*MSigDB* canonical pathways, directly into your analysis pipeline alongside
remote resources.

## CLAMPfull

Finally, we refine the base model by integrating pathway priors using 
`CLAMPfull`, which applies cross-validation to optimize latent variable
regularization. In this new version, `CLAMPfull` incorporates 
**variable priors** that adjust the influence of each pathway adaptively, 
improving convergence and stability across heterogeneous datasets.

```{r CLAMPfull, message=FALSE}
wb_fullRes <- CLAMPfull(
    wb_Y_z,
    priorMat = matchedPathsWB,
    clamp.base.result = wb_baseRes,
    svdres = wb_svd,
    clamp_k = wb_clamp_k,
    use_cpp = TRUE
)
```

## Display significant latent variables

```{r summary-table}
# Display significant latent variables
wb_summary_df <- as.data.frame(wb_fullRes$summary) %>%
    dplyr::filter(FDR < 0.05 & AUC > 0.7) %>%
    dplyr::arrange(FDR) %>%
    dplyr::select(LV, pathway, FDR, AUC)

datatable(
    wb_summary_df,
    filter = "top",
    options = list(
        pageLength = 10,
        autoWidth  = TRUE
    ),
    rownames = FALSE,
    class = "stripe hover compact"
) %>%
    formatSignif(c("AUC", "FDR"), 3)
```

The recovered LVs are biologically coherent for whole blood. LV13 aligns with neutrophil
signatures, LV10 with platelets, LV14 with erythrocytes, LV12 with NK cells,
and LV11 with plasma cells, covering the major cellular constituents of whole
blood. Together, these results indicate that CLAMP successfully decomposes the
bulk transcriptomic signal into its dominant blood-cell-type components.

# Example 2: File-Backed Matrix (Alzheimer’s Brain)

This example uses data from Alzheimer’s brain samples from a *Neurobiology of Disease* study (Barbash et al., 2017; DOI: https://doi.org/10.1016/j.nbd.2017.06.008). It demonstrates the on‑disk workflow with a file‑backed FBM to handle large‑scale transcriptomic datasets.

## Cleanup old FBM files

```{r fbm-cleanup, message=FALSE}
output_dir <- here("output", "alzFBM")
fbm_base <- file.path(output_dir, "FBMalz")
bk_paths <- paste0(fbm_base, c(".bk", "_preproc.bk", "_preproc_filtered.bk"))
file.remove(bk_paths[file.exists(bk_paths)])
```

### Computing CPM on a File-Backed Matrix (FBM)

For file-backed matrices (FBMs), you can compute counts-per-million (CPM) 
in-place—without loading the entire dataset into RAM—using the 
`cpmCLAMPFBM()` function from **CLAMP**:

## HDF5 schema

CLAMP HDF5 files follow a fixed layout. You can inspect the expected structure with `clamp_h5_schema()`:

```{r alz-schema}
clamp_h5_schema()
```

## Load HDF5 expression

`read_clamp_alz_expression()` downloads the file via ExperimentHub, validates it
against the schema, and returns a genes × samples matrix with row and column names ready to use.

```{r hdf5-load, message=FALSE}
expr_mat <- read_clamp_alz_expression()
genes <- rownames(expr_mat)
```

## Construct file‑backed FBM

```{r fbm-construct, message=FALSE}
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
alzFBM <- FBM(
    nrow = nrow(expr_mat), ncol = ncol(expr_mat),
    backingfile = fbm_base
)
blk <- 1000

for (i in seq_len(ceiling(nrow(expr_mat) / blk))) {
    rows <- ((i - 1) * blk + 1):min(i * blk, nrow(expr_mat))
    alzFBM[rows, ] <- expr_mat[rows, , drop = FALSE]
}
```

## CPM, preprocess and z‑score FBM

```{r fbm-preprocess, message=FALSE}
prep_alz <- preprocessCLAMPFBM(
    fbm = alzFBM,
    mean_cutoff = 0.5,
    var_cutoff = 0.1
)

alz_fbm_filt <- prep_alz$fbm_filtered
alz_rowStats <- prep_alz$rowStats
zscoreCLAMPFBM(alz_fbm_filt, alz_rowStats)
alz_genes <- genes[prep_alz$kept_rows]
```

## Compute SVD and infer k

For file-backed matrices, `compute_svd()` dispatches to `bigstatsr::big_SVD()` automatically, avoiding loading the entire matrix into RAM.

```{r fbm-svd, message=FALSE}
# Select SVD rank and compute SVD (dispatches to bigstatsr for FBM)
alz_svd_k   <- select_svd_k(alz_fbm_filt)
alz_svd     <- compute_svd(alz_fbm_filt, k = alz_svd_k)

# Select clamp_k (elbow method by default)
alz_clamp_k <- select_clamp_k(alz_svd, n_samples = ncol(alz_fbm_filt),
                              svd_k = alz_svd_k)
alz_clamp_k
```

## CLAMPbase

```{r fbm-CLAMPbase, message=FALSE}
alz_baseRes <- CLAMPbase(
    Y = alz_fbm_filt,
    svdres = alz_svd,
    clamp_k = alz_clamp_k
)
```

## Prepare pathway priors

```{r fbm-priors-download, eval=FALSE}
# How to fetch the libraries; not run during vignette build.
enrichr_url <- "https://maayanlab.cloud/Enrichr/geneSetLibrary"
alz_gmtList <- list(
    GTEx_Tissues = getGMT(
        paste0(enrichr_url, "?mode=text&libraryName=GTEx_Tissues_V8_2023")
    ),
    BP = getGMT(
        paste0(enrichr_url, "?mode=text&libraryName=GO_Biological_Process_2025")
    ),
    MSigDB = getGMT(
        paste0(enrichr_url, "?mode=text&libraryName=MSigDB_Hallmark_2020")
    )
)
```

```{r fbm-priors, message=FALSE}
alz_gmtList <- list(
    GTEx_Tissues = readRDS(
        system.file("extdata", "GTEx_Tissues_V8_2023.rds", package = "CLAMP")
    ),
    BP = readRDS(
        system.file(
            "extdata", "GO_Biological_Process_2025.rds",
            package = "CLAMP"
        )
    ),
    MSigDB = readRDS(
        system.file("extdata", "MSigDB_Hallmark_2020.rds", package = "CLAMP")
    )
)

alz_pathMat <- gmtListToSparseMat(alz_gmtList)
alz_matched <- getMatchedPathwayMat(alz_pathMat, alz_genes)
```

## CLAMPfull

```{r fbm-CLAMPfull, message=FALSE}
alz_fullRes <- CLAMPfull(
    alz_fbm_filt,
    priorMat = alz_matched,
    clamp.base.result = alz_baseRes,
    svdres = alz_svd,
    clamp_k = alz_clamp_k,
    use_cpp = TRUE
)
```

## Display significant latent variables

```{r fbm-summary}
alz_summary_df <- as.data.frame(alz_fullRes$summary) %>%
    dplyr::filter(FDR < 0.05 & AUC > 0.7) %>%
    dplyr::arrange(FDR) %>%
    dplyr::select(LV, pathway, FDR, AUC)

datatable(
    alz_summary_df,
    filter = "top",
    options = list(
        pageLength = 10,
        autoWidth  = TRUE
    ),
    rownames = FALSE,
    class = "stripe hover compact"
) %>%
    formatSignif(c("AUC", "FDR"), 3)
```

The significant LVs align with brain-relevant transcriptional programs implicated in Alzheimer’s disease.
LV3 and LV10 are enriched for GTEx brain-region signatures, including spinal cord, substantia nigra,
frontal cortex, and cortex, suggesting that these axes capture genuine neural transcriptional variation.
LV1 further supports disease relevance through enrichment for mitochondrial respiration and oxidative phosphorylation pathways,
which are linked to impaired brain energy metabolism in Alzheimer’s disease.

# Example 3: Tab-Delimited Count File (Pancreatic Islets)

In this example, we apply the in‑memory CLAMP workflow to RNA‑Seq count data from GEO accession **GSE164416** 
(Wigger *et al.* 2021; “Multi‑omics profiling of living human pancreatic islet donors reveals heterogeneous beta-cell
trajectories towards type 2 diabetes”, DOI: 10.1038/s42255-021-00420-9). After preprocessing the raw counts and fitting
the CLAMP model, we perform a differential analysis of latent‑variable activities to compare non‑diabetic (ND) and type 2
diabetic (T2D) samples.

## Load count data and map gene symbols

```{r islet-load, message=FALSE}
islet_df <- read_islet_counts()

islet_df$symbol <- mapIds(org.Hs.eg.db,
    keys = islet_df$ensembl,
    column = "SYMBOL",
    keytype = "ENSEMBL",
    multiVals = "first"
)

islet_df <- islet_df[!is.na(islet_df$symbol), ]
```

## Aggregate counts by gene symbol

```{r islet-aggregate, message=FALSE}
# Sum counts per symbol
setDT(islet_df)
num_cols <- names(islet_df)[sapply(islet_df, is.numeric)]
expr <- islet_df[, lapply(.SD, sum), by = symbol, .SDcols = num_cols]
expr <- as.data.frame(expr)
rownames(expr) <- expr$symbol
expr$symbol <- NULL
expr <- as.matrix(expr)
```

## CPM, preprocess and z-score

```{r islet-preprocess, message=FALSE}
prep_is <- preprocessCLAMP(
    Y = expr,
    mean_cutoff = 0.5,
    var_cutoff = 0.1
)

iso_Yf <- prep_is$Y_filtered
iso_rowS <- prep_is$rowStats

iso_Yz <- zscoreCLAMP(
    Y_filtered = iso_Yf,
    rowStats = iso_rowS
)
```

## Compute SVD and infer k

```{r islet-svd, message=FALSE}
# Select SVD rank and compute SVD
islet_svd_k   <- select_svd_k(iso_Yz)
islet_svd     <- compute_svd(iso_Yz, k = islet_svd_k)

# Select clamp_k (elbow method by default)
islet_clamp_k <- select_clamp_k(islet_svd, n_samples = ncol(iso_Yz),
                                svd_k = islet_svd_k)
islet_clamp_k
```

## CLAMPbase

```{r islet-CLAMPbase, message=FALSE}
islet_baseRes <- CLAMPbase(
    Y = iso_Yz,
    svdres = islet_svd,
    clamp_k = islet_clamp_k
)
```

## Prepare pathway priors

```{r islet-priors-download, eval=FALSE}
# How to fetch the libraries; not run during vignette build.
enrichr_url <- "https://maayanlab.cloud/Enrichr/geneSetLibrary"
islet_gmtList <- list(
    GTEx_Tissues = getGMT(
        paste0(enrichr_url, "?mode=text&libraryName=GTEx_Tissues_V8_2023")
    ),
    Diabetes_Perturbations = getGMT(
        paste0(
            enrichr_url,
            "?mode=text&libraryName=Diabetes_Perturbations_GEO_2022"
        )
    ),
    MSigDB_Hallmark = getGMT(
        paste0(enrichr_url, "?mode=text&libraryName=MSigDB_Hallmark_2020")
    )
)
```

```{r islet-priors, message=FALSE}
islet_gmtList <- list(
    GTEx_Tissues = readRDS(
        system.file("extdata", "GTEx_Tissues_V8_2023.rds", package = "CLAMP")
    ),
    Diabetes_Perturbations = readRDS(
        system.file(
            "extdata", "Diabetes_Perturbations_GEO_2022.rds",
            package = "CLAMP"
        )
    ),
    MSigDB_Hallmark = readRDS(
        system.file("extdata", "MSigDB_Hallmark_2020.rds", package = "CLAMP")
    )
)

islet_pathMat <- gmtListToSparseMat(islet_gmtList)
islet_matched <- getMatchedPathwayMat(islet_pathMat, rownames(iso_Yz))
islet_chatObj <- getChat(islet_matched)
```

## CLAMPfull

```{r islet-CLAMPfull, message=FALSE}
islet_fullRes <- CLAMPfull(
    iso_Yz,
    priorMat = islet_matched,
    clamp.base.result = islet_baseRes,
    svdres = islet_svd,
    clamp_k = islet_clamp_k,
    use_cpp = TRUE
)
```

## Display significant latent variables

```{r islet-summary}
islet_summary_df <- as.data.frame(islet_fullRes$summary) %>%
    dplyr::filter(FDR < 0.05 & AUC > 0.7) %>%
    dplyr::arrange(FDR) %>%
    dplyr::select(LV, pathway, FDR, AUC)

datatable(
    islet_summary_df,
    filter = "top",
    options = list(
        pageLength = 10,
        autoWidth  = TRUE
    ),
    rownames = FALSE,
    class = "stripe hover compact"
) %>%
    formatSignif(c("AUC", "FDR"), 3)
```

The significant LVs reflect key biological processes relevant to type 2 diabetes. 
LV20 and LV16 capture alpha- and beta-cell identity programs, highlighting pancreatic islet endocrine biology. 
LV21 aligns with pancreas-specific GTEx tissue signatures, supporting tissue relevance, 
while LV1 captures oxidative phosphorylation, protein secretion, and beta-cell-related programs. 

## Differential latent-variable expression between conditions

Rows of the B matrix correspond to LVs and columns to samples. By grouping 
samples by condition (ND vs T2D), `differentialLVActivity()` computes average
LV expression per group and tests for LVs that differ between healthy and
diabetic islets.

```{r islet-diff, message=FALSE}
islet_metadata <- read_islet_metadata()

lv_stats_all_vs_nd <- differentialLVActivity(
    islet_fullRes,
    metadata = islet_metadata,
    sample_col = "id",
    group_col = "type",
    reference = "ND"
)

sig_lv_all_vs_nd <- lv_stats_all_vs_nd %>%
    dplyr::filter(FDR < 0.1)

sig_pathway <- islet_summary_df %>%
    dplyr::filter(FDR < 0.05 & AUC > 0.7) %>%
    dplyr::filter(LV %in% sig_lv_all_vs_nd$LV) %>%
    dplyr::arrange(FDR) %>%
    dplyr::select(LV, pathway, FDR, AUC)

datatable(
    sig_pathway,
    filter = "top",
    options = list(
        pageLength = 10,
        autoWidth  = TRUE
    ),
    rownames = FALSE,
    class = "stripe hover compact"
) %>%
    formatSignif(c("AUC", "FDR"), 3)
```

The top differentially active LVs highlight biological axes distinguishing T2D from ND islets. 
LV7 links to diabetic adipose tissue and TNF-alpha signaling via NF-kB, consistent with inflammation and metabolic dysfunction.
LV20 and LV16 map to alpha- and beta-cell programs. 
LV9 and LV3 are associated with islet perturbation and diabetic mouse islet signatures, 
supporting disease-relevant changes in islet transcriptional states,
while LV10 suggests a vascular component relevant to T2D.

# Projection: Applying one CLAMP model to another dataset

`projectCLAMP()` reuses the gene loadings (**Z**) from a fitted CLAMP model
and estimates latent-variable activities (**B**) for a new expression matrix.
Projection uses the same genes in the same order; when both matrices have row
names, `projectCLAMP()` aligns the common genes automatically before solving
for **B**.

Here we project the whole-blood expression matrix from Example 1 into the
full latent-variable space learned from the pancreatic islet model in Example 3.

```{r project-islet-model-to-whole-blood, message=TRUE}
islet_model_genes <- rownames(islet_fullRes$Z)
wb_project_genes <- rownames(wb_Y_z)

common_genes <- intersect(islet_model_genes, wb_project_genes)
cat(
    "Overlapping genes:", length(common_genes), "/", length(islet_model_genes),
    "islet model genes",
    sprintf(
        "(%.1f%%)\n",
        100 * length(common_genes) / length(islet_model_genes)
    )
)

# projectCLAMP aligns common row names in the model's gene order
wb_projected_B <- projectCLAMP(islet_fullRes, wb_Y_z)

dim(wb_projected_B)
wb_projected_B[
    seq_len(min(5, nrow(wb_projected_B))),
    seq_len(min(5, ncol(wb_projected_B))),
    drop = FALSE
]
```

# Choosing the Number of Latent Variables (CLAMP_K)

CLAMP_K controls how many latent variables the model learns.  Too few and
biologically distinct signals merge; too many and noise is absorbed into
spurious components.  `select_clamp_k()` is the unified interface: it takes
the SVD result, the number of samples, the SVD truncation rank, and an
optional `method` argument, and returns a list with `$clamp_k` (number of LVs)
and `$scale` (regularization scale used downstream).

## Elbow method (default)

The elbow heuristic fits a smoothing spline to the singular-value scree plot
and returns the index at which curvature is maximised.  This is the fastest
option and works well when the signal-to-noise boundary is clear.

```{r k-elbow}
select_clamp_k(
    wb_svd,
    n_samples = ncol(wb_Y_z),
    svd_k     = wb_svd_k,
    method    = "elbow"
)
```

## Permutation method

The permutation approach shuffles each row of the input matrix independently
`B` times and recomputes the SVD to build a null distribution of singular
values.  The number of components whose observed singular value exceeds the
95th percentile of the null is returned.  This is more conservative and
slower, but robust to smooth scree plots.

```{r k-permutation, eval=FALSE}
select_clamp_k(
    wb_svd,
    n_samples = ncol(wb_Y_z),
    svd_k     = wb_svd_k,
    method    = "permutation",
    data      = wb_Y_z,
    B         = 2
)
```

## Gavish–Donoho optimal hard threshold (PCAtools)

The Gavish–Donoho threshold (Gavish & Donoho, 2014) identifies the
singular-value cutoff below which components are statistically
indistinguishable from noise, given matrix dimensions and an estimate of the
noise level. **PCAtools** implements this via `chooseGavishDonoho()`.

```{r k-gavish, eval=FALSE}
select_clamp_k(
    wb_svd,
    n_samples = ncol(wb_Y_z),
    svd_k     = wb_svd_k,
    method    = "gavish_donoho",
    data      = wb_Y_z
)
```

# Visualization

CLAMP provides dedicated plotting functions built on **ggplot2**, prefixed
`CLAMPplot` or `CLAMPdotplot`.  The examples below use the whole-blood result
`wb_fullRes` computed in Example 1.

## Pathway–LV association heatmap (`CLAMPplotU`)

`CLAMPplotU` displays the pathway loading matrix **U** after filtering by AUC
and FDR.  Only the top-`top` pathways per LV are shown, making it easy to
scan which pathways drive each latent variable.

```{r plot-U, message=FALSE, fig.width=14, fig.height=10}
CLAMPplotU(
    wb_fullRes,
    auc.cutoff = 0.6,
    fdr.cutoff = 0.05,
    top        = 3
)
```

## Top-gene loading plot (`CLAMPplotTopZ`)

`CLAMPplotTopZ` ranks genes by their Z loading for each selected LV and plots
the top genes as loading-versus-rank scatter plots.  The highest-loading genes
are labelled directly.

```{r plot-topZ, message=FALSE, fig.width=9, fig.height=7}
# Use the first few LVs that have pathway support
lv_with_paths <- wb_fullRes$withPrior[
    seq_len(min(4, length(wb_fullRes$withPrior)))
]

CLAMPplotTopZ(
    wb_fullRes,
    top       = 50,
    label.top = 10,
    index     = lv_with_paths
)
```

Only one LV:
```{r plot-topZ1, message=FALSE, fig.width=5, fig.height=3}
# Use the first few LVs that have pathway support
lv_with_paths <- wb_fullRes$withPrior[1]

CLAMPplotTopZ(
    wb_fullRes,
    top       = 50,
    label.top = 10,
    index     = lv_with_paths
)
```
## Single-LV pathway dot plot (`CLAMPdotplot`)

`CLAMPdotplot` shows the top pathways for one selected LV as a lollipop chart.
Dot size encodes AUC; dot colour encodes `-log10(FDR)`.  Use `x.axis` and
`order.by` to choose whether the x-axis and pathway ranking use AUC or
`-log10(FDR)`.

Plot order by AUC:

```{r plot-dot, message=FALSE, fig.width=7, fig.height=4, out.width='100%'}
CLAMPdotplot(
    wb_fullRes,
    lv         = "LV2",
    top        = 15,
    auc.cutoff = 0.6,
    fdr.cutoff = 0.1,
    x.axis     = "AUC",
    order.by   = "AUC"
)
```

Plot order by FDR:

```{r plot-dot-fdr, message=FALSE, fig.width=7, fig.height=4, out.width='100%'}
CLAMPdotplot(
    wb_fullRes,
    lv         = "LV2",
    top        = 15,
    auc.cutoff = 0.6,
    fdr.cutoff = 0.1,
    x.axis     = "-log10(FDR)",
    order.by   = "-log10(FDR)"
)
```

## All-LV pathway dot plot (`CLAMPdotplotAll`)

`CLAMPdotplotAll` gives a compact overview of all significant pathway–LV
associations across every latent variable. Dot size encodes AUC and dot colour
encodes `-log10(FDR)`.

```{r plot-dotAll, message=FALSE, fig.width=10, fig.height=7.5, out.width='100%'}
CLAMPdotplotAll(
    wb_fullRes,
    auc.cutoff = 0.65,
    fdr.cutoff = 0.05,
    top.per.lv = 5
)
```

# Parallelization in CLAMP

CLAMP supports multi-core parallelization for computationally intensive 
operations, particularly when working with large datasets 
and file-backed matrices (FBMs). The `ncores` parameter can be used in several 
key functions to speed up processing.

The following CLAMP functions accept an `ncores` parameter:

- `CLAMPbase()` 
- `CLAMPfull()`
- `projectCLAMP()` 
- `preprocessCLAMPFBM()`
- `zscoreCLAMPFBM()`
- `cpmCLAMPFBM()`

# Session Information

```{r session-info}
sessionInfo()
```