---
title: "From query to analysis"
author: "Sean Davis"
date: today
format:
  html:
    toc: true
    toc-depth: 3
    fig-width: 7
    fig-height: 4.5
knitr:
  opts_chunk:
    comment: "#>"
    collapse: true
    message: false
    fig.retina: 1
vignette: >
  %\VignetteIndexEntry{From query to analysis}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---

```{r}
#| label: setup
#| include: false
## Chunks that download from Synapse run only when a token is available.
have_token <- nzchar(Sys.getenv("SYNAPSE_AUTH_TOKEN"))
fmt <- function(x) format(x, big.mark = ",", scientific = FALSE)
pct <- function(x, digits = 0) paste0(round(100 * x, digits), "%")
## An analysis result quoted in the prose: `value` when the Synapse steps
## ran, `otherwise` (text that does not depend on them) when they did not.
## `value` is only evaluated when there is a token.
result <- function(value, otherwise = "") {
    if (have_token) value else otherwise
}
code_list <- function(x) paste0("`", x, "`", collapse = ", ")
```

Package: [*HTAN*](https://bioconductor.org/packages/HTAN)
`r packageVersion("HTAN")`
([source](https://github.com/seandavi/HTAN)).

This vignette follows one HTAN dataset from a metadata query to a
biological result. We pick a small, open spatial transcriptomics section
of the colon, confirm where it came from, download it, read it into a
`SpatialExperiment`, wrap it with clinical data in a
`MultiAssayExperiment`, and run a short analysis: quality control,
normalisation, clustering and spatial plots. Each step comes with the
reason it is done that way.

It assumes you have read the [introduction](HTAN.html) and the
[data model vignette](data-model.html). The metadata steps run anywhere.
The download and every step after it need a free Synapse account and a
personal access token in `SYNAPSE_AUTH_TOKEN` (see the introduction).
Where the vignette is built without a token, as on the Bioconductor build
machines, that code is shown but not run.

```{r}
#| label: packages
#| warning: false
library(SpatialExperiment)
library(scater)
library(scrapper)
library(HTAN)
library(dplyr)
library(ggplot2)
con <- htan_con("r7")
```

The order matters: the Bioconductor packages attach functions named
`count()`, `desc()`, `slice()` and `rename()` that mask dplyr's, so dplyr
is attached last.

# Choosing a dataset

## The question

The Vanderbilt atlas (HTA11) studied colorectal precancer and cancer, and
a section of colon makes a good first spatial dataset. The colonic mucosa
has an obvious architecture (epithelial crypts, the stroma of the lamina
propria, immune cells). If the analysis is right, its expression clusters
should match recognisable tissue compartments, and agree with what the
metadata say the tissue is.

## Finding candidate files

10x Visium measures the transcriptome in a grid of 55 µm spots across a
tissue section. Level 4 Visium files are processed objects, ready to
analyse. These are the open ones, by center and format:

```{r}
#| label: visium-candidates
visium <- htan_files(
    assay == "10xVisiumSpatialTranscriptomics-RNA-seq", level == 4,
    access == "open-synapse", con = con)
visium |>
    count(HTAN_Center, extension) |>
    arrange(desc(n)) |>
    collect()
```

Three things guide the choice. The first is format: `.h5ad` (AnnData)
files follow a common layout that anndataR reads reliably, whereas `.rds`
files hold whatever R object the center saved. The second is provenance.
We want a file from a single participant, measured on a biospecimen whose
tissue type and site are recorded, so the result can be read against
them. The third is size, small enough to download in seconds. The
biospecimen fields come from the `biospecimen` table, joined on the
assayed biospecimen:

```{r}
#| label: candidates
biospecimens <- htan_tbl("biospecimen", con) |>
    select(HTAN_Assayed_Biospecimen_ID = HTAN_Biospecimen_ID,
        Tumor_Tissue_Type, Site_of_Resection_or_Biopsy)
candidates <- visium |>
    filter(extension == ".h5ad", n_participants == 1) |>
    inner_join(biospecimens, by = "HTAN_Assayed_Biospecimen_ID") |>
    filter(Tumor_Tissue_Type != "Not Otherwise Specified",
        Site_of_Resection_or_Biopsy != "Not Reported") |>
    arrange(File_Size) |>
    select(HTAN_Data_File_ID, Filename, File_Size, Tumor_Tissue_Type,
        Site_of_Resection_or_Biopsy) |>
    collect()
nrow(candidates)
candidates |>
    transmute(file = basename(Filename), MB = round(File_Size / 1e6),
        Tumor_Tissue_Type, site = Site_of_Resection_or_Biopsy) |>
    head(10)
```

One property that matters is missing from the metadata: the version of
AnnData that wrote each file. Files written before AnnData 0.8 (released
in 2022) use an older on-disk layout that anndataR does not read.
`read_htan()` reads those with a fallback that keeps the matrices, `obs`,
`var` and `obsm`, but drops `uns`. Every file in this list smaller than
`7794_AS_3_unfiltered.h5ad` is in the older layout (we checked by opening
each one). So we use `7794_AS_3`, the smallest file in the current
layout, which anndataR reads in full:

```{r}
#| label: pick
visium_id <- "HTA11_99999933430_8422134111"
stopifnot(visium_id %in% candidates$HTAN_Data_File_ID)
file <- htan_files(HTAN_Data_File_ID == visium_id, con = con) |> collect()
glimpse(select(file, HTAN_Data_File_ID, Filename, File_Size, md5,
    entityId, HTAN_Participant_ID, HTAN_Assayed_Biospecimen_ID, access))
```

The `entityId` is the file's Synapse ID. `access` is `open-synapse`, so
you need a token but not a data access request.

## What the metadata say about it

The assay table holds the processing details:

```{r}
#| label: workflow
workflow <- htan_tbl(
    "10xvisium_spatialtranscriptomics_scRNAseq_level4_metadata", con) |>
    filter(HTAN_Data_File_ID == visium_id) |>
    select(Workflow_Version, Visium_Workflow_Type,
        Visium_Workflow_Parameters_Description, Run_ID, Data_Release) |>
    collect()
glimpse(workflow)
```

The file was made with Space Ranger against the GRCh38-2020-A reference,
with the human whole-transcriptome probe set, so this is the probe-based
Visium chemistry. The record disagrees with itself about the Space Ranger
version (`Workflow_Version` is `r workflow$Workflow_Version`, but
`Visium_Workflow_Type` says `r workflow$Visium_Workflow_Type`). These
fields were filled in by hand at submission, so read them sceptically.
"Unfiltered" in the file name suggests the file holds every spot on the
capture area, including those off the tissue; if so, the analysis has to
remove them.

# Tracing lineage

Before using a processed file, check what it was made from:

```{r}
#| label: lineage
lineage <- htan_lineage(visium_id, con = con) |>
    arrange(depth, related_file_id) |>
    collect()
lineage |>
    transmute(depth, Component = sub("10xVisiumSpatialTranscriptomics-",
        "", Component), file = basename(Filename))
```

The chain runs from FASTQ (Level 1, controlled) through aligned BAM
(Level 2, controlled) and the Space Ranger output (Level 3, open) to this
object (Level 4). All were measured on one biospecimen. Here are its
record and the participant's:

```{r}
#| label: biospecimen
biospecimen <- htan_tbl("biospecimen", con) |>
    filter(HTAN_Biospecimen_ID == !!file$HTAN_Assayed_Biospecimen_ID) |>
    select(HTAN_Biospecimen_ID, HTAN_Parent_ID, Acquisition_Method_Type,
        Preservation_Method, Tumor_Tissue_Type, Site_of_Resection_or_Biopsy,
        Histologic_Morphology_Code) |>
    collect()
glimpse(biospecimen)

htan_tbl("clinical_tier1_demographics", con) |>
    filter(HTAN_Participant_ID == !!file$HTAN_Participant_ID) |>
    select(HTAN_Participant_ID, Gender, Race, Vital_Status) |>
    collect()

diagnosis <- htan_tbl("clinical_tier1_diagnosis", con) |>
    filter(HTAN_Participant_ID == !!file$HTAN_Participant_ID) |>
    select(Primary_Diagnosis, Morphology, Tissue_or_Organ_of_Origin,
        Tumor_Grade) |>
    collect()
diagnosis
```

The section is a `r tolower(biospecimen$Acquisition_Method_Type)` of the
`r tolower(biospecimen$Site_of_Resection_or_Biopsy)`, with tissue type
"`r biospecimen$Tumor_Tissue_Type`", HTAN's term for tissue from a primary
tumour (as opposed to premalignant, normal or metastatic tissue). Nothing
more specific is recorded: the morphology code is
`r biospecimen$Histologic_Morphology_Code` (not reported), and the
participant's primary diagnosis is "`r diagnosis$Primary_Diagnosis`".
So the metadata do not say whether the lesion is an adenoma or a
carcinoma, or of what type, and we have to read the analysis without
that. The participant has other data too, which could extend the
analysis later:

```{r}
#| label: other-data
htan_files(participants = file$HTAN_Participant_ID, con = con) |>
    count(Component, access) |>
    arrange(Component) |>
    collect()
```

# Downloading

`htan_download()` checks the access tier and total size, asks Synapse for
a short-lived download URL, fetches the file, and verifies its md5 against
the metadata. The file is kept in the package's download cache (a
BiocFileCache under `tools::R_user_dir("HTAN", "cache")`), so a second
call finds the verified copy instead of downloading it again. To get a
copy in a directory of your choice as well, pass `dest`.

```{r}
#| label: download
#| eval: !expr have_token
dl <- htan_download(file, con = con)
dl
```

# Reading into a SpatialExperiment

`read_htan()` reads `.h5ad` files with anndataR. When the AnnData object
has a `spatial` embedding (the layout scanpy writes for Visium), the result
is a `SpatialExperiment` with the spot positions as `spatialCoords()`:

```{r}
#| label: read
#| eval: !expr have_token
spe <- read_htan(dl)
spe
```

Level 4 objects differ between centers, so look before you compute. We
need to know which assay holds raw counts, whether there are spot
coordinates, and which spots are under tissue.

```{r}
#| label: inspect
#| eval: !expr have_token
assayNames(spe)
names(colData(spe))
names(rowData(spe))
head(rownames(spe))
is(spe, "SpatialExperiment")
```

Marker genes are easier to read as symbols. If the row names are Ensembl
IDs, we switch to a symbol column from `rowData()` (scanpy and Space
Ranger name it differently), keeping the names unique. Some probe-set
features are named `DEPRECATED_` followed by an Ensembl ID. We count them
and leave them as they are.

```{r}
#| label: symbols
#| eval: !expr have_token
if (mean(grepl("^ENSG", rownames(spe))) > 0.5) {
    symbol_col <- intersect(c("gene_symbols", "gene_name", "feature_name",
        "symbol", "gene_ids"), names(rowData(spe)))[1]
    rownames(spe) <- make.unique(as.character(rowData(spe)[[symbol_col]]))
}
head(rownames(spe))
sum(grepl("^DEPRECATED_", rownames(spe)))
```

anndataR maps AnnData's main matrix `X` to an assay named `X`. Whether `X`
holds raw counts or normalised values is up to the center, so we test it.
Raw counts are whole numbers.

```{r}
#| label: counts
#| eval: !expr have_token
x <- assay(spe, "X")
sample_values <- as.numeric(x[, seq_len(min(200, ncol(x)))])
all(sample_values == round(sample_values))
assay(spe, "counts") <- x
assay(spe, "X") <- NULL
```

If this were `FALSE`, the counts would have to come from a layer of the
AnnData object or from the Level 3 Space Ranger matrix listed in the
lineage.

If the object has no spatial embedding, the spot positions are in the
Level 3 `tissue_positions.csv` from the lineage table, which can be
downloaded and joined by barcode. This chunk does that only when it has
to:

```{r}
#| label: coordinates
#| eval: !expr have_token
if (!is(spe, "SpatialExperiment")) {
    pos_id <- lineage$related_file_id[
        grepl("tissue_positions", lineage$Filename)]
    pos_dl <- htan_download(pos_id, con = con)
    ## Space Ranger 2.0 writes a header row; 1.x did not.
    has_header <- startsWith(readLines(pos_dl$path, n = 1), "barcode")
    pos <- read.csv(pos_dl$path, header = has_header)
    names(pos) <- c("barcode", "in_tissue", "array_row", "array_col",
        "pxl_row", "pxl_col")
    pos <- pos[match(colnames(spe), pos$barcode), ]
    colData(spe)$in_tissue <- pos$in_tissue
    spe <- as(spe, "SpatialExperiment")
    spatialCoords(spe) <- as.matrix(pos[, c("pxl_col", "pxl_row")])
}
dim(spatialCoords(spe))
```

Finally, we keep only spots under tissue. Spots outside it capture
ambient RNA and would form a spurious cluster of their own:

```{r}
#| label: in-tissue
#| eval: !expr have_token
if ("in_tissue" %in% names(colData(spe))) {
    print(table(in_tissue = spe$in_tissue))
    spe <- spe[, spe$in_tissue == 1]
}
dim(spe)
```

# Adding clinical context: a MultiAssayExperiment

`htan_experiment()` builds a `MultiAssayExperiment` from objects named by
their `HTAN_Data_File_ID`. It looks up each file's participant in
`id_provenance`, maps every spot to that participant in the `sampleMap`,
and fills `colData` with the participant's demographics and diagnosis.
Biospecimen annotation is per sample rather than per participant, so it is
added to the spots' own `colData` first:

```{r}
#| label: mae
#| eval: !expr have_token
for (field in c("HTAN_Biospecimen_ID", "Tumor_Tissue_Type",
    "Site_of_Resection_or_Biopsy", "Histologic_Morphology_Code")) {
    colData(spe)[[field]] <- biospecimen[[field]]
}
mae <- htan_experiment(setNames(list(spe), visium_id), con = con)
mae
MultiAssayExperiment::colData(mae)[, c("Gender", "Race", "Vital_Status",
    "Primary_Diagnosis")]
head(MultiAssayExperiment::sampleMap(mae))
```

With one file the container is small. It pays off once more data go
in, such as the participant's single-cell data or other participants'
Visium sections: every spot and cell stays linked to the right clinical
record. For a file pooled from several participants, `htan_experiment()`
needs the `colData` column that names each cell's participant
(`participant_col`); see the [data model vignette](data-model.html#a-file-from-several-participants).

The analysis below works on the `SpatialExperiment` inside it:

```{r}
#| label: get-experiment
#| eval: !expr have_token
spe <- MultiAssayExperiment::experiments(mae)[[1]]
```

# Analysis

The steps follow the standard Bioconductor workflow for spot-based
spatial transcriptomics, as in the book *Orchestrating Spatial Transcriptomics
Analysis with Bioconductor* ([OSTA](https://bioconductor.org/books/release/OSTA/))
and its single-cell counterpart
([scrapbook](https://bioconductor.org/books/release/scrapbook/)). Like
both books, we compute with
[scrapper](https://bioconductor.org/packages/scrapper/), which implements
the methods of scran and scuttle in C++ libraries (libscran) shared with
the Python and JavaScript ports, and plot with scater. Each scrapper
function ending in `.se` takes a `SummarizedExperiment` (here the
`SpatialExperiment`) and returns it with its results added.

## Quality control

The per-spot metrics are total counts (library size) and the number of
genes detected. In most RNA data a third metric is the share of counts
from the mitochondrially encoded genes (`MT-` genes), which rises in
damaged cells. Poly(A) capture picks up whatever transcripts are present,
but a probe-based assay measures only the genes its panel targets, and
the panel used for this section has no probes for the `MT-` genes. So the
code looks for them and adds the mitochondrial share only if it finds
any. (Despite the name, genes such as `MTRNR2L12` are nuclear, and the
`^MT-` pattern rightly skips them.)

```{r}
#| label: qc
#| eval: !expr have_token
is_mito <- grepl("^MT-", rownames(spe), ignore.case = TRUE)
sum(is_mito)
qc_subsets <- if (any(is_mito)) list(mito = is_mito) else list()
spe <- quickRnaQc.se(spe, subsets = qc_subsets)
summary(spe$sum)
summary(spe$detected)
if (any(is_mito)) summary(spe$subset.proportion.mito)
```

Plotted on the tissue, the metrics show whether low-quality spots are
scattered (technical) or form a region (biological). Image row
coordinates increase downwards, hence the reversed y axis.

```{r}
#| label: fig-qc
#| eval: !expr have_token
#| fig-height: 4
#| fig-cap: "Library size of each spot under tissue."
spot_coords <- function(spe) {
    xy <- as.data.frame(spatialCoords(spe))
    names(xy) <- c("x", "y")
    xy
}
spot_plot <- function(spe, values, name) {
    ggplot(cbind(spot_coords(spe), value = values),
        aes(x, y, colour = value)) +
        geom_point(size = 1) +
        scale_y_reverse() +
        coord_fixed() +
        labs(x = NULL, y = NULL, colour = name) +
        theme_void()
}
spot_plot(spe, log10(spe$sum + 1), "log10\ncounts") +
    scale_colour_viridis_c()
```

Outlier spots are flagged by their distance from the median in median
absolute deviations (MADs), not by fixed thresholds, because a sensible
threshold depends on the tissue, the probe panel and the sequencing
depth. `quickRnaQc.se()` has already done this. It set each threshold
three MADs from the median, testing library size and detected genes on
the low side only (and the mitochondrial share, when there is one, on the
high side), and stored the verdict in `keep`. The thresholds are in the
object's metadata:

```{r}
#| label: qc-filter
#| eval: !expr have_token
thresholds <- metadata(spe)$qc$thresholds
unlist(thresholds)
n_before <- ncol(spe)
n_discard <- sum(!spe$keep)
c(low_sum = sum(spe$sum < thresholds$sum),
    low_detected = sum(spe$detected < thresholds$detected),
    discard = n_discard)
spe <- spe[, spe$keep]
dim(spe)
```

A large discarded fraction, or discarded spots forming a region, would
call for a closer look before dropping them`r result(paste0("; here ",
fmt(n_discard), " of the ", fmt(n_before), " spots ",
if (n_discard == 1) "is" else "are", " discarded"))`.

## Normalisation

Spots differ in library size for technical reasons (capture efficiency,
tissue thickness) and biological ones (cell density). The usual choice
for Visium is library-size normalisation followed by a log transform. The
pooling-based size factors common in single-cell work assume that most
genes are not differentially expressed between pools of cells, an
assumption that is less safe for spots mixing several cell types.
`normalizeRnaCounts.se()` uses the library sizes as size factors by
default, scaled to average 1, and adds the `logcounts` assay.

```{r}
#| label: normalise
#| eval: !expr have_token
spe <- normalizeRnaCounts.se(spe)
summary(spe$sizeFactor)
```

## Feature selection and dimensionality reduction

Highly variable genes carry the structure; the rest add noise. We model
the variance of each gene's log-expression against its mean and mark the
2,000 genes with the largest residual above that trend in
`rowData(spe)$hvg`. Mitochondrial genes, where there are any, are then
dropped from the set so that tissue damage does not drive the
clustering. This panel has none, so here that step does nothing.

```{r}
#| label: hvgs
#| eval: !expr have_token
spe <- chooseRnaHvgs.se(spe, top = 2000)
if (any(is_mito)) {
    rowData(spe)$hvg <- rowData(spe)$hvg & !is_mito
}
sum(rowData(spe)$hvg)
hvg_stats <- rowData(spe)[rowData(spe)$hvg, ]
head(rownames(hvg_stats)[order(hvg_stats$residuals, decreasing = TRUE)], 20)
```

```{r}
#| label: pca
#| eval: !expr have_token
spe <- runPca.se(spe, features = rowData(spe)$hvg, number = 30)
pca <- metadata(spe)$PCA
percent_var <- 100 * pca$variance.explained / pca$total.variance
round(head(percent_var, 10), 1)
```

We keep thirty components, a common default and a generous one for a
tissue with a handful of domains. No single component explains much of
the variance`r result(paste0(" (the first ", round(percent_var[1], 1),
"%, the second ", round(percent_var[2], 1), "%)"))`, which is usual for
spots that each average several cells.

## Clustering

Spots are clustered on the principal components with a shared
nearest-neighbour graph (Jaccard-weighted) and the Leiden algorithm at
resolution 0.5, the settings OSTA uses for its colorectal Visium workflow.
The clustering sees expression alone, which turns the spatial plot that
follows into a test: if clusters computed without coordinates form
coherent regions, they reflect real tissue structure. Spatially aware
methods (for example BayesSpace or Banksy) smooth over neighbours and
would give cleaner domains, but we would lose that independent check.

```{r}
#| label: cluster
#| eval: !expr have_token
spe <- clusterGraph.se(spe, method = "leiden", resolution = 0.5,
    more.build.args = list(weight.scheme = "jaccard"),
    output.name = "cluster")
table(spe$cluster)
```

```{r}
#| label: fig-clusters
#| eval: !expr have_token
#| fig-height: 4.5
#| fig-cap: "Expression clusters on the tissue section."
spot_plot(spe, spe$cluster, "Cluster") +
    guides(colour = guide_legend(override.aes = list(size = 3)))
```

```{r}
#| label: fig-umap
#| eval: !expr have_token
#| fig-height: 4
#| fig-cap: "UMAP of spots, coloured by cluster."
spe <- runUmap.se(spe)
plotReducedDim(spe, "UMAP", colour_by = "cluster")
```

## What the clusters are

Marker genes tell us what each cluster is. `scoreMarkers.se()` compares
each cluster with every other one and summarises the pairwise effect
sizes for each gene: Cohen's d (the difference in mean log-expression in
standard deviations), the AUC (the probability that a spot in the cluster
has higher expression than a spot in the other), and the differences in
mean expression and in the proportion of spots where the gene is
detected. Genes are ranked by Cohen's d averaged over the comparisons.
Here effect sizes are a better summary than p-values. The clusters were
defined from these same expression values, so a test for differences
between them is circular, and its p-values would not mean what they
usually do.

```{r}
#| label: markers
#| eval: !expr have_token
markers <- scoreMarkers.se(spe, groups = spe$cluster)
top_markers <- vapply(markers, function(m) {
    paste(head(rownames(m), 8), collapse = ", ")
}, "")
tibble::tibble(cluster = names(top_markers), top_markers = top_markers)
```

Some clusters are led by members of multi-copy gene families, such as
`SERF1A` and `SERF1B`, `FAM153B`, the `USP17L` genes and `MTRNR2L12`,
rather than by genes that mark a cell type. Each of these families has
several near-identical copies in the genome, so a probe, or a read, may
not be specific to one of them, and the `MTRNR2L` genes are also close
copies of the mitochondrial rRNA gene `MT-RNR2`. Their signal may
therefore be partly technical (this is unverified). The pattern is
easier to see in a per-cluster summary of the number of spots, their
median depth, the median share of each spot's counts that goes to these
families, and the mean log-expression of an epithelial and a stromal
marker:

```{r}
#| label: cluster-summary
#| eval: !expr have_token
multicopy <- grepl("^(SERF1|USP17L|FAM153|MTRNR2L|TAF11L|GAGE|CT45A)",
    rownames(spe))
sum(multicopy)
clusters <- factor(spe$cluster)
per_cluster <- function(v, f = median) as.vector(tapply(v, clusters, f))
cluster_summary <- tibble::tibble(
    cluster = levels(clusters),
    spots = per_cluster(spe$sum, length),
    median_counts = per_cluster(spe$sum),
    median_genes = per_cluster(spe$detected),
    multicopy_share = round(per_cluster(
        colSums(counts(spe)[multicopy, ]) / spe$sum), 2),
    EPCAM = round(per_cluster(logcounts(spe)["EPCAM", ], mean), 2),
    COL1A1 = round(per_cluster(logcounts(spe)["COL1A1", ], mean), 2)
)
cluster_summary
epi <- cluster_summary$cluster[which.max(cluster_summary$EPCAM)]
stroma <- cluster_summary$cluster[which.max(cluster_summary$COL1A1)]
shallow <- cluster_summary$cluster[which.min(cluster_summary$median_counts)]
deepest <- cluster_summary$cluster[which.max(cluster_summary$median_counts)]
```

The top of the ranking for the most epithelial cluster (highest mean
`EPCAM`) shows the statistics behind it: the mean log-expression and
detected proportion in the cluster, then the mean Cohen's d and AUC
against the other clusters.

```{r}
#| label: markers-epithelial
#| eval: !expr have_token
previewMarkers(markers[[epi]], rows = 8, columns = c("mean", "detected",
    d = "cohens.d.mean", auc = "auc.mean"))
```

To pin the clusters down, we plot known markers of the colonic mucosa:
`EPCAM` for epithelium, `MUC2` for goblet cells and their mucin, `CA1` for
the absorptive colonocytes that line normal colonic epithelium, `COL1A1`
for stroma and `PTPRC` for immune cells. `SERF1A` stands in for the
multi-copy families. Their expression levels differ widely, so each gene
is coloured relative to its own maximum:

```{r}
#| label: fig-genes
#| eval: !expr have_token
#| fig-height: 6
#| fig-width: 8
#| fig-cap: "Marker genes on the tissue, each relative to its maximum."
genes <- intersect(c("EPCAM", "MUC2", "CA1", "COL1A1", "PTPRC", "SERF1A"),
    rownames(spe))
long <- do.call(rbind, lapply(genes, function(g) {
    v <- logcounts(spe)[g, ]
    cbind(spot_coords(spe), gene = g, relative = v / max(v))
}))
ggplot(long, aes(x, y, colour = relative)) +
    geom_point(size = 0.6) +
    scale_colour_viridis_c(labels = scales::percent) +
    scale_y_reverse() +
    coord_fixed() +
    facet_wrap(~ gene) +
    labs(x = NULL, y = NULL, colour = "Share of\nmaximum") +
    theme_void()
```

```{r}
#| label: fig-cluster-genes
#| eval: !expr have_token
#| fig-height: 4
#| fig-cap: "Marker expression by cluster."
plotDots(spe, features = genes, group = "cluster")
```

```{r}
#| label: interpretation
#| eval: !expr have_token
#| include: false
detected_in <- function(g) mean(counts(spe)[g, ] > 0)
colonocyte <- intersect(c("CA1", "FABP1", "SLC26A3", "AQP8"), rownames(spe))
colonocyte_max <- max(vapply(colonocyte, detected_in, 0))
cluster_n <- function(k) paste("cluster", k)
top_genes <- function(k, n = 6) code_list(head(rownames(markers[[k]]), n))
share <- function(k) {
    pct(cluster_summary$multicopy_share[cluster_summary$cluster == k])
}
```

Read together, the plots and the summary suggest three kinds of spot:

- Epithelium. The cluster with the highest `EPCAM`
  `r result(paste0("(", cluster_n(epi), ") "))`has epithelial and
  secretory genes among its top markers`r result(paste0(" (",
  top_genes(epi, 8), ")", if (epi == deepest) ", and the deepest spots"))`.
- Stroma. The cluster with the highest `COL1A1`
  `r result(paste0("(", cluster_n(stroma), ") "))`is led by collagens and
  other stromal and extracellular-matrix genes`r result(paste0(" (",
  top_genes(stroma), ")"))`.
- Shallow spots. The cluster with the fewest counts per spot
  `r result(paste0("(", cluster_n(shallow), ") "))`is led by the multi-copy
  families`r result(paste0(" (", top_genes(shallow), ")"))`, which take
  `r result(share(shallow), "a large share")` of a typical spot's counts
  there. These spots lie along the edges of the section, where the
  library-size map is darkest, and the cluster may reflect capture depth
  more than a cell type.

The section as a whole does not look like normal colonic mucosa. `MUC2`
is detected throughout`r result(paste0(" (in ",
pct(detected_in("MUC2")), " of spots)"))`, while the markers of absorptive
colonocytes (`r code_list(c("CA1", "FABP1", "SLC26A3", "AQP8"))`) are
almost absent`r result(paste0(" (each detected in at most ",
pct(colonocyte_max), " of spots)"))`, and so is the immune marker
`PTPRC``r result(paste0(" (", pct(detected_in("PTPRC")), ")"))`. That
would fit a mucin-rich lesion, but the histology is not recorded, so it
remains a hypothesis: the H&E image of the section (an auxiliary file of a
few gigabytes) is the place to check it.

So the expression clusters do form regions on the tissue, which suggests
real structure. But they separate spots by depth and by a few gene
families as well as by tissue compartment. Two next steps would sharpen
them: drop the multi-copy families from the highly variable genes and
recluster, and compare the regions with the H&E image.

# Where to go from here

Other Vanderbilt Visium files cover more colon lesions (`candidates`
above). You can read several, give each its file ID, and combine them in
one `MultiAssayExperiment`, then integrate them with batch correction
(for example scrapper's `correctMnn.se()`) before clustering them
together.

The participant also has single-cell RNA-seq (Level 3 count matrices).
Cell types from those can be used to deconvolve the spots, for example
with
[SPOTlight](https://bioconductor.org/packages/release/bioc/html/SPOTlight.html).

For spatially smoothed domains there are
[BayesSpace](https://bioconductor.org/packages/release/bioc/html/BayesSpace.html)
and [Banksy](https://bioconductor.org/packages/release/bioc/html/Banksy.html);
for spatially variable genes,
[nnSVG](https://bioconductor.org/packages/release/bioc/html/nnSVG.html).

Beyond this participant, `htan_files()` and the clinical tables find
comparable sections from other atlases, for example premalignant versus
invasive tissue.

# Notes on reproducibility

The metadata release is pinned (`htan_con("r7")`), and the file is
identified by its HTAN ID and checked against its md5. The random steps
(PCA, clustering, UMAP) use scrapper's fixed default seeds, so on a given
platform and set of package versions they give the same result on every
run without `set.seed()`. Across versions and platforms the details can
change: builds of this vignette with different scrapper releases and
numerical libraries have given five or six clusters. That is why every
result quoted in the text above is computed from the output rather than
written in by hand. The software versions are listed below.

The file's content on Synapse is not pinned. If the center replaced the
file, the md5 check would fail rather than silently give a different
result.

# Session information

```{r}
#| label: session-info
sessionInfo()
```