## -----------------------------------------------------------------------------
#| label: setup
#| include: false
library(HTAN)
library(dplyr)
fmt <- function(x) format(x, big.mark = ",", scientific = FALSE)
pct <- function(x, digits = 1) paste0(format(round(100 * x, digits),
    nsmall = digits), "%")


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


## -----------------------------------------------------------------------------
#| label: schema
schema <- htan_tbl("schema", con) |> collect()
nrow(schema)
schema |>
    filter(Attribute == "HTAN Parent Data File ID") |>
    select(Attribute, Parent, Description) |>
    as.list()


## -----------------------------------------------------------------------------
#| label: table-list
#| echo: false
desc <- c(
    clinical_tier1_demographics =
        "Sex, race, ethnicity, vital status; one row per participant",
    clinical_tier1_diagnosis =
        "Primary diagnosis, site, stage, grade; one or more rows per participant",
    clinical_tier1_therapy = "Treatments, with days from index",
    clinical_tier1_followup = "Follow-up visits, progression, comorbidities",
    clinical_tier1_exposure = "Smoking, alcohol and other exposures",
    clinical_tier1_familyhistory = "Family history of cancer",
    clinical_tier1_moleculartest =
        "Clinical molecular tests (e.g. receptor status, mutations)",
    clinical_tier2 = "HTAN disease-agnostic clinical extensions",
    biospecimen = "Every biospecimen and its parent, type, site, preservation",
    id_provenance =
        "One row per file-parent-biospecimen link: the lineage table",
    gc_drs_map = "HTAN file ID to CRDC General Commons DRS URI",
    imaging_channel_metadata =
        "Channel (antibody/marker) lists for imaging files",
    tissueOrOrganOfOrigin_to_wholeOrgan_mapping =
        "Fine-grained site of origin to whole organ",
    phase1_longitudinal_data =
        "Per-participant timeline of biopsies, treatment and vital status",
    schema = "The Phase 1 data model: attributes, descriptions, valid values"
)
htan_tables(con) |>
    filter(name %in% names(desc)) |>
    mutate(description = desc[name]) |>
    select(table = name, rows, description) |>
    arrange(match(table, names(desc))) |>
    knitr::kable(format.args = list(big.mark = ","))


## -----------------------------------------------------------------------------
#| label: components
#| echo: false
key <- function(x) tolower(gsub("[^A-Za-z0-9]", "", x))
comp <- htan_files(con = con) |>
    count(Component, assay, level) |>
    collect()
schema |>
    mutate(key = key(Attribute)) |>
    select(key, Description) |>
    distinct(key, .keep_all = TRUE) |>
    right_join(mutate(comp, key = key(Component)), by = "key") |>
    arrange(assay, level) |>
    select(Component, level, files = n, Description) |>
    knitr::kable(format.args = list(big.mark = ","))


## -----------------------------------------------------------------------------
#| label: parse-ids
parse_htan_id <- function(id) {
    m <- regmatches(id, regexec("^(HTA[0-9]+)_([0-9]+)(?:_(.+))?$", id))
    tibble::tibble(
        id = id,
        center = vapply(m, `[`, "", 2L),
        participant = vapply(m, function(x) paste(x[2:3], collapse = "_"),
            ""),
        rest = vapply(m, `[`, "", 4L)
    )
}
parse_htan_id(c("HTA11_7862", "HTA11_7862_2000001011",
    "HTA11_7862_200000101134111", "HTA5_0000_1003255"))


## -----------------------------------------------------------------------------
#| label: prefix-center
htan_files(con = con) |>
    mutate(prefix = regexp_extract(HTAN_Data_File_ID, "^HTA[0-9]+")) |>
    count(prefix, HTAN_Center) |>
    arrange(as.integer(substr(prefix, 4L, 6L))) |>
    collect()


## -----------------------------------------------------------------------------
#| label: prefix-vs-participant
f <- htan_files(con = con) |>
    select(HTAN_Data_File_ID, HTAN_Participant_ID, n_participants) |>
    collect() |>
    mutate(id_participant = sub("^(HTA[0-9]+_[0-9]+).*$", "\\1",
        HTAN_Data_File_ID))
single <- filter(f, n_participants == 1)
mismatch <- filter(single, id_participant != HTAN_Participant_ID)
nrow(mismatch)
mismatch |>
    count(prefix = sub("_.*", "", HTAN_Data_File_ID), sort = TRUE)
head(mismatch, 3)


## -----------------------------------------------------------------------------
#| label: biospecimens
bs <- htan_tbl("biospecimen", con) |>
    select(HTAN_Biospecimen_ID, HTAN_Parent_ID, Biospecimen_Type,
        Acquisition_Method_Type, Preservation_Method, Tumor_Tissue_Type,
        Site_of_Resection_or_Biopsy) |>
    collect()
participants <- htan_tbl("clinical_tier1_demographics", con) |>
    distinct(HTAN_Participant_ID) |>
    pull()

bs |>
    mutate(parent_is = case_when(
        HTAN_Parent_ID == HTAN_Biospecimen_ID ~ "itself",
        HTAN_Parent_ID %in% participants ~ "a participant",
        HTAN_Parent_ID %in% bs$HTAN_Biospecimen_ID ~ "a biospecimen",
        .default = "an ID in neither table"
    )) |>
    count(parent_is)


## -----------------------------------------------------------------------------
#| label: walk-biospecimen
biospecimen_path <- function(id, bs) {
    parent <- setNames(bs$HTAN_Parent_ID, bs$HTAN_Biospecimen_ID)
    path <- id
    while (id %in% names(parent) && !parent[[id]] %in% path) {
        id <- parent[[id]]
        path <- c(path, id)
    }
    path
}
biospecimen_path("HTA5_1430_1022", bs)


## -----------------------------------------------------------------------------
#| label: walk-biospecimen-table
bs |>
    filter(HTAN_Biospecimen_ID %in% biospecimen_path("HTA5_1430_1022", bs)) |>
    arrange(HTAN_Biospecimen_ID) |>
    select(HTAN_Biospecimen_ID, HTAN_Parent_ID, Preservation_Method,
        Site = Site_of_Resection_or_Biopsy)


## -----------------------------------------------------------------------------
#| label: prov-columns
htan_columns("id_provenance", con) |>
    select(column, description) |>
    mutate(description = substr(description, 1, 70))


## -----------------------------------------------------------------------------
#| label: prov-size
prov <- htan_tbl("id_provenance", con)
prov |>
    summarise(
        rows = n(),
        files = n_distinct(HTAN_Data_File_ID),
        participants = n_distinct(HTAN_Participant_ID),
        biospecimens = n_distinct(HTAN_Assayed_Biospecimen_ID)
    ) |>
    collect()


## -----------------------------------------------------------------------------
#| label: prov-multiplicity
per_file <- prov |>
    group_by(HTAN_Data_File_ID) |>
    summarise(
        parents = n_distinct(HTAN_Parent_Data_File_ID),
        biospecimens = n_distinct(HTAN_Assayed_Biospecimen_ID),
        participants = n_distinct(HTAN_Participant_ID)
    ) |>
    collect()
bucket <- function(x) {
    cut(x, c(-Inf, 0, 1, 2, 10, 100, Inf),
        labels = c("0", "1", "2", "3-10", "11-100", ">100"))
}
tibble::tibble(
    count = levels(bucket(1)),
    parents = as.vector(table(bucket(per_file$parents))),
    biospecimens = as.vector(table(bucket(per_file$biospecimens))),
    participants = as.vector(table(bucket(per_file$participants)))
)


## -----------------------------------------------------------------------------
#| label: visium-file
visium_id <- "HTA11_7862_200000101134111"
htan_files(HTAN_Data_File_ID == visium_id, con = con) |>
    select(Component, Filename, File_Size, HTAN_Participant_ID,
        HTAN_Assayed_Biospecimen_ID, access) |>
    collect() |>
    glimpse()


## -----------------------------------------------------------------------------
#| label: visium-lineage
htan_lineage(visium_id, con = con) |>
    arrange(depth, related_file_id) |>
    collect() |>
    transmute(depth, related_file_id, file = basename(Filename))


## -----------------------------------------------------------------------------
#| label: visium-biospecimen
biospecimen_path("HTA11_7862_2000001011", bs)
bs |>
    filter(HTAN_Biospecimen_ID == "HTA11_7862_2000001011") |>
    select(-HTAN_Parent_ID) |>
    glimpse()


## -----------------------------------------------------------------------------
#| label: visium-diagnosis
htan_tbl("clinical_tier1_diagnosis", con) |>
    filter(HTAN_Participant_ID == "HTA11_7862") |>
    select(Primary_Diagnosis, Tissue_or_Organ_of_Origin, Age_at_Diagnosis) |>
    collect()


## -----------------------------------------------------------------------------
#| label: visium-assay-table
htan_tbl("10xvisium_spatialtranscriptomics_scRNAseq_level4_metadata",
    con) |>
    filter(HTAN_Data_File_ID == visium_id) |>
    pull(HTAN_Parent_Data_File_ID)


## -----------------------------------------------------------------------------
#| label: visium-prov-rows
prov |>
    filter(HTAN_Data_File_ID == visium_id) |>
    select(HTAN_Parent_Data_File_ID, HTAN_Assayed_Biospecimen_ID) |>
    collect()


## -----------------------------------------------------------------------------
#| label: comma-parents
count_commas <- function(table) {
    htan_tbl(table, con) |>
        summarise(
            files = n(),
            several_parents = sum(as.integer(
                grepl(",", HTAN_Parent_Data_File_ID)), na.rm = TRUE)
        ) |>
        collect() |>
        mutate(table = table, .before = 1)
}
bind_rows(lapply(c("scRNAseq_level3_metadata", "scRNAseq_level4_metadata",
    "imaging_level4_metadata",
    "10xvisium_spatialtranscriptomics_scRNAseq_level4_metadata"),
    count_commas))


## -----------------------------------------------------------------------------
#| label: descendants
htan_lineage("HTA12_251_1000", direction = "descendants", con = con) |>
    arrange(depth, related_file_id) |>
    select(depth, related_file_id, Component) |>
    collect()


## -----------------------------------------------------------------------------
#| label: imaging-lineage
htan_lineage("HTA7_1000_1003", con = con) |>
    arrange(depth) |>
    collect() |>
    transmute(depth, related_file_id, Component, file = basename(Filename))


## -----------------------------------------------------------------------------
#| label: pooled
pooled_id <- "HTA5_0000_1003255"
htan_files(HTAN_Data_File_ID == pooled_id, con = con) |>
    select(Filename, n_participants, HTAN_Participant_ID) |>
    collect() |>
    glimpse()
prov |>
    filter(HTAN_Data_File_ID == pooled_id) |>
    distinct(HTAN_Parent_Data_File_ID, HTAN_Participant_ID,
        HTAN_Assayed_Biospecimen_ID) |>
    collect()


## -----------------------------------------------------------------------------
#| label: duplicates
dup_count <- function(table, id) {
    htan_tbl(table, con) |>
        summarise(rows = n(), distinct_ids = n_distinct(.data[[id]])) |>
        collect() |>
        mutate(table = table, .before = 1)
}
bind_rows(
    dup_count("clinical_tier1_demographics", "HTAN_Participant_ID"),
    dup_count("clinical_tier1_diagnosis", "HTAN_Participant_ID"),
    dup_count("biospecimen", "HTAN_Biospecimen_ID")
)


## -----------------------------------------------------------------------------
#| label: missingness
htan_tbl("biospecimen", con) |>
    count(Tumor_Tissue_Type) |>
    arrange(desc(n)) |>
    collect()


## -----------------------------------------------------------------------------
#| label: text-numbers
htan_tbl("clinical_tier1_diagnosis", con) |>
    head(1) |>
    select(Age_at_Diagnosis, Days_to_Last_Follow_up, Year_of_Diagnosis) |>
    collect() |>
    vapply(class, "")


## -----------------------------------------------------------------------------
#| label: no-provenance
htan_files(n_participants == 0, con = con) |>
    count(Component) |>
    arrange(desc(n)) |>
    collect()


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