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