kinome_template/pathway_analysis.R at main · CogDisResLab/kinome_template · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
# Pathway Analysis

suppressPackageStartupMessages({
  library(tidyverse)
  library(enrichR)
  library(writexl)
})

load_reference_data <- function() {
  load_specific_reference <- function(chip) {
    id_file <- file.path(
      "reference_data",
      str_glue("{str_to_lower(chip)}_id_map.Rds")
    )
    hgnc_file <- file.path(
      "reference_data",
      str_glue("{str_to_lower(chip)}_hgnc_map.Rds")
    )

    id <- read_rds(id_file)
    hgnc <- read_rds(hgnc_file)

    peptide_map <- id |>
      inner_join(hgnc)
  }

  stk <- load_specific_reference("STK")
  ptk <- load_specific_reference("PTK")

  stk |>
    bind_rows(ptk)
}

perform_pathway_enrichment <- function(filename, peptide_map,
                                       threshold = 0.10, databases = NULL) {
  dpp_data <- read_csv(file.path("results", filename)) |>
    group_by(Peptide) |>
    filter(abs(LFC) == max(abs(LFC))) |>
    ungroup() |>
    select(Peptide, LFC) |>
    inner_join(peptide_map)

  lower_bound <- quantile(dpp_data[["LFC"]], threshold)
  upper_bound <- quantile(dpp_data[["LFC"]], 1 - threshold)

  selected_genes <- dpp_data |>
    filter(LFC >= upper_bound | LFC <= lower_bound) |>
    pull(Gene)

  if (is.null(databases)) {
    databases <- c(
      "GO_Molecular_Function_2025",
      "GO_Cellular_Component_2025",
      "GO_Biological_Process_2025",
      "KEGG_2021_Human",
      "Reactome_Pathways_2024"
    )
  }

  enriched <- enrichr(selected_genes, databases)
}

peptide_map <- load_reference_data()


dpp_files <- list.files("results", "dpp") |>
  set_names(~ .x |> str_extract("-dpp_(.*)\\.csv", 1L))

dpp_names <- names(dpp_files)

run_prefix <- dpp_files |>
  basename() |>
  str_extract("(.*)-dpp.*", 1L)

outfiles <- file.path(
  "results",
  str_glue_data(
    list(
      run_prefix = run_prefix,
      dpp_names = dpp_names
    ),
    "{run_prefix}-{dpp_names}-pathways.xlsx"
  )
)

results <- dpp_files |>
  map(
    ~ perform_pathway_enrichment(.x, peptide_map)
  ) |>
  map2(
    outfiles,
    ~ write_xlsx(
      .x,
      .y
    )
  )