---
title: "CRC"
author: "Di Xiao"
output:
  BiocStyle::html_document
vignette: >
  %\VignetteIndexEntry{CRC}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


This vignette demonstrates the SwitchClass workflow using a colorectal cancer (CRC) plasma proteomics dataset. It illustrates how to preprocess the data, perform label-switch classification to distinguish baseline-aligned and perturbation-aligned molecular responses, and visualize results using UMAP, quadrant plots, boxplots, and pathway enrichment.

## Environment setup
```{r setup, warning=FALSE, message=FALSE}
library(SwitchClass)
library(ggplot2)
library(ggpubr)

library(limma)
library(PhosR)
```

## 1. Load and Prepare CRC Dataset
Load the preprocessed expression matrix and associated sample metadata from the CRC plasma proteome study. The dataset includes healthy control (HC), pre-treatment (Pre), and post-treatment (Post) groups.  
Apply quality filtering to retain proteins consistently detected across samples, followed by imputation of missing values using PhosR utilities to obtain a complete expression matrix.
```{r load-crc, fig.width=10, fig.height=5, warning=FALSE, message=FALSE}
data("CRC_subset")

# log2 transform
expr.log <- log2(as.matrix(CRC_subset$expr) + 1)
rownames(expr.log) <- rownames(CRC_subset$expr)

# filtering
group <- CRC_subset$group
crc.filtered <- PhosR::selectGrps(expr.log, group, percent = 0.1, n = 3)
tmp <- apply(crc.filtered, 2, as.numeric)
rownames(tmp) <- rownames(crc.filtered)
dim(tmp)

# imputation
crc.imputed <- PhosR::tImpute(
  PhosR::scImpute(mat = tmp, percent = 0.5, grps = group),
  m = 1.8,
  s = 0.3
)

labels <- factor(colnames(crc.imputed), levels = colnames(crc.imputed)[order(group)])

p1 <- PhosR::plotQC(crc.filtered, panel = "pca", grps = group, labels = labels) +
  ggplot2::scale_color_manual(values = c(HC = "#5AAF59", Pre = "#D3542D", Post = "#E5CA4E")) +
  ggplot2::theme_classic() +
  ggplot2::theme(
    panel.background = ggplot2::element_blank(),
    plot.background = ggplot2::element_blank(),
    panel.border = ggplot2::element_blank(),
    panel.grid = ggplot2::element_blank()
  )

p2 <- PhosR::plotQC(crc.imputed, panel = "pca", grps = group, labels = labels) +
  ggplot2::scale_color_manual(values = c(HC = "#5AAF59", Pre = "#D3542D", Post = "#E5CA4E")) +
  ggplot2::theme_classic() +
  ggplot2::theme(
    panel.background = ggplot2::element_blank(),
    plot.background = ggplot2::element_blank(),
    panel.border = ggplot2::element_blank(),
    panel.grid = ggplot2::element_blank()
  )

ggpubr::ggarrange(p1, p2, ncol=2, nrow=1, common.legend = TRUE, legend="right")
```

## 2. Differential Expression Analysis

Perform differential expression analysis between healthy controls and pre-treatment CRC samples using limma. Identify features showing significant fold-change and adjusted p-value thresholds.
```{r de-analysis, warning=FALSE, message=FALSE}

pp <- group[group %in% c("HC","Pre")]
design_DE <- model.matrix(~ 0 + pp)

cancer.pp <- crc.imputed[,group %in% c("HC","Pre")]
fit_DE  <- limma::lmFit(cancer.pp, design_DE)

cont_DE <- limma::makeContrasts(pre_vs_HC = ppPre - ppHC, levels = colnames(design_DE))
fit_DE2 <- limma::eBayes(limma::contrasts.fit(fit_DE, cont_DE))

DE.table <- limma::topTable(fit_DE2, number = Inf)
if (nrow(DE.table) == 0L) {
  stop("No differential expression results were returned by limma.")
}
DE.table$gene <- rownames(DE.table)
DE.table <- DE.table[order(DE.table$adj.P.Val, DE.table$P.Value), ]

DE.genes <- DE.table$gene[which((DE.table[, "adj.P.Val"] < 0.05) & (abs(DE.table[,"logFC"]) > 0.5849625))]
length(DE.genes)

```

## 3. Label-Switch Classification
Use the SwitchClass algorithm to train random forest classifiers under baseline-aligned and perturbation-aligned label schemes, then compute feature-level directional importance scores.
```{r switchclass, warning=FALSE, message=FALSE}
# Use DE genes to focus the classification
X_bc <- crc.imputed[DE.genes, , drop = FALSE]

# Label schemes: baseline-aligned vs perturbation-aligned
y_rev <- paste0("c", ifelse(group=="HC",1, ifelse(group=="Pre",2,1)))
y_per <- paste0("c", ifelse(group=="HC",1, ifelse(group=="Pre",2,2)))

# run SwitchClass
res <- SwitchClass::label_switch_classify(
    X_bc, y_reverse = y_rev, y_persist = y_per,
    ntimes = 3, ntree = 100, seed = 1
)
delta <- res$delta

```

## 4. UMAP Visualization of Feature Subsets
Visualize sample clustering based on UMAP embeddings computed using baseline-aligned, perturbation-aligned, or all features. Consistent color mapping highlights group separation.
```{r umap, warning=FALSE, message=FALSE, fig.width=12, fig.height=4}

res_trip <- umap_triptych_by_delta(
  X = X_bc,
  delta = delta,
  groups = group,
  top_n = 100,
  palette = c(HC = "#5AAF59", Pre = "#D3542D", Post = "#E5CA4E")
)

# View the combined plot
res_trip$plot


```

## 5. Quadrant Assignment and Delta–Fold-Change Scatter Plot
Combine directional importance scores with log fold-changes to assign each feature to one of four quadrants, reflecting distinct molecular response behaviors. Plot the scatter to reveal major feature trends.
```{r quadrants, warning=FALSE, message=FALSE, fig.width=6, fig.height=6}
Hm <- matrixStats::rowMeans2(as.matrix(X_bc[, group == "HC", drop = FALSE]))
Dm <- matrixStats::rowMeans2(as.matrix(X_bc[, group == "Pre", drop = FALSE]))
lfc <- Hm - Dm   

out <- plot_delta_vs_fc(delta, lfc, delta_thresh = 0.15, label_top = 8)

df_quad <- out$data  # tidy table with quadrants
out$plot  # the scatter plot

```

## 6. Expression Patterns of Representative Features
Generate boxplots for selected representative features from each quadrant to illustrate characteristic expression changes across groups.
```{r boxplots, warning=FALSE, message=FALSE, fig.width=10, fig.height=5}
genes_order <- c("RASA4","RASA4B","SMYD5","MRPS25","BLM","CPSF6","HBB","MSL1")
palette_crc <- c(HC = "#5AAF59", Pre = "#D3542D", Post = "#E5CA4E")
group <- factor(group, levels = c("HC","Pre", "Post"))

p <- plot_feature_boxplots(
  X = X_bc,
  group = group,
  features = genes_order,
  feature_order = genes_order,
  palette = palette_crc,
  ncol = 4,
  title = "Top example genes across groups"
)
p

```
## 7. Heatmap of Quadrant-Defined Features
Create a heatmap showing standardized expression of genes across samples, grouped by both quadrants and sample type (HC, Pre, Post).
```{r heatmap, fig.width=8, fig.height=5, warning=FALSE, message=FALSE}
plot_heatmap_quadrants(
  X = crc.imputed,
  group = group,
  df_quad = df_quad,
  cutoff = 0.15,
  palette_group = c(HC = "#5AAF59", Pre = "#D3542D", Post = "#E5CA4E")
)

```

## 8. Pathway Enrichment Analysis (Reactome)
Perform over-representation analysis using Reactome gene sets from msigdbr to identify biological pathways enriched within each quadrant. Visualize the top pathways per quadrant as -log10(p-value) barplots.
```{r enrichment, fig.width=8, fig.height=6, warning=FALSE, message=FALSE, eval=FALSE}
out <- enrich_reactome_quadrants(
  df = df_quad,
  universe = rownames(CRC_subset$expr),
  top_terms = 5
)

out$plot      # the combined 4-panel barplot
# out$results$Q1  # raw over-representation table for Q1
# out$barplot_df  # tidy plotted values
    

```

## Session information

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