---
title: "Getting Started with conformeR"
author: "Justine Leclerc and Christof Seiler"
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
    number_sections: true
vignette: >
 %\VignetteIndexEntry{Getting Started with conformeR}
 %\VignetteEngine{knitr::rmarkdown}
 %\VignetteEncoding{UTF-8}
---

```{r options, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    dev = "png",
    fig.width = 12,
    fig.height = 7,
    out.width = "100%",
    warning = FALSE,
    message = FALSE
)
```

```{css, echo = FALSE}
.callout {
border-left: 4px solid #4472A8;
background-color: #F2F6FA;
padding: 0.6em 1em;
margin: 1em 0;
border-radius: 2px;
}
.callout p {
margin: 0.3em 0;
}
```

# Introduction

## Overview

LEMUR fits a low-dimensional embedding and assigns cells to neighborhoods 
based on similar responses to a given contrast and gene. 
Because these neighborhood assignments are themselves derived from a fitted 
model, they are subject to sampling variability: refitting the model on a 
different sample of cells can shift the neighborhood boundaries.

conformeR quantifies this uncertainty by wrapping the neighborhood construction
step in a conformal prediction framework. It is designed to integrate naturally
into existing Bioconductor workflows for single-cell data analysis, operating 
downstream of SingleCellExperiment-based preprocessing and LEMUR analyses to 
provide uncertainty-aware inference for neighborhood assignments.

It provides two complementary views of that uncertainty:

- **Conformal selection** asks, for a given gene, *which cells can we
confidently say belong to the affected neighborhood?* It returns a
subset of neighborhood cells with a guaranteed bound on the
proportion of false inclusions.
- **Conformal clustering** asks, *for each
cell, how uncertain are we about its neighborhood assignment?* It
returns a prediction *set* of plausible neighborhoods per cell,
which can contain zero, one, or several labels.

<div class="callout">
<p>**Why this matters.** A single LEMUR fit gives you one particular
neighborhood assignment, with no way to tell how much of it would
survive resampling. `conformeR` turns that single assignment into a
statement with a formal, finite-sample coverage guarantee.</p>
</div>

# Installation
Install conformeR running the following commands.
```{r installation, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}

BiocManager::install("juslecl/conformeR")
```


```{r setup and theme}
library(conformeR)
library(SingleCellExperiment)
library(dplyr)
library(ggplot2)
library(tibble)
library(knitr)

theme_set(
    theme_minimal(base_size = 13) +
        theme(
            panel.grid.minor = element_blank(),
            plot.title = element_text(face = "bold"),
            legend.position = "bottom"
        )
)
conformer_palette <- c(
    "inside" = "#4472A8",
    "outside" = "#D9822B",
    "none" = "#B0B0B0"
)
```

# Use-case of conformeR
## Example dataset

This vignette runs `conformeR` on a reduced version of the single-cell
RNA-seq dataset from Zhao, Dovas, Spinazzi, et al. (2021),
*"Deconvolution of cell type-specific drug responses in human tumor
tissue with single-cell RNA-seq,"* published in *Genome Medicine*
(GEO accession `GSE148842`). As with any input to `conformeR`, the
data have already been preprocessed and quality-checked upstream.

```{r load-data}
data("sce_reduced")
sce_reduced
rownames(sce_reduced) <- rowData(sce_reduced)$gene
head(colData(sce_reduced))
```

We restrict the analysis to three genes with distinct expected
behavior, to illustrate a range of neighborhood sizes and
uncertainty patterns: `HIST3H2A`, `PTPRZ1`, and `HBEGF`.

```{r genes-of-interest}
pan_affect_genes <- rowData(sce_reduced) |>
    as_tibble() |>
    dplyr::select(gid, gene) |>
    filter(gene %in% c("HIST3H2A", "PTPRZ1", "HBEGF"))

kable(pan_affect_genes, caption = "Genes of interest for this vignette")
```

## Running conformeR

`conformeR()` fits `LEMUR` and layers both conformal procedures on
top of it in one call. The key arguments are:

- `design_lemur` — the LEMUR design formula; here we adjust for
`patient_id` and estimate the effect of `condition`.
- `contrast_column` — which column of the design defines the
contrast of interest.
- `genes_of_interest` — restrict the (more expensive) conformal
layers to a gene subset; LEMUR itself is still fit genome-wide.
- `what` — which conformal procedure(s) to run: `"conf_selection"`,
`"conf_clustering"`, or both.

By default the target coverage is `cp_alpha`, and the data are split
50% training / 50% test, with the calibration set taking
25% of the training set (`cp_train_frac` and `cp_cal_frac`
set these proportions).

```{r run-conformer}
set.seed(2002)
res <- conformeR(
    sce_reduced,
    design_lemur = ~ patient_id + condition,
    contrast_column = "condition",
    genes_of_interest = pan_affect_genes$gene,
    what = c("conf_selection", "conf_clustering"),
    n_cores = 1,
    verbose = TRUE
)
```

## Conformal selection

For each gene, conformal selection returns the set of cells flagged
as belonging to the neighborhood affected by `condition`, at the
target coverage level. The table below reports both the raw
neighborhood size and its share of all cells in the test set.

```{r selection-summary}
fit <- res$fit_lemur
cp_out <- res$conf_results |>
    dplyr::mutate(
        inside_cs = purrr::map(inside_cs, ~ colnames(fit) %in% .x),
        inside_cc = purrr::map(inside_cc, ~ colnames(fit) %in% .x),
        outside_cc = purrr::map(outside_cc, ~ colnames(fit) %in% .x)
    ) |>
    tidyr::unnest(c(inside_cs, inside_cc, outside_cc)) |>
    dplyr::mutate(cell = rep(colnames(fit), nrow(fit)))

selection_summary <- cp_out |>
    group_by(gene) |>
    summarize(
        `Neighborhood size` = sum(inside_cs),
        `% of cells` = round(sum(inside_cs) / ncol(res$fit_lemur), 2) * 100,
        .groups = "drop"
    )

kable(selection_summary,
    caption =
        "Conformal selection: neighborhood size per gene"
)
```

## Conformal clustering

Conformal clustering instead reports, per cell, a *prediction set*:
the neighborhoods that cell could plausibly belong to at the target
coverage. A prediction set of size zero or two means the procedure could not
confidently assign that cell to any neighborhood. It is ambiguous.

```{r clustering-summary}
clustering_summary <- cp_out |>
    group_by(gene) |>
    summarize(
        `Prediction set = inside` =
            sum(inside_cc + outside_cc == 1 & inside_cc == 1),
        `Prediction set = outside` =
            sum(inside_cc + outside_cc == 1 & outside_cc == 1),
        `Prediction set size 2 or empty` = sum(inside_cc + outside_cc != 1),
        .groups = "drop"
    )

kable(clustering_summary,
    caption =
        "Conformal clustering: prediction set composition per gene"
)
```

## Visualizing results on the LEMUR embedding

### Conformal selection

```{r plot-selection, fig.cap = "Conformal selection overlaid on LEMUR"}
plotter_conformal_selection(res, pan_affect_genes$gene)
```

Cells shaded as part of the selected neighborhood are those the
procedure can confidently attribute to the condition effect for that
gene; unshaded cells fall outside the guaranteed selection at the
chosen coverage level.

### Conformal clustering

```{r plot-clustering, fig.cap = "Conformal prediction sets overlaid on LEMUR"}
plotter_conformal_clustering(res, pan_affect_genes$gene)
```

Here color reflects the prediction set itself rather than a single
hard assignment, making it easy to spot the boundary regions where
neighborhood membership is genuinely uncertain versus the core
regions where every resampling agrees.

# Session info

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