## -----------------------------------------------------------------------------
#| 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")


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


## -----------------------------------------------------------------------------
#| 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")


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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))


## -----------------------------------------------------------------------------
#| 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 = ","))


## -----------------------------------------------------------------------------
#| 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 = ","))


## -----------------------------------------------------------------------------
#| 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")


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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 = ","))


## -----------------------------------------------------------------------------
#| 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 = ","))


## -----------------------------------------------------------------------------
#| 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())


## -----------------------------------------------------------------------------
#| 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()


## -----------------------------------------------------------------------------
#| 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 = ","))


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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")


## -----------------------------------------------------------------------------
#| 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 = ","))


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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")


## -----------------------------------------------------------------------------
#| 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 = ","))


## -----------------------------------------------------------------------------
#| 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))


## -----------------------------------------------------------------------------
#| 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 = ","))


## -----------------------------------------------------------------------------
#| 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 = ","))


## -----------------------------------------------------------------------------
#| label: tier3
htan_tables(con) |>
    filter(grepl("tier3", name)) |>
    select(name, rows)


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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)


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


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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)


## -----------------------------------------------------------------------------
#| 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")


## -----------------------------------------------------------------------------
#| 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))


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


## -----------------------------------------------------------------------------
#| label: session-info
sessionInfo()