---
title: "Exploring the HTAN Phase 1 atlas"
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{Exploring the HTAN Phase 1 atlas}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---

```{r}
#| label: setup
#| include: false
fmt <- function(x) format(x, big.mark = ",", scientific = FALSE)
pct <- function(x, digits = 0) paste0(round(100 * x, digits), "%")
tb <- function(x) paste(format(round(x / 1e12, 1), nsmall = 1), "TB")
```

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

This vignette describes HTAN Phase 1 as of Release 7, the final Phase 1
release: who the participants are, where their samples came from, which
assays were run on them, at which data levels the results were released,
and how much clinical annotation comes with them. Read it before you plan
an analysis, to judge what the data can and cannot support.

Every table and figure below, and nearly every count in the text, is
computed from the Release 7 metadata snapshot. The [data model
vignette](data-model.html) explains the tables used here, and its warnings
about duplicates, comma-separated lists and inconsistent spelling apply
throughout.

We open one connection to the snapshot and use it for everything below.

```{r}
#| label: load
library(HTAN)
library(dplyr)
library(ggplot2)
con <- htan_con("r7")
```

```{r}
#| label: theme
#| include: false
theme_set(theme_minimal(base_size = 11) +
    theme(panel.grid.minor = element_blank(),
        panel.grid.major.y = element_blank(),
        plot.title.position = "plot"))
blue <- "#2a78d6"
tier_colours <- c(`open-synapse` = "#2a78d6", `open-gc` = "#eb6834",
    controlled = "#1baf7a")
```

# Building the working tables

Most of what follows works from three tables. `files` has one row per file
from `htan_files()`, plus a coarse *modality*. `file_participants` links
files to participants through `id_provenance`, one row per file and
participant, so a pooled file counts once for every participant in it.
`imaging_types` gives the imaging technology of each Level 2 image. We need
it because `Component` does not record the technology: every multiplexed
image, H&E slide and MERFISH run is `ImagingLevel2`.

```{r}
#| label: working-tables
imaging_types <- htan_tbl("imaging_level2_metadata", con) |>
    select(HTAN_Data_File_ID, Imaging_Assay_Type) |>
    collect()

modality_of <- function(assay, imaging_type) {
    case_when(
        assay %in% c("ScRNA-seq", "ScATAC-seq", "ScDNA-seq") ~
            "Single-cell sequencing",
        grepl("Visium|Xenium|GeoMx|Slide-seq", assay) |
            imaging_type == "MERFISH" | assay == "ExSeqMinimal" ~
            "Spatial transcriptomics",
        imaging_type == "H&E" ~ "H&E",
        assay %in% c("Imaging", "SRRSImaging") ~ "Multiplexed imaging",
        grepl("^Bulk|HI-C", assay) ~ "Bulk sequencing",
        assay == "ElectronMicroscopy" ~ "Electron microscopy",
        assay %in% c("MassSpectrometry", "RPPA") ~
            "Proteomics and metabolomics",
        .default = "Other"
    )
}

files <- htan_files(con = con) |>
    select(HTAN_Data_File_ID, HTAN_Center, Component, assay, level,
        extension, File_Size, access, n_participants) |>
    collect() |>
    left_join(imaging_types, by = "HTAN_Data_File_ID") |>
    mutate(modality = modality_of(assay, Imaging_Assay_Type))

file_participants <- htan_tbl("id_provenance", con) |>
    distinct(HTAN_Data_File_ID, HTAN_Participant_ID) |>
    collect() |>
    inner_join(select(files, HTAN_Data_File_ID, HTAN_Center, assay, level,
        modality, Imaging_Assay_Type), by = "HTAN_Data_File_ID")

dim(files)
dim(file_participants)
```

A Level 3 or 4 imaging file has the imaging type of the Level 2 image it
was computed from, and the join above does not trace that link. So the
modality is exact for Level 2 images and approximate for derived imaging
files, which fall under "Multiplexed imaging" unless their assay says
otherwise. Participant counts below use every level.

# The release at a glance

```{r}
#| label: glance
#| include: false
demo <- htan_tbl("clinical_tier1_demographics", con) |> collect()
bios <- htan_tbl("biospecimen", con) |> collect()
n_part_files <- n_distinct(file_participants$HTAN_Participant_ID)
n_part_demo <- n_distinct(demo$HTAN_Participant_ID)
n_bios <- n_distinct(bios$HTAN_Biospecimen_ID)
n_bios_assayed <- htan_tbl("id_provenance", con) |>
    summarise(n = n_distinct(HTAN_Assayed_Biospecimen_ID)) |>
    pull()
size_by_tier <- files |>
    group_by(access) |>
    summarise(files = n(), bytes = sum(File_Size, na.rm = TRUE))
```

Release 7 lists `r fmt(nrow(files))` files totalling
`r tb(sum(files$File_Size, na.rm = TRUE))`. They were derived from
`r fmt(n_bios_assayed)` assayed biospecimens and `r fmt(n_part_files)`
participants. The clinical demographics table has records for
`r fmt(n_part_demo)` participants and the biospecimen table lists
`r fmt(n_bios)` biospecimens. Both include participants and samples with
no released files, so to count what has data, count from `id_provenance`,
as we do here. (The DCC's 2025 description of Phase 1 counted 2,042
participants and 8,425 biospecimens, before Release 7 added data from
another 422 participants.)

By access tier:

```{r}
#| label: access-tiers
size_by_tier |>
    mutate(share_of_bytes = pct(bytes / sum(bytes)), size = tb(bytes)) |>
    select(access, files, size, share_of_bytes) |>
    knitr::kable(format.args = list(big.mark = ","))
```

Open files on Synapse, the processed data most analyses start from, are
`r pct(size_by_tier$files[size_by_tier$access == "open-synapse"] / nrow(files))`
of the files but only
`r pct(size_by_tier$bytes[size_by_tier$access == "open-synapse"] / sum(size_by_tier$bytes))`
of the bytes. Most of the volume is raw sequencing (controlled) and Level 2
images in the CRDC General Commons.

# Centers

```{r}
#| label: centers
centers <- file_participants |>
    group_by(HTAN_Center) |>
    summarise(participants = n_distinct(HTAN_Participant_ID)) |>
    left_join(
        files |>
            group_by(HTAN_Center) |>
            summarise(files = n(),
                size_tb = round(sum(File_Size, na.rm = TRUE) / 1e12, 2),
                assays = n_distinct(assay)),
        by = "HTAN_Center") |>
    arrange(desc(participants))
knitr::kable(centers, format.args = list(big.mark = ","))
```

```{r}
#| label: centers-facts
#| include: false
top <- centers$HTAN_Center[1]
dcis_top <- htan_tbl("clinical_tier1_diagnosis", con) |>
    filter(HTAN_Center == !!top,
        Primary_Diagnosis == "Ductal carcinoma in situ NOS") |>
    summarise(n = n_distinct(HTAN_Participant_ID)) |>
    pull()
most_files <- centers$HTAN_Center[which.max(centers$files)]
em_share <- mean(files$assay[files$HTAN_Center == most_files] ==
    "ElectronMicroscopy")
```

The centers differ by orders of magnitude, and which center looks largest
depends on what you count. `r most_files` has the most files
(`r fmt(max(centers$files))`), `r pct(em_share)` of them electron
microscopy, which produces one file per imaged plane or tile, but only
`r centers$participants[centers$HTAN_Center == most_files]` participants.
`r top` has by far the most participants (`r fmt(centers$participants[1])`),
`r fmt(dcis_top)` of them with a diagnosis of ductal carcinoma in situ.
File counts say more about how a center packaged its data than about how
much biology it covers. To compare scope, count participants or
biospecimens.

```{r}
#| label: fig-centers
#| fig-height: 4.5
#| fig-cap: "Participants with at least one released file, by center."
ggplot(centers, aes(participants, reorder(HTAN_Center, participants))) +
    geom_col(fill = blue, width = 0.7) +
    geom_text(aes(label = fmt(participants)), hjust = -0.15, size = 3.2,
        colour = "grey25") +
    scale_x_continuous(expand = expansion(mult = c(0, 0.12))) +
    labs(x = "Participants", y = NULL)
```

# Organs and cancer types

The diagnosis table records each participant's tissue or organ of origin
in ICD-O-3 site terms ("Upper-outer quadrant of breast", "Sigmoid colon").
The snapshot's `tissueOrOrganOfOrigin_to_wholeOrgan_mapping` table
collapses those to whole organs. Not every term in Release 7 is in the
mapping ("Upper lobe lung" and several soft-tissue sites are missing), so
we match terms case-insensitively, assign the unmapped ones by keyword and
keep anything left as it is.

```{r}
#| label: organs
organ_map <- htan_tbl("tissueOrOrganOfOrigin_to_wholeOrgan_mapping", con) |>
    collect() |>
    mutate(key = tolower(Tissue_or_Organ_of_Origin)) |>
    distinct(key, .keep_all = TRUE)
dx <- htan_tbl("clinical_tier1_diagnosis", con) |>
    select(HTAN_Participant_ID, HTAN_Center, Tissue_or_Organ_of_Origin,
        Primary_Diagnosis) |>
    collect() |>
    mutate(key = tolower(Tissue_or_Organ_of_Origin)) |>
    left_join(select(organ_map, key, Whole_Organ), by = "key") |>
    mutate(organ = case_when(
            !is.na(Whole_Organ) ~ Whole_Organ,
            grepl("lung", key) ~ "Lung",
            grepl("^skin", key) ~ "Skin",
            grepl("connective|soft tissue|paraspinal", key) ~ "Soft Tissue",
            grepl("bone", key) ~ "Bone",
            .default = Tissue_or_Organ_of_Origin),
        organ = if_else(tolower(organ) %in% c("not reported", "unknown",
            "unknown primary site"), "Not reported", organ))
with_files <- unique(file_participants$HTAN_Participant_ID)
organ_center <- dx |>
    filter(HTAN_Participant_ID %in% with_files) |>
    distinct(HTAN_Participant_ID, organ, HTAN_Center)
organs <- organ_center |>
    distinct(HTAN_Participant_ID, organ) |>
    count(organ, sort = TRUE)
main_center <- organ_center |>
    count(organ, HTAN_Center) |>
    group_by(organ) |>
    summarise(main_center = HTAN_Center[which.max(n)],
        main_center_share = pct(max(n) / sum(n)))
organs |>
    left_join(main_center, by = "organ") |>
    head(15)
```

```{r}
#| label: fig-organs
#| fig-height: 4.5
#| fig-cap: "Participants with released data, by organ of origin (top 15)."
organs |>
    slice_head(n = 15) |>
    ggplot(aes(n, reorder(organ, n))) +
    geom_col(fill = blue, width = 0.7) +
    geom_text(aes(label = fmt(n)), hjust = -0.15, size = 3.2,
        colour = "grey25") +
    scale_x_continuous(expand = expansion(mult = c(0, 0.12))) +
    labs(x = "Participants", y = NULL)
```

Breast dominates, and the `main_center` column shows why: most breast
participants come from one center's cohort. Lung, colorectal, skin
(melanoma) and pancreas follow. The ranking says more about which centers
enrolled large cohorts than about the burden of each cancer.

"Not reported" is large. Many participants, especially in the pre-cancer
atlases, have no usable organ in the diagnosis table, and their organ has
to be read from the biospecimen table (`Site_of_Resection_or_Biopsy`)
instead. A participant with several diagnosis records can appear under
more than one organ.

The most common primary diagnoses, again among participants with data:

```{r}
#| label: diagnoses
dx |>
    filter(HTAN_Participant_ID %in% with_files) |>
    distinct(HTAN_Participant_ID, Primary_Diagnosis) |>
    count(Primary_Diagnosis, sort = TRUE) |>
    slice_head(n = 20) |>
    knitr::kable(format.args = list(big.mark = ","))
```

HTAN's emphasis on transitions shows here. Ductal carcinoma in situ, a
non-invasive breast lesion, is the most common single diagnosis.
Precursor states such as familial adenomatous polyposis appear further
down the list, alongside metastatic and recurrent cancers. The `NOS`
("not otherwise specified") terms are ICD-O-3 morphology names.

# Assays and data levels

## Which assays, at which levels

File counts by assay and data level:

```{r}
#| label: assay-levels
assay_levels <- files |>
    mutate(level = if_else(is.na(level), "none", paste("Level", level))) |>
    count(modality, assay, level) |>
    tidyr::pivot_wider(names_from = level, values_from = n,
        values_fill = 0) |>
    select(modality, assay, `Level 1`, `Level 2`, `Level 3`, `Level 4`,
        none) |>
    arrange(modality, assay)
knitr::kable(assay_levels, format.args = list(big.mark = ","))
```

Sequencing assays have all four levels, but Levels 1 and 2 (reads) are
controlled, and many datasets stop at Level 3 (count matrices). Level 4
single-cell objects exist for a subset of datasets.

Imaging starts at Level 2. Raw tiles (Level 1) were rarely submitted; the
released image is the stitched OME-TIFF. Segmentation (Level 3) and cell
tables (Level 4) exist for a fraction of images.

Some assays have no levels at all. Xenium experiments are single bundles,
Visium auxiliary files and "Other Assay" files sit outside the scheme, and
ExSeq uses a minimal metadata template.

## Who ran what

Counting participants rather than files, center by modality:

```{r}
#| label: fig-heatmap
#| fig-height: 5
#| fig-width: 8
#| fig-cap: "Participants profiled with each modality, by center (any level)."
heat <- file_participants |>
    distinct(HTAN_Center, modality, HTAN_Participant_ID) |>
    count(HTAN_Center, modality) |>
    mutate(HTAN_Center = factor(HTAN_Center,
        levels = rev(centers$HTAN_Center)))
ggplot(heat, aes(modality, HTAN_Center, fill = n)) +
    geom_tile(colour = "white", linewidth = 0.6) +
    geom_text(aes(label = n, colour = log10(n) > 0.55 * log10(max(n))),
        size = 3) +
    scale_fill_gradient(low = "#e3eefb", high = "#0d3d75", trans = "log10",
        name = "Participants") +
    scale_colour_manual(values = c(`TRUE` = "white", `FALSE` = "grey15"),
        guide = "none") +
    labs(x = NULL, y = NULL) +
    theme(axis.text.x = element_text(angle = 35, hjust = 1),
        panel.grid = element_blank())
```

```{r}
#| label: heat-facts
#| include: false
centers_with <- function(m) unique(heat$HTAN_Center[heat$modality == m])
both <- intersect(centers_with("Single-cell sequencing"),
    centers_with("Multiplexed imaging"))
st_centers <- heat |>
    filter(modality == "Spatial transcriptomics") |>
    arrange(desc(n))
prot <- file_participants |>
    filter(modality == "Proteomics and metabolomics") |>
    summarise(n = n_distinct(HTAN_Participant_ID)) |>
    pull()
```

Rows are ordered by participant count. Single-cell sequencing and
multiplexed imaging are common across the network:
`r length(both)` of the `r nlevels(heat$HTAN_Center)` centers have both.
Spatial transcriptomics is less widespread: `r length(st_centers$HTAN_Center)`
centers have it, led by `r paste(head(as.character(st_centers$HTAN_Center), 3), collapse = ", ")`.
Electron microscopy is OHSU's alone, and proteomics and metabolomics cover
only `r prot` participants. The darkest tiles are bulk sequencing and H&E
at Duke. Its large breast cohort was characterised mainly by bulk RNA and
exome sequencing and histology, with single-cell sequencing for only a few
participants.

## Multiplexed imaging technologies

Level 2 images by imaging technology, with the participants and centers
behind each:

```{r}
#| label: imaging-types
imaging_participants <- file_participants |>
    filter(!is.na(Imaging_Assay_Type)) |>
    group_by(Imaging_Assay_Type) |>
    summarise(participants = n_distinct(HTAN_Participant_ID),
        centers = paste(sort(unique(sub("HTAN ", "", HTAN_Center))),
            collapse = ", ")) |>
    arrange(desc(participants))
imaging_files <- count(imaging_types, Imaging_Assay_Type, name = "images")
imaging_participants |>
    left_join(imaging_files, by = "Imaging_Assay_Type") |>
    select(Imaging_Assay_Type, images, participants, centers) |>
    knitr::kable(format.args = list(big.mark = ","))
```

H&E is the most widely shared modality. It is a common reference against
which images from different technologies can be registered and compared.
The multiplexed technologies (cyclic immunofluorescence in its CyCIF,
t-CyCIF and MxIF forms, CODEX, mIHC, MIBI, IMC and Orion) measure tens of
protein markers per cell in intact tissue. Their Level 4 cell tables are
what `read_htan()` turns into `SpatialExperiment` objects. MERFISH is
listed under imaging because it is image-based, but it measures RNA, so
we count it as spatial transcriptomics.

# Single-cell and spatial coverage

Much of HTAN's value is in participants profiled more than one way. We
give each participant the set of modalities, out of five, that they have
data for:

```{r}
#| label: combinations
combos <- file_participants |>
    filter(modality %in% c("Single-cell sequencing", "Spatial transcriptomics",
        "Multiplexed imaging", "H&E", "Bulk sequencing")) |>
    distinct(HTAN_Participant_ID, modality) |>
    group_by(HTAN_Participant_ID) |>
    summarise(combo = paste(sort(modality), collapse = " + ")) |>
    count(combo, sort = TRUE)
head(combos, 12)
```

```{r}
#| label: fig-combinations
#| fig-height: 4.5
#| fig-width: 8
#| fig-cap: >-
#|   Participants by the combination of modalities they were profiled with
#|   (12 most common combinations).
combos |>
    slice_head(n = 12) |>
    ggplot(aes(n, reorder(combo, n))) +
    geom_col(fill = blue, width = 0.7) +
    geom_text(aes(label = fmt(n)), hjust = -0.15, size = 3.2,
        colour = "grey25") +
    scale_x_continuous(expand = expansion(mult = c(0, 0.12))) +
    labs(x = "Participants", y = NULL)
```

```{r}
#| label: coverage
#| include: false
has <- function(m) {
    unique(file_participants$HTAN_Participant_ID[
        file_participants$modality == m])
}
sc <- has("Single-cell sequencing")
st <- has("Spatial transcriptomics")
mx <- has("Multiplexed imaging")
```

`r fmt(length(sc))` participants have single-cell sequencing,
`r fmt(length(st))` spatial transcriptomics and `r fmt(length(mx))`
multiplexed imaging. `r fmt(length(intersect(sc, st)))` have both
single-cell and spatial transcriptomics, and
`r fmt(length(intersect(sc, mx)))` have both single-cell sequencing and
multiplexed imaging. In these participants, cell types defined from
dissociated cells can be located in tissue.
`r fmt(sum(combos$n[!grepl("\\+", combos$combo)]))` participants
(`r pct(sum(combos$n[!grepl("\\+", combos$combo)]) / sum(combos$n))` of those
with any of these five modalities) were profiled with only one, so
cross-modal analyses draw on a much smaller set than the headline
participant count.

# What is ready to analyse

Open Level 3 and 4 files are the starting point for most reanalyses. By
modality and file type:

```{r}
#| label: analysis-ready
files |>
    filter(access == "open-synapse", level %in% c(3, 4)) |>
    count(modality, level, extension) |>
    group_by(modality, level) |>
    slice_max(n, n = 4, with_ties = FALSE) |>
    ungroup() |>
    arrange(modality, level, desc(n)) |>
    knitr::kable(format.args = list(big.mark = ","))
```

Single-cell Level 4 objects come mostly as `.rds` (R objects, often Seurat)
and `.h5ad` (AnnData), with some delimited tables; Level 3 is count
matrices in 10x, MatrixMarket and delimited formats. Imaging Level 4 is
almost entirely CSV cell tables, and Level 3 is segmentation masks.
`read_htan()` covers the `.rds`, `.h5ad`, 10x `.h5`, Space Ranger and CSV
cases.

File sizes span several orders of magnitude. That matters when you plan
downloads.

```{r}
#| label: fig-sizes
#| fig-height: 4
#| fig-cap: >-
#|   Sizes of open Level 3 and 4 files on Synapse, by modality (log scale).
files |>
    filter(access == "open-synapse", level %in% c(3, 4), File_Size > 0,
        modality %in% c("Single-cell sequencing", "Spatial transcriptomics",
            "Multiplexed imaging", "Bulk sequencing")) |>
    ggplot(aes(File_Size, modality)) +
    geom_boxplot(outlier.size = 0.4, outlier.alpha = 0.3, fill = "#e3eefb",
        colour = blue, width = 0.6) +
    facet_wrap(~ paste("Level", level)) +
    scale_x_log10(labels = scales::label_bytes(units = "auto_si")) +
    labs(x = "File size", y = NULL)
```

```{r}
#| label: size-summary
files |>
    filter(access == "open-synapse", level %in% c(3, 4), File_Size > 0,
        modality %in% c("Single-cell sequencing", "Spatial transcriptomics",
            "Multiplexed imaging", "Bulk sequencing")) |>
    group_by(modality, level) |>
    summarise(files = n(),
        median = scales::label_bytes(units = "auto_si")(median(File_Size)),
        largest = scales::label_bytes(units = "auto_si")(max(File_Size)),
        .groups = "drop")
```

Level 3 spatial transcriptomics files are small because each Space Ranger
output is split into its pieces (barcodes, features, matrix, positions).
Level 4 single-cell objects are the largest, since many combine several
samples into one object. `htan_download()` refuses downloads above 20 GB
in total unless you raise `max_bytes`.

# Biospecimens

Biospecimens by tumour tissue type, with the various missing-value codes
folded into "Not reported":

```{r}
#| label: biospecimen-types
norm <- function(x) {
    x <- trimws(x)
    x[is.na(x) | tolower(x) %in% c("not reported", "not recorded",
        "unknown", "not applicable", "")] <- "Not reported"
    paste0(toupper(substr(x, 1, 1)), substr(x, 2, nchar(x)))
}
bios |>
    mutate(Tumor_Tissue_Type = norm(Tumor_Tissue_Type)) |>
    count(Tumor_Tissue_Type, sort = TRUE) |>
    knitr::kable(format.args = list(big.mark = ","))
```

```{r}
#| label: tissue-facts
#| include: false
tt <- norm(bios$Tumor_Tissue_Type)
share_primary <- mean(tt == "Primary")
share_normal <- mean(grepl("^Normal", tt))
share_pre <- mean(grepl("^Premalignant|^Atypia", tt))
```

Primary tumour is the largest category (`r pct(share_primary)` of
biospecimens). Metastatic, premalignant and atypical tissue together are
a large share (premalignant and atypia alone are `r pct(share_pre)`), and
that is the material the study of transitions needs. Normal and
normal-adjacent tissue are `r pct(share_normal)`, so comparisons against
normal tissue are possible only for some atlases.

By preservation method:

```{r}
#| label: preservation
bios |>
    mutate(Preservation_Method = norm(Preservation_Method)) |>
    count(Preservation_Method, sort = TRUE) |>
    slice_head(n = 8) |>
    knitr::kable(format.args = list(big.mark = ","))
```

Formalin-fixed paraffin-embedded (FFPE) tissue is the majority. FFPE suits
imaging and some spatial methods; fresh or frozen tissue is needed for most
single-cell sequencing. So a biospecimen's preservation method largely
predicts which assays could have been run on it.

# Clinical variables

## What is recorded

Tier 1 clinical data are spread over seven tables. Here are their sizes
and the number of participants each covers, with the tier 2 table
alongside:

```{r}
#| label: clinical-tables
clin <- c("clinical_tier1_demographics", "clinical_tier1_diagnosis",
    "clinical_tier1_therapy", "clinical_tier1_followup",
    "clinical_tier1_exposure", "clinical_tier1_familyhistory",
    "clinical_tier1_moleculartest", "clinical_tier2")
bind_rows(lapply(clin, function(t) {
    htan_tbl(t, con) |>
        summarise(rows = n(), participants = n_distinct(HTAN_Participant_ID)) |>
        collect() |>
        mutate(table = t, .before = 1)
})) |>
    knitr::kable(format.args = list(big.mark = ","))
```

The tier 3 tables add disease-specific fields for a few hundred
participants:

```{r}
#| label: tier3
htan_tables(con) |>
    filter(grepl("tier3", name)) |>
    select(name, rows)
```

## How complete it is

A column can exist and still be mostly empty. For each diagnosis field we
compute the share of records with an informative value: not missing, and
not "Not Reported", "unknown", "Not Applicable" or a similar code.

```{r}
#| label: completeness
uninformative <- c("", "not reported", "not recorded", "unknown",
    "not applicable", "not available", "na", "unknown primary site")
completeness <- function(table) {
    x <- htan_tbl(table, con) |> collect()
    x <- x[setdiff(names(x), c("Component", "Id", "entityId", "Manifest_Id",
        "Manifest_Version", "HTAN_Center", "HTAN_Participant_ID"))]
    informative <- vapply(x, function(v) {
        v <- tolower(trimws(as.character(v)))
        mean(!is.na(v) & !v %in% uninformative)
    }, numeric(1))
    tibble::tibble(table = table, field = names(informative),
        informative = informative)
}
dx_complete <- completeness("clinical_tier1_diagnosis") |>
    arrange(desc(informative))
head(dx_complete, 10)
```

```{r}
#| label: fig-completeness
#| fig-height: 6
#| fig-cap: >-
#|   Share of diagnosis records with an informative value, for the 30
#|   best-filled fields.
dx_complete |>
    slice_head(n = 30) |>
    ggplot(aes(informative, reorder(field, informative))) +
    geom_col(fill = blue, width = 0.7) +
    scale_x_continuous(labels = scales::label_percent(),
        limits = c(0, 1), expand = expansion(mult = c(0, 0.02))) +
    labs(x = "Records with an informative value", y = NULL)
```

```{r}
#| label: completeness-summary
#| include: false
n_fields <- nrow(dx_complete)
n_half <- sum(dx_complete$informative >= 0.5)
```

Of the `r n_fields` diagnosis fields, `r n_half` are informative in at least
half the records. Primary diagnosis, site, age and the main staging fields
are usable for most participants. The long tail of disease-specific fields
(Gleason patterns, neuroblastoma classifications, Breslow thickness) is
filled only where it applies, so it is mostly empty overall. Before you
stratify an analysis by a clinical variable, compute its completeness in
the subset you will analyse rather than in the whole release.

## Demographics and outcome

Gender, vital status and race, one row per participant with data:

```{r}
#| label: demographics
demo_clean <- demo |>
    filter(HTAN_Participant_ID %in% with_files) |>
    distinct(HTAN_Participant_ID, .keep_all = TRUE) |>
    mutate(across(c(Gender, Race, Ethnicity, Vital_Status), norm))
count(demo_clean, Gender, sort = TRUE)
count(demo_clean, Vital_Status, sort = TRUE)
count(demo_clean, Race, sort = TRUE)
```

Most participants are female, because of the breast atlases. Vital status
is known for `r pct(mean(demo_clean$Vital_Status != "Not reported"))` of
participants with data, but outcome analyses also need follow-up times,
whose completeness varies by atlas; check it in your subset.

Age at diagnosis is recorded in days since birth, and two things about the
values matter before you use them. First, `0` is a placeholder, not an age
(see the [data model vignette](data-model.html#pitfalls)). Second, ages are
modified for privacy: the demographics field `Age_Is_Obfuscated` says
whether a participant's age "has been modified for compliance reasons",
and the most common value is exactly 6,573 days (18 years), which looks
like a floor applied to younger participants.

```{r}
#| label: ages
ages <- htan_tbl("clinical_tier1_diagnosis", con) |>
    select(HTAN_Participant_ID, HTAN_Center, Age_at_Diagnosis) |>
    collect() |>
    filter(HTAN_Participant_ID %in% with_files) |>
    distinct(HTAN_Participant_ID, .keep_all = TRUE) |>
    mutate(age_days = suppressWarnings(as.numeric(Age_at_Diagnosis)))
ages |>
    summarise(
        participants = n(),
        missing_or_code = sum(is.na(age_days)),
        zero = sum(age_days == 0, na.rm = TRUE),
        at_18_years = sum(age_days == 6573, na.rm = TRUE)
    )
count(demo_clean, Age_Is_Obfuscated)
ages |>
    filter(age_days == 6573) |>
    count(HTAN_Center, sort = TRUE)
```

```{r}
#| label: fig-age
#| fig-height: 3.5
#| fig-cap: >-
#|   Age at diagnosis of participants with released data (zero values
#|   removed). The spike at 18 years holds every participant reported at
#|   the privacy floor.
ages |>
    filter(age_days > 0) |>
    ggplot(aes(age_days / 365.25)) +
    geom_histogram(binwidth = 2, boundary = 0, fill = blue,
        colour = "white") +
    labs(x = "Age at diagnosis (years)", y = "Participants")
```

`r fmt(sum(ages$age_days > 0, na.rm = TRUE))` participants have a non-zero
age, with a median of
`r round(median(ages$age_days[ages$age_days > 0], na.rm = TRUE) / 365.25)`
years. The spike at 18 includes the CHOP pediatric participants, so HTAN's
age data cannot separate children from young adults. Read every age as
approximate.

# Time: longitudinal sampling

HTAN set out to sample tumours across time. The
`phase1_longitudinal_data` table puts each participant's biopsies,
diagnoses, treatments and follow-up on a single timeline (days from the
index date):

```{r}
#| label: longitudinal
long <- htan_tbl("phase1_longitudinal_data", con) |> collect()
count(long, Event_Type, sort = TRUE)
biopsies <- long |>
    filter(Event_Type == "Biopsy or Resection") |>
    group_by(HTAN_Participant_ID) |>
    summarise(timepoints = n_distinct(Start))
table(timepoints = pmin(biopsies$timepoints, 4))
```

```{r}
#| label: longitudinal-summary
#| include: false
n_multi <- sum(biopsies$timepoints >= 2)
```

`r fmt(n_multi)` participants have biopsies or resections on two or more
different days (the `4` row is four or more). They are the material for
before-and-after comparisons. The rest were sampled on a single day,
sometimes at several sites or as several sections; that supports comparisons
across space but not across time.

# Summary

- Release 7 covers `r fmt(n_part_files)` participants with data, dominated
  in number by breast and lung cohorts and in bytes by raw sequencing and
  images.
- Single-cell sequencing and multiplexed imaging are the most widely used
  modalities; spatial transcriptomics and cross-modal profiling cover a
  minority of participants.
- Most processed data are open, on Synapse, as `.rds`, `.h5ad`, 10x and CSV
  files that `read_htan()` can read.
- Clinical data are rich in fields but uneven in completeness; check the
  variables you need in the subset you study.

The [analysis vignette](analysis-workflow.html) takes one participant from
this atlas through download and analysis.

# Session information

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