---
title: "cellpaintr Vignette"
author: 
- name: Christof Seiler
  affiliation: 
  - Center of Experimental Rheumatology, Department of Rheumatology, University Hospital Zurich, University of Zurich, Switzerland
- name: Phatthamon Laphanuwat
  affiliation: 
  - Center of Experimental Rheumatology, Department of Rheumatology, University Hospital Zurich, University of Zurich, Switzerland
- name: Caroline Ospelt
  affiliation: 
  - Center of Experimental Rheumatology, Department of Rheumatology, University Hospital Zurich, University of Zurich, Switzerland
date: "`r Sys.Date()`"
output: 
  BiocStyle::html_document
vignette: >
  %\VignetteIndexEntry{cellpaintr Vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
    echo = TRUE, warning = FALSE, message = FALSE,
    fig.retina = 2, dpi = 96,
    fig.width = 7.2916667, fig.asp = 0.6178571
)
```

# Introduction

High-content imaging can capture detailed morphological and functional features on a single-cell level. It has been widely used in drug screen assays and for generating large-scale image datasets for machine learning models. Software such as CellProfiler extract features from high-content images for downstream analysis. CellProfiler implements a suite of traditional image processing algorithms to segment cells and extract about 1,000 morphological and texture features.

To simplify data analysis of CellProfiler features in R---many options exist in Python, such as PyCytominer---we developed an R/Bioconductor package. Our package leverages Bioconductor objects for batch correction and dimension reduction, lowering the learning curve for users familiar with single-cell RNA sequencing workflows.

# Installation

Install this package.

```{r install-package, eval=FALSE}
if (!require("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("cellpaintr")
```

# Workflow

## Preparation

Load packages.

```{r load-libs}
library(cellpaintr)
library(scater)
library(scrapper)
library(purrr)
library(dplyr)
library(ggrepel)
```

## Load Data

Load `CellProfiler` output and create a `SingleCellExperiment` object.

```{r read-data}
set.seed(23)
cell_file <- generate_data()
sce <- loadData(cell_file)
```

## Data Cleaning and Transformation

Prepare the data for machine learning.

```{r data-preprocessing}
# remove cells with missing features
sce <- removeNAs(sce)

# number of cells
plotCellsPerImage(sce)

# remove outlier images
sce <- removeOutliers(sce, min = 0, max = 300)

# remove outlier cells
stats <- perCellQCMetrics(sce, assay.type = "features")
sce <- sce[, !isOutlier(stats$sum)]

# remove features with zero spread
sce <- removeLowVariance(sce, robust = TRUE)

# remove zero-inflated features
sce <- removeZeroInflation(sce)
```

Transform non-negatives features and scale all features.

```{r transform-data}
sce <- transformLogScale(sce, robust = TRUE)
```

## Unsupervised Analysis

Pseudo-bulk over images.

```{r aggreate}
aggr <- scrapper::aggregateAcrossCells(
    assay(sce, "tfmfeatures"),
    factors = colData(sce)[, c("ImageNumber", "Drug", "Patient")]
)
sce_aggr <- SingleCellExperiment(
    assays = list(sums = aggr$sums),
    colData = DataFrame(aggr$combinations, ncells = aggr$counts)
)
```

Use `scater` for exploratory data analysis. PCA on pseudo-bulk features.

```{r compute-pca}
sce_aggr <- scater::runPCA(sce_aggr,
    assay.type = "sums",
    ncomponents = 10
)
scater::plotReducedDim(sce_aggr, dimred = "PCA", colour_by = "Patient")
scater::plotReducedDim(sce_aggr, dimred = "PCA", colour_by = "Drug")
```

What is driving PC1?

```{r plot-pca}
plotPCACor(sce_aggr, filter_by = 1, assay_type = "sums")
```

UMAP.

```{r umap}
sce_aggr <- scater::runUMAP(sce_aggr, exprs_values = "sums")
scater::plotUMAP(sce_aggr, colour_by = "Patient")
scater::plotUMAP(sce_aggr, colour_by = "Drug")
```

## Supervised Analysis

Compare positive control to negative control.

```{r predict-single-drug}
set.seed(23)

sce$Drug <- as.factor(sce$Drug)
sce$Drug <- relevel(sce$Drug, ref = "D1")

types <- c("AreaShape", "Intensity", "Texture")

sce_single <- predictLOO(
    sce,
    target = "Drug", group = "Patient",
    interest_level = "D7", reference_level = "D1",
    types = types,
    n_threads = 1
)
```

Summarize feature subgroups in a volcano plot.

```{r plot-volcano}
volcanoPlot(sce_single,
    target = "Drug", group = "Patient",
    p_cutoff = 0.05, fc_cutoff = 0.5
)
```

Individual scores.

```{r plot-scores}
plotLOO(sce_single, target = "Drug", group = "Patient")
plotAUC(sce_single, target = "Drug", group = "Patient")
plotROC(sce_single, target = "Drug", group = "Patient")
```

Compare a list of drugs to negative control.

```{r predict-multiple-drugs}
set.seed(23)

reference_level <- "D1"
interest_levels <- setdiff(levels(sce$Drug), reference_level)

sce_list <- map(interest_levels, function(interest_level) {
    predictLOO(
        sce,
        target = "Drug", group = "Patient",
        interest_level = interest_level, reference_level = reference_level,
        types = types,
        n_threads = 1
    )
}, .progress = "drug screening")
```

Custom volcano plot for summarizing drugs faceted by features.

```{r plot-volcano-custom}
# prepare stats table
p_adj_cutoff <- 0.05
fc_cutoff <- 0.5
tb_stats <- lapply(sce_list, calculateStats,
    target = "Drug", group = "Patient"
) |>
    bind_rows() |>
    mutate(pvalue_adj = p.adjust(pvalue, method = "BH")) |>
    mutate(Selection = ifelse(
        pvalue_adj < p_adj_cutoff & log2FoldChange > fc_cutoff,
        Target, ""
    ))

# volcano plot
tb_stats |>
    ggplot(aes(log2FoldChange, -log10(pvalue_adj), label = Selection)) +
    geom_vline(xintercept = c(0, fc_cutoff), alpha = 0.5, linetype = "dashed") +
    geom_hline(
        yintercept = c(0, -log10(p_adj_cutoff)), alpha = 0.5,
        linetype = "dashed"
    ) +
    geom_point(aes(color = Target)) +
    geom_text_repel(aes(color = Target), max.overlaps = Inf) +
    xlab("log2 fold change") +
    ylab("-log10 p-value (BH adjusted)") +
    facet_wrap(~Feature) +
    ggtitle(paste0("False Discovery Rate: ", 100 * p_adj_cutoff, "%"))
```

# Session Info {-}

```{r session_info}
sessionInfo()
```