knitr::kable(search_gpmap("SLC44A")[, c("name", "sample_size", "num_study_extractions", "num_coloc_groups", "call")])

knitr::kable(search_gpmap("22:37042914"))

proxy_variants <- search_gpmap("rs17078078")
knitr::kable(proxy_variants)

hemoglobin_trait <- trait(931)
names(hemoglobin_trait)
#> [1] "trait"             "coloc_groups"      "study_extractions"
nrow(hemoglobin_trait$coloc_groups)
#> [1] 27999
nrow(hemoglobin_trait$study_extractions)
#> [1] 1342

slc44a1_coloc_pairs <- gene("SLC44A1", include_coloc_pairs = TRUE)
nrow(slc44a1_coloc_pairs$coloc_pairs)
#> [1] 1096
knitr::kable(head(slc44a1_coloc_pairs$coloc_pairs[, c("study_extraction_a_id", "study_extraction_b_id", "h4")]))

study_extractions_tbl <- dplyr::bind_rows(slc44a1_no_trans$study_extractions)
study_extractions_subset <- dplyr::select(study_extractions_tbl, id, trait_name, min_p)

hydrated_coloc_pairs <- slc44a1_coloc_pairs$coloc_pairs |>
  dplyr::left_join(slc44a1_coloc_pairs$study_extractions, by = c("study_extraction_a_id" = "id")) |>
  dplyr::rename(trait_a = trait_name) |>
  dplyr::left_join(slc44a1_coloc_pairs$study_extractions, by = c("study_extraction_b_id" = "id")) |>
  dplyr::rename(trait_b = trait_name)

knitr::kable(head(hydrated_coloc_pairs[, c("trait_a", "trait_b", "h4")]))

ad_search <- search_gpmap("alzheimer")
ad_trait_ids <- ad_search$type_id

ad_coloc <- dplyr::bind_rows(lapply(ad_trait_ids, function(tid) {
  tryCatch(trait(tid)$coloc_groups, error = function(e) NULL)
})) |>
  dplyr::filter(!grepl("methyl", .data$data_type, ignore.case = TRUE))

ad_genes <- ad_coloc |>
  dplyr::filter(!is.na(.data$gene_id)) |>
  dplyr::distinct(.data$gene_id, .data$gene)

nrow(ad_genes)
#> [1] 228

all_pleiotropy <- get_all_gene_pleiotropies()

ad_genes_ranked <- ad_genes |>
  dplyr::inner_join(
    all_pleiotropy[, c("gene_id", "distinct_trait_categories", "distinct_protein_coding_genes")],
    by = "gene_id"
  ) |>
  dplyr::arrange(.data$distinct_trait_categories)

top_genes <- head(ad_genes_ranked, 40)

gene_evidence <- lapply(top_genes$gene, function(g) {
  gi <- tryCatch(gene(g, include_associations = TRUE), error = function(e) NULL)
  if (is.null(gi)) return(NULL)

  cg <- gi$coloc_groups
  rr <- gi$rare_results

  ad_cg <- if (!is.null(cg) && nrow(cg) > 0) {
    cg[cg$trait_id %in% ad_trait_ids, ]
  } else {
    data.frame()
  }

  ad_rr <- if (!is.null(rr) && nrow(rr) > 0) {
    rr[rr$trait_id %in% ad_trait_ids, ]
  } else {
    data.frame()
  }

  n_tissues <- if (!is.null(cg) && nrow(cg) > 0 && "tissue" %in% names(cg)) {
    dplyr::n_distinct(cg$tissue[!is.na(cg$tissue)])
  } else {
    0L
  }

  data.frame(
    gene = g,
    n_coloc_groups = dplyr::n_distinct(ad_cg$coloc_group_id),
    n_tissues = n_tissues,
    n_rare_groups = if (nrow(ad_rr) > 0) dplyr::n_distinct(ad_rr$rare_result_group_id) else 0L
  )
}) |> dplyr::bind_rows()

ad_summary <- ad_genes_ranked |>
  dplyr::filter(.data$gene %in% gene_evidence$gene) |>
  dplyr::left_join(gene_evidence, by = "gene") |>
  dplyr::select("gene", "distinct_trait_categories", "distinct_protein_coding_genes",
                "n_tissues", "n_coloc_groups", "n_rare_groups") |>
  dplyr::arrange(.data$distinct_trait_categories, .data$gene)

ad_summary_long <- ad_summary |>
  tidyr::pivot_longer(
    cols = c("distinct_trait_categories", "distinct_protein_coding_genes",
             "n_tissues", "n_coloc_groups", "n_rare_groups"),
    names_to = "metric",
    values_to = "value"
  ) |>
  dplyr::mutate(
    metric = factor(.data$metric,
      levels = c("distinct_trait_categories", "distinct_protein_coding_genes",
                 "n_tissues", "n_coloc_groups", "n_rare_groups"),
      labels = c("Trait categories", "Protein-coding genes",
                 "Distinct tissues", "Colocalization groups", "Rare result groups")
    ),
    gene = factor(.data$gene, levels = rev(ad_summary$gene))
  )

ad_summary_plot <- ggplot2::ggplot(ad_summary_long, ggplot2::aes(x = value, y = gene, fill = metric)) +
  ggplot2::geom_col(show.legend = FALSE, width = 0.5) +
  ggplot2::facet_wrap(~ metric, scales = "free_x", nrow = 1) +
  ggplot2::labs(x = NULL, y = NULL, title = "Alzheimer's disease gene targets: pleiotropy and evidence summary") +
  ggplot2::theme_bw(base_size = 9) +
  ggplot2::theme(
    strip.text = ggplot2::element_text(face = "bold"),
    panel.grid.minor = ggplot2::element_blank(),
    panel.grid.major.y = ggplot2::element_blank(),
    axis.text.y = ggplot2::element_text(size = 5)
  )

ad_summary_plot

trem2 <- gene("TREM2", include_associations = TRUE)
trem2$gene$distinct_trait_categories
#> [1] 1
trem2$gene$distinct_protein_coding_genes
#> [1] 2

trem2_ad_coloc_group_ids <- unique(
  trem2$coloc_groups[trem2$coloc_groups$trait_id %in% ad_trait_ids, ]$coloc_group_id
)
trem2_ad_rare_group_ids <- if (!is.null(trem2$rare_results) && nrow(trem2$rare_results) > 0) {
  unique(trem2$rare_results[trem2$rare_results$trait_id %in% ad_trait_ids, ]$rare_result_group_id)
} else {
  character(0)
}

trem2_ad_groups <- trem2$coloc_groups[
  trem2$coloc_groups$coloc_group_id %in% trem2_ad_coloc_group_ids,
]
trem2_ad_rare <- if (length(trem2_ad_rare_group_ids) == 0L) {
  trem2$rare_results[0, ]
} else {
  trem2$rare_results[
    trem2$rare_results$rare_result_group_id %in% trem2_ad_rare_group_ids,
  ]
}

trem2_all_rows <- dplyr::bind_rows(
  trem2_ad_groups,
  if (!is.null(trem2_ad_rare) && nrow(trem2_ad_rare) > 0L) trem2_ad_rare else NULL
)

variant_maf <- trem2_all_rows |>
  dplyr::filter(!is.na(.data$eaf)) |>
  dplyr::mutate(maf = pmin(.data$eaf, 1 - .data$eaf)) |>
  dplyr::group_by(.data$variant_id) |>
  dplyr::summarise(maf = min(.data$maf, na.rm = TRUE), .groups = "drop")

make_panel_label <- function(display_snp, variant_id, maf_lookup) {
  snp <- dplyr::coalesce(as.character(display_snp), paste0("variant ", as.character(variant_id)))
  maf <- maf_lookup$maf[match(variant_id, maf_lookup$variant_id)]
  maf_str <- ifelse(is.na(maf), "", paste0(" (MAF=", formatC(maf, format = "f", digits = 3), ")"))
  paste0(snp, maf_str)
}

trem2_forest_coloc <- trem2_ad_groups |>
  dplyr::filter(!is.na(.data$beta), !is.na(.data$se), .data$se > 0) |>
  dplyr::mutate(
    lo = .data$beta - 1.96 * .data$se,
    hi = .data$beta + 1.96 * .data$se,
    row_label = forest_row_label(.data$trait_name, .data$data_type, .data$tissue),
    panel_id = make_panel_label(.data$display_snp, .data$variant_id, variant_maf),
    data_type_label = label_data_type(.data$data_type)
  ) |>
  dplyr::group_by(.data$coloc_group_id) |>
  dplyr::arrange(.data$beta, .by_group = TRUE) |>
  dplyr::mutate(row_label = factor(.data$row_label, levels = unique(.data$row_label))) |>
  dplyr::ungroup()

trem2_forest_rare <- if (is.null(trem2_ad_rare) || nrow(trem2_ad_rare) == 0L) {
  NULL
} else {
  rr <- trem2_ad_rare
  if (!"tissue" %in% names(rr)) rr$tissue <- NA_character_
  if (!"eaf" %in% names(rr)) rr$eaf <- NA_real_
  rr |>
    dplyr::filter(!is.na(.data$beta), !is.na(.data$se), .data$se > 0) |>
    dplyr::mutate(
      lo = .data$beta - 1.96 * .data$se,
      hi = .data$beta + 1.96 * .data$se,
      row_label = forest_row_label(.data$trait_name, .data$data_type, .data$tissue),
      panel_id = make_panel_label(.data$display_snp, .data$variant_id, variant_maf),
      data_type_label = label_data_type(.data$data_type)
    ) |>
    dplyr::group_by(.data$panel_id) |>
    dplyr::arrange(.data$beta, .by_group = TRUE) |>
    dplyr::mutate(row_label = factor(.data$row_label, levels = unique(.data$row_label))) |>
    dplyr::ungroup()
}

trem2_forest <- dplyr::bind_rows(trem2_forest_coloc, trem2_forest_rare) |>
  dplyr::mutate(panel_id = factor(.data$panel_id, levels = unique(.data$panel_id)))

cbf_palette <- c(
  "Gene Expression" = "#0072B2",
  "Splice Variant"  = "#E69F00",
  "Protein"         = "#009E73",
  "Plasma Protein"  = "#56B4E9",
  "Transcript"      = "#CC79A7",
  "Methylation"     = "#F0E442",
  "Metabolome"      = "#D55E00",
  "Cell Trait"      = "#999999",
  "Phenotype"       = "#000000"
)

ggplot2::ggplot(trem2_forest, ggplot2::aes(x = beta, y = row_label, colour = data_type_label)) +
  ggplot2::geom_vline(xintercept = 0, linetype = 2, colour = "grey55") +
  ggplot2::geom_errorbarh(
    ggplot2::aes(xmin = lo, xmax = hi), height = 0, linewidth = 0.35
  ) +
  ggplot2::geom_point(size = 2) +
  ggplot2::facet_wrap(~ panel_id, scales = "free_y", ncol = 1) +
  ggplot2::scale_colour_manual(values = cbf_palette) +
  ggplot2::labs(
    x = "Effect (beta) with 95% CI", y = NULL,
    colour = "Data type", title = "TREM2 \u2014 Alzheimer's disease"
  ) +
  ggplot2::theme_bw() +
  ggplot2::theme(
    strip.text = ggplot2::element_text(face = "bold"),
    panel.grid.minor = ggplot2::element_blank(),
    legend.position = "bottom",
    legend.direction = "vertical"
  ) +
  ggplot2::guides(colour = ggplot2::guide_legend(ncol = 1))
#> Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
#> ℹ Please use the `orientation` argument of `geom_errorbar()` instead.
#> This warning is displayed once per session.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.
#> `height` was translated to `width`.