---
title: "3. Case Study: A Disease-Associated Monocyte State in COVID-19"
author:
    - name: Anthony Christidis
      affiliation:
        - &core_affiliation Core for Computational Biomedicine, Harvard Medical School
    - name: Andrew Ghazi
      affiliation:
        - *core_affiliation
    - name: Smriti Chawla
      affiliation:
        - *core_affiliation
    - name: Nitesh Turaga
      affiliation:
        - *core_affiliation
    - name: Ludwig Geistlinger
      affiliation:
        - *core_affiliation
    - name: Robert Gentleman
      affiliation:
        - *core_affiliation
package: scDiagnostics
output:
  BiocStyle::html_document:
    toc: true
    toc_float: true
vignette: >
  %\VignetteIndexEntry{3. Case Study: A Disease-Associated Monocyte State in COVID-19}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options:
  markdown:
    wrap: 72
---

```{r setup, include = FALSE}
knitr::knit_hooks$set(pngquant = knitr::hook_pngquant)

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  dev = "ragg_png",
  dpi = 72,
  fig.retina = 2,
  fig.align = "center",
  out.width = "100%",
  pngquant = "--speed=1 --quality=1-5"
)
```

# Purpose

This vignette walks through the project-detect-characterize workflow
introduced in
[vignette 1](https://ccb-hms.github.io/scDiagnostics/articles/Introduction.html)
on a real disease case study: PBMC scRNA-seq data from healthy donors
(reference) and donors with severe COVID-19 (query), from Stephenson et
al. (2021). The goal is to find out whether CD14 monocytes in the
severe-COVID query look like a distinct state relative to the healthy
reference, and if so, what distinguishes them.

```{r, message=FALSE}
library(scDiagnostics)
library(SingleCellExperiment)

set.seed(100)
```

# The data

`covid_reference_data` (healthy donors) and `covid_query_data` (severe
COVID-19 donors) are downsampled subsets of the Stephenson et al. (2021)
PBMC atlas, restricted to 5 shared cell types and a gene panel that
always includes a 25-gene interferon-response signature (Yoshida et
al.); see `?covid_reference_data` for the full processing details.

```{r, message=FALSE}
data("covid_reference_data")
data("covid_query_data")

table(covid_reference_data$author_cell_type_merged)
table(covid_query_data$azimuth_celltype_l1_merged)
```

The reference's cell type column (`author_cell_type_merged`) reflects
the original authors' annotation; the query's (`azimuth_celltype_l1_merged`)
comes from Azimuth reference mapping. Both were computed independently
of `scDiagnostics` - we are auditing an annotation transfer that has
already happened, not producing one.

# Step 1: project

`plotCellTypePCA()` projects the query onto the reference's PCA space
and compares the distributions of each cell type along the leading PCs:

```{r, fig.height=6, fig.width=10}
shared_cell_types <- c("CD14 mono", "CD4 T", "CD8 T", "B cell", "NK_16hi")

plotCellTypePCA(
    query_data = covid_query_data,
    reference_data = covid_reference_data,
    cell_types = shared_cell_types,
    query_cell_type_col = "azimuth_celltype_l1_merged",
    ref_cell_type_col = "author_cell_type_merged",
    pc_subset = 1:3)
```

# Step 2: detect

Focusing specifically on CD14 monocytes, `detectAnomaly()` builds an
Isolation Forest on the reference's PCA projection and scores how
anomalous each query cell looks relative to it:

```{r, fig.height=5, fig.width=10}
anomaly_output <- detectAnomaly(
    reference_data = covid_reference_data,
    query_data = covid_query_data,
    ref_cell_type_col = "author_cell_type_merged",
    query_cell_type_col = "azimuth_celltype_l1_merged",
    cell_types = "CD14 mono",
    pc_subset = 1:5,
    n_tree = 500)

mean(anomaly_output[["CD14 mono"]]$query_anomaly)
```

```{r, fig.height=5, fig.width=10}
plot(anomaly_output, cell_type = "CD14 mono", pc_subset = 1:3, data_type = "query")
```

In this downsampled dataset, a substantial fraction of the query's CD14
monocytes are flagged as anomalous relative to the healthy reference.
Unlike the rare/withheld-cell-type scenarios in [vignette
2](https://ccb-hms.github.io/scDiagnostics/articles/ZeiselBenchmarking.html),
this is not necessarily a small, rare subpopulation - a disease process
can plausibly shift a large fraction of a cell type's expression
profile, and that is a hypothesis worth checking directly rather than
assuming anomaly detection here means the same thing it did there.

# Step 3: characterize

`calculateGeneShifts()` tests each gene in a specified panel for a
distributional shift between reference and query, optionally comparing
only the non-anomalous reference cells against the anomalous query
cells (`anomaly_comparison = TRUE`). We focus this on the 25-gene
Yoshida et al. interferon-response signature already included in the
gene panel of both objects:

```{r, message=FALSE}
yoshida_ifn_signature <- c(
    "BST2", "CMPK2", "EIF2AK2", "EPSTI1", "HERC5", "IFI35", "IFI44L",
    "IFI6", "IFIT3", "ISG15", "LY6E", "MX1", "MX2", "OAS1", "OAS2",
    "PARP9", "PLSCR1", "SAMD9", "SAMD9L", "SP110", "STAT1", "TRIM22",
    "UBE2L6", "XAF1", "IRF7")

gene_shifts <- calculateGeneShifts(
    query_data = covid_query_data[yoshida_ifn_signature, ],
    reference_data = covid_reference_data[yoshida_ifn_signature, ],
    query_cell_type_col = "azimuth_celltype_l1_merged",
    ref_cell_type_col = "author_cell_type_merged",
    cell_types = "CD14 mono",
    pc_subset = 1:5,
    n_top_loadings = 25,
    detect_anomalies = TRUE,
    anomaly_comparison = TRUE)

head(gene_shifts$PC1[order(gene_shifts$PC1$p_adjusted), ], 10)
```

The genes with the smallest adjusted p-values here (e.g. `IFI6`, `LY6E`,
`BST2`) are all canonical interferon-stimulated genes, each with
substantially higher mean expression in the anomalous query cells than
in the non-anomalous reference cells. `plot()` visualizes this as a
heatmap of per-gene z-scores:

```{r, fig.height=6, fig.width=8}
plot(gene_shifts, cell_type = "CD14 mono", pc_subset = 1:5,
    plot_type = "heatmap", plot_by = "p_adjusted", n_genes = 15,
    show_anomalies = TRUE, pseudo_bulk = TRUE, cluster_cols = TRUE)
```

or as fold-changes relative to the reference:

```{r, fig.height=6, fig.width=8}
plot(gene_shifts, cell_type = "CD14 mono", pc_subset = 1:5,
    plot_type = "barplot", plot_by = "p_adjusted", n_genes = 15,
    show_anomalies = TRUE, pseudo_bulk = TRUE)
```

Together, these three steps give a concrete, checkable answer: yes,
CD14 monocytes in the severe-COVID query look different from the
healthy reference in PCA space, `detectAnomaly()` flags a large fraction
of them accordingly, and the genes distinguishing the flagged cells are
specifically interferon-response genes - consistent with a known
biological interferon-activated monocyte state in severe COVID-19,
rather than an artifact of the annotation transfer itself. The original
manuscript further shows this same interferon signature recovered
regardless of which of four independent annotation tools (Azimuth,
SingleR, CellTypist, scVI) produced the query labels; that
cross-tool comparison is not reproduced in this vignette.

------------------------------------------------------------------------

# R Session Info

```{r SessionInfo, echo=FALSE, message=FALSE, warning=FALSE, comment=NA}
options(width = 80)
sessionInfo()
```