library(MotrpacRatTraining6moData)
library(MotrpacRatTraining6moWATData)
library(MotrpacRatTraining6moWAT)
#> Loading required package: Biobase
#> Loading required package: BiocGenerics
#> 
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#> 
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
#>     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#>     Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
#>     tapply, union, unique, unsplit, which.max, which.min
#> Welcome to Bioconductor
#> 
#>     Vignettes contain introductory material; view with
#>     'browseVignettes()'. To cite Bioconductor, see
#>     'citation("Biobase")', and for packages 'citation("pkgname")'.
#> 
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following object is masked from 'package:Biobase':
#> 
#>     combine
#> The following objects are masked from 'package:BiocGenerics':
#> 
#>     combine, intersect, setdiff, union
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(tidyr)
library(tibble)
library(purrr)
library(fgsea)
# Entrez to gene symbol conversion vector for leading edge
entrez_to_symbol <- pluck(PROT_DA, "MvF_SED") %>%
  filter(!is.na(entrez_gene)) %>%
  dplyr::select(entrez_gene, gene_symbol) %>%
  deframe()
head(entrez_to_symbol)
#>     26194     26195     26198     26196     26204     26199 
#>  "Mt-nd2"  "Mt-co1"  "Mt-co2" "Mt-atp8"  "Mt-co3"  "Mt-nd3"

# Human to rat gene conversion
human_to_rat <- RAT_TO_HUMAN_GENE %>%
  dplyr::select(HUMAN_ORTHOLOG_SYMBOL, RAT_NCBI_GENE_ID) %>%
  distinct() %>%
  deframe()
head(human_to_rat)
#>   A1BG   NAT2    ADA   CDH2   AKT3   MED6 
#> 140656 116632  24165  83501  29414 299180

# Convert human gene symbols to rat Entrez IDs
PROT_MITOCARTA <- MITOCARTA_HS %>%
  mutate(rat_entrez = map(human_genes,
                          ~ as.character(na.omit(human_to_rat[.x]))),
         set_size = lengths(rat_entrez),
         rat_entrez = map(rat_entrez, intersect, names(entrez_to_symbol)),
         set_size_post = lengths(rat_entrez),
         ratio = set_size_post / set_size) %>%
  filter(ratio >= 0.85,
         set_size <= 300, # Same as MSIGDB_PATHWAYS filter
         set_size_post >= 5)
# List of pathways to test
MITOCARTA_PATHWAYS <- PROT_MITOCARTA %>%
  dplyr::select(pathway, rat_entrez) %>%
  deframe()
length(MITOCARTA_PATHWAYS) # 68
#> [1] 68
head(MITOCARTA_PATHWAYS)
#> $`mt-tRNA modifications`
#>  [1] "362766" "252827" "291978" "290633" "63864"  "300852" "361191" "314548"
#>  [9] "304567" "64016"  "315550" "303067" "362586" "288914" "304012" "362754"
#> [17] "362976" "312616" "319113"
#> 
#> $`mt-rRNA modifications`
#>  [1] "295985" "305845" "502632" "304323" "360569" "307050" "499191" "298426"
#>  [9] "360796" "362416" "315550" "308140" "362586" "366012"
#> 
#> $Translation
#>   [1] "301254"    "498179"    "361184"    "361005"    "498000"    "681219"   
#>   [7] "295238"    "304919"    "367172"    "361641"    "362759"    "306879"   
#>  [13] "301463"    "288916"    "297113"    "361974"    "360821"    "114017"   
#>  [19] "294672"    "305256"    "312054"    "305317"    "307491"    "364070"   
#>  [25] "292028"    "363172"    "313867"    "297082"    "311748"    "305845"   
#>  [31] "315141"    "304323"    "360569"    "289491"    "691075"    "293666"   
#>  [37] "303746"    "299938"    "301250"    "297799"    "293754"    "171061"   
#>  [43] "292244"    "297372"    "301240"    "680747"    "309140"    "287302"   
#>  [49] "64360"     "295224"    "287635"    "497876"    "300974"    "301352"   
#>  [55] "291206"    "100363539" "290632"    "56281"     "303685"    "684304"   
#>  [61] "363023"    "287962"    "296551"    "299743"    "309440"    "301552"   
#>  [67] "287656"    "293054"    "294963"    "293149"    "309176"    "362517"   
#>  [73] "297601"    "361037"    "362388"    "299628"    "287356"    "691814"   
#>  [79] "303673"    "310653"    "363187"    "499185"    "292758"    "289143"   
#>  [85] "298517"    "688912"    "288621"    "301249"    "294230"    "289469"   
#>  [91] "362094"    "689432"    "683519"    "360594"    "498406"    "297459"   
#>  [97] "362216"    "361883"    "689025"    "294767"    "290850"    "296995"   
#> [103] "287126"    "297727"    "294696"    "296134"    "100360017" "113958"   
#> [109] "301371"    "311903"    "363289"    "100910751" "296462"    "305606"   
#> [115] "684274"    "683897"    "686234"    "361473"    "293128"    "499191"   
#> [121] "298426"    "691393"    "313429"    "690214"    "310856"    "500199"   
#> [127] "362681"    "309911"    "297969"    "307235"    "292268"    "292759"   
#> [133] "688717"    "360645"    "310672"    "308140"    "307210"    "679068"   
#> [139] "293481"    "309596"    "690654"    "287924"   
#> 
#> $`Mitochondrial ribosome`
#>  [1] "361005"    "295238"    "288916"    "289491"    "691075"    "293666"   
#>  [7] "303746"    "299938"    "301250"    "297799"    "293754"    "171061"   
#> [13] "292244"    "297372"    "301240"    "680747"    "309140"    "287302"   
#> [19] "64360"     "295224"    "287635"    "497876"    "300974"    "301352"   
#> [25] "291206"    "100363539" "290632"    "56281"     "303685"    "684304"   
#> [31] "363023"    "287962"    "296551"    "299743"    "309440"    "301552"   
#> [37] "287656"    "293054"    "294963"    "293149"    "309176"    "362517"   
#> [43] "297601"    "361037"    "362388"    "299628"    "287356"    "691814"   
#> [49] "303673"    "310653"    "363187"    "499185"    "292758"    "289143"   
#> [55] "298517"    "688912"    "288621"    "301249"    "294230"    "289469"   
#> [61] "362094"    "689432"    "683519"    "360594"    "498406"    "297459"   
#> [67] "362216"    "361883"    "689025"    "294767"    "290850"    "296995"   
#> [73] "287126"    "297727"    "294696"    "296134"    "100360017" "113958"   
#> [79] "301371"    "500199"   
#> 
#> $`Translation factors`
#>  [1] "498179" "498000" "681219" "114017" "294672" "303673" "311903" "305606"
#>  [9] "684274" "686234" "361473" "360645" "307210" "679068" "293481"
#> 
#> $`mt-tRNA synthetases`
#>  [1] "301254" "361184" "304919" "361641" "306879" "297113" "361974" "360821"
#>  [9] "307491" "364070" "292028" "363172" "293128" "313429" "309911" "297969"
#> [17] "292759" "310672" "309596" "690654" "287924"
## FGSEA
PROT_MITOCARTA_FGSEA <- map(PROT_DA, function(res_i) {
  stats <- rank_genes(res_i, genes = "entrez_gene")

  map(names(stats), function(contrast_i) {
    message(contrast_i)
    set.seed(0)
    fgseaMultilevel(pathways = MITOCARTA_PATHWAYS,
                    stats = stats[[contrast_i]],
                    nPermSimple = 10000, nproc = 1) %>%
      mutate(contrast = contrast_i)
  }) %>%
    data.table::rbindlist() %>%
    mutate(contrast = factor(contrast, levels = unique(contrast)),
           padj = p.adjust(pval, method = "BH"),
           leadingEdge_genes = map(leadingEdge,
                                   ~ na.omit(entrez_to_symbol[.x])),
           leadingEdge = map(leadingEdge, as.numeric)) %>% # consistency
    left_join(dplyr::select(PROT_MITOCARTA, pathway, hierarchy),
              by = "pathway") %>%
    relocate(hierarchy, .after = pathway) %>%
    relocate(contrast, .before = leadingEdge)
})
# Save
usethis::use_data(PROT_MITOCARTA_FGSEA, internal = FALSE,
                  overwrite = TRUE, version = 3, compress = "bzip2")