library(HTAN)
library(dplyr)
library(ggplot2)
con <- htan_con("r7")

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)
#> [1] 221688     11
dim(file_participants)
#> [1] 245119      7

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

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

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)

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)
#> # A tibble: 15 × 4
#>    organ                           n main_center     main_center_share
#>    <chr>                       <int> <chr>           <chr>            
#>  1 Breast                       1083 HTAN Duke       76%              
#>  2 Lung                          336 HTAN BU         54%              
#>  3 Colorectal                    194 HTAN HMS        52%              
#>  4 Not reported                  191 HTAN Vanderbilt 92%              
#>  5 Other and Ill-defined Sites    83 HTAN HMS        53%              
#>  6 Skin                           73 HTAN DFCI       60%              
#>  7 Pancreas                       64 HTAN WUSTL      95%              
#>  8 Bone Marrow                    62 HTAN CHOP       48%              
#>  9 Adrenal Gland                  34 HTAN HTAPP      85%              
#> 10 Ovary                          27 HTAN HTAPP      67%              
#> 11 Head and Neck                  25 HTAN WUSTL      84%              
#> 12 Brain                          24 HTAN HTAPP      83%              
#> 13 Cervix                         20 HTAN WUSTL      100%             
#> 14 Kidney                         15 HTAN WUSTL      87%              
#> 15 Soft Tissue                    13 HTAN HTAPP      92%

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)

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

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

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

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

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)
#> # A tibble: 12 × 2
#>    combo                                                                       n
#>    <chr>                                                                   <int>
#>  1 Bulk sequencing + H&E                                                     574
#>  2 Bulk sequencing + Single-cell sequencing                                  417
#>  3 Single-cell sequencing                                                    248
#>  4 H&E                                                                       227
#>  5 Bulk sequencing                                                           143
#>  6 Multiplexed imaging                                                       135
#>  7 Bulk sequencing + H&E + Multiplexed imaging                               117
#>  8 H&E + Multiplexed imaging                                                 104
#>  9 Bulk sequencing + Multiplexed imaging + Single-cell sequencing + Spati…    79
#> 10 Bulk sequencing + Single-cell sequencing + Spatial transcriptomics         67
#> 11 Multiplexed imaging + Spatial transcriptomics                              58
#> 12 Bulk sequencing + Multiplexed imaging + Single-cell sequencing             52

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)

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

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)

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")
#> # A tibble: 7 × 5
#>   modality                level files median largest
#>   <chr>                   <int> <int> <chr>  <chr>  
#> 1 Bulk sequencing             3  1656 3 MB   432 MB 
#> 2 Multiplexed imaging         3  4626 2 MB   8 GB   
#> 3 Multiplexed imaging         4  3539 4 MB   18 GB  
#> 4 Single-cell sequencing      3  4711 18 MB  5 GB   
#> 5 Single-cell sequencing      4  1075 389 MB 40 GB  
#> 6 Spatial transcriptomics     3  1646 294 kB 124 MB 
#> 7 Spatial transcriptomics     4   175 158 MB 365 MB

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

bios |>
    mutate(Preservation_Method = norm(Preservation_Method)) |>
    count(Preservation_Method, sort = TRUE) |>
    slice_head(n = 8) |>
    knitr::kable(format.args = list(big.mark = ","))

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

htan_tables(con) |>
    filter(grepl("tier3", name)) |>
    select(name, rows)
#> # A tibble: 5 × 2
#>   name                       rows
#>   <chr>                     <dbl>
#> 1 clinical_tier3_brain         82
#> 2 clinical_tier3_breast       129
#> 3 clinical_tier3_colorectal   304
#> 4 clinical_tier3_lung          77
#> 5 clinical_tier3_melanoma      61

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)
#> # A tibble: 10 × 3
#>    table                    field                             informative
#>    <chr>                    <chr>                                   <dbl>
#>  1 clinical_tier1_diagnosis Age_at_Diagnosis                        0.971
#>  2 clinical_tier1_diagnosis Tissue_or_Organ_of_Origin               0.915
#>  3 clinical_tier1_diagnosis Days_to_Last_Follow_up                  0.887
#>  4 clinical_tier1_diagnosis Days_to_Last_Known_Disease_Status       0.887
#>  5 clinical_tier1_diagnosis Primary_Diagnosis                       0.864
#>  6 clinical_tier1_diagnosis Site_of_Resection_or_Biopsy             0.853
#>  7 clinical_tier1_diagnosis Last_Known_Disease_Status               0.761
#>  8 clinical_tier1_diagnosis Progression_or_Recurrence               0.7  
#>  9 clinical_tier1_diagnosis Year_of_Diagnosis                       0.576
#> 10 clinical_tier1_diagnosis Tumor_Grade                             0.564

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)

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)
#> # A tibble: 3 × 2
#>   Gender           n
#>   <chr>        <int>
#> 1 Female        1734
#> 2 Male           524
#> 3 Not reported    94
count(demo_clean, Vital_Status, sort = TRUE)
#> # A tibble: 3 × 2
#>   Vital_Status     n
#>   <chr>        <int>
#> 1 Alive         1389
#> 2 Dead           483
#> 3 Not reported   480
count(demo_clean, Race, sort = TRUE)
#> # A tibble: 6 × 2
#>   Race                          n
#>   <chr>                     <int>
#> 1 White                      1611
#> 2 Not reported                347
#> 3 Black or african american   316
#> 4 Asian                        53
#> 5 Other                        19
#> 6 Not allowed to collect        6

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)
    )
#> # A tibble: 1 × 4
#>   participants missing_or_code  zero at_18_years
#>          <int>           <int> <int>       <int>
#> 1         2273              66   287         318
count(demo_clean, Age_Is_Obfuscated)
#> # A tibble: 3 × 2
#>   Age_Is_Obfuscated     n
#>   <chr>             <int>
#> 1 False               146
#> 2 True               1746
#> 3 <NA>                460
ages |>
    filter(age_days == 6573) |>
    count(HTAN_Center, sort = TRUE)
#> # A tibble: 10 × 2
#>    HTAN_Center          n
#>    <chr>            <int>
#>  1 HTAN HTAPP          82
#>  2 HTAN DFCI           61
#>  3 HTAN CHOP           48
#>  4 HTAN HMS            44
#>  5 HTAN WUSTL          37
#>  6 HTAN MSK            29
#>  7 HTAN Stanford        9
#>  8 HTAN BU              6
#>  9 HTAN TNP - TMA       1
#> 10 HTAN TNP SARDANA     1

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

long <- htan_tbl("phase1_longitudinal_data", con) |> collect()
count(long, Event_Type, sort = TRUE)
#> # A tibble: 13 × 2
#>    Event_Type                      n
#>    <chr>                       <int>
#>  1 Biopsy or Resection          1731
#>  2 Diagnosis                    1492
#>  3 Last Followup_from_followup  1321
#>  4 Vital Status                  803
#>  5 Recurrence                    571
#>  6 Recurrence_from_followup      542
#>  7 Therapy                       265
#>  8 Molecular Test                238
#>  9 Last Follow Up                217
#> 10 Last Known Disease Status     182
#> 11 Progression_from_followup     160
#> 12 Progression                    57
#> 13 Progression Free               41
biopsies <- long |>
    filter(Event_Type == "Biopsy or Resection") |>
    group_by(HTAN_Participant_ID) |>
    summarise(timepoints = n_distinct(Start))
table(timepoints = pmin(biopsies$timepoints, 4))
#> timepoints
#>    1 
#> 1731

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.3 dplyr_1.2.1   HTAN_0.99.1  
#> 
#> loaded via a namespace (and not attached):
#>  [1] utf8_1.2.6          generics_0.1.4      tidyr_1.3.2        
#>  [4] RSQLite_3.53.3      digest_0.6.39       magrittr_2.0.5     
#>  [7] evaluate_1.0.5      grid_4.6.1          RColorBrewer_1.1-3 
#> [10] fastmap_1.2.0       blob_1.3.0          jsonlite_2.0.0     
#> [13] DBI_1.3.0           purrr_1.2.2         scales_1.4.0       
#> [16] httr2_1.3.0         duckdb_1.5.6        cli_3.6.6          
#> [19] rlang_1.3.0         dbplyr_2.6.0        bit64_4.8.6        
#> [22] withr_3.0.3         cachem_1.1.0        yaml_2.3.12        
#> [25] otel_0.2.0          tools_4.6.1         memoise_2.0.1      
#> [28] filelock_1.0.3      curl_8.1.0          buildtools_1.0.0   
#> [31] vctrs_0.7.3         R6_2.6.1            BiocFileCache_3.3.0
#> [34] lifecycle_1.0.5     bit_4.6.0           pkgconfig_2.0.3    
#> [37] pillar_1.11.1       gtable_0.3.6        glue_1.8.1         
#> [40] xfun_0.62           tibble_3.3.1        tidyselect_1.2.1   
#> [43] sys_3.4.3           knitr_1.52          farver_2.1.2       
#> [46] htmltools_0.5.9     rmarkdown_2.32      maketools_1.3.2    
#> [49] labeling_0.4.3      compiler_4.6.1      S7_0.2.2