Using fgsea with edgeR results
2025-01-24
Last updated: 2025-01-24
Checks: 7 0
Knit directory: muse/
This reproducible R Markdown analysis was created with workflowr (version 1.7.1). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20200712) was run prior to running
the code in the R Markdown file. Setting a seed ensures that any results
that rely on randomness, e.g. subsampling or permutations, are
reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version eef0287. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for
the analysis have been committed to Git prior to generating the results
(you can use wflow_publish or
wflow_git_commit). workflowr only checks the R Markdown
file, but you know if there are other scripts or data files that it
depends on. Below is the status of the Git repository when the results
were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: data/pbmc3k.csv
Ignored: data/pbmc3k.csv.gz
Ignored: data/pbmc3k/
Ignored: r_packages_4.4.0/
Ignored: r_packages_4.4.1/
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were
made to the R Markdown (analysis/fgsea_edger.Rmd) and HTML
(docs/fgsea_edger.html) files. If you’ve configured a
remote Git repository (see ?wflow_git_remote), click on the
hyperlinks in the table below to view the files as they were in that
past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | eef0287 | Dave Tang | 2025-01-24 | Using fgsea with edgeR results |
Installation
First install fgsea.
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if (!require("fgsea", quietly = TRUE))
BiocManager::install("fgsea")
library(fgsea)edgeR results
An example differential gene expression results table.
edger_res <- readr::read_csv("https://raw.githubusercontent.com/davetang/muse/refs/heads/main/data/13970886_edger_res.csv", show_col_types = FALSE)
head(edger_res)# A tibble: 6 × 6
ensembl_gene_id logFC logCPM F PValue adjusted_pvalue
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 ENSG00000000003 2.73 4.83 4.28 0.0684 0.109
2 ENSG00000000005 -7.00 0.541 17.6 0.00216 0.0138
3 ENSG00000000419 0.120 5.34 0.114 0.743 0.776
4 ENSG00000000457 -0.708 5.31 3.35 0.0993 0.145
5 ENSG00000000460 -0.897 3.95 2.66 0.136 0.186
6 ENSG00000000938 1.54 5.60 1.86 0.205 0.258 Add ranking metric.
edger_res |>
dplyr::mutate(rank_metric = logFC * -log10(PValue)) -> edger_resConvert to Entrez Gene IDs
Use {org.Hs.eg.db}.
if (!require("org.Hs.eg.db", quietly = TRUE))
BiocManager::install("org.Hs.eg.db")
Attaching package: 'BiocGenerics'The following objects are masked from 'package:stats':
IQR, mad, sd, var, xtabsThe 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, saveRDS, setdiff,
table, tapply, union, unique, unsplit, which.max, which.minWelcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
Attaching package: 'S4Vectors'The following object is masked from 'package:utils':
findMatchesThe following objects are masked from 'package:base':
expand.grid, I, unnamesuppressPackageStartupMessages(library(org.Hs.eg.db))Convert to Entrez Gene IDs.
AnnotationDbi::select(
org.Hs.eg.db,
keys = edger_res$ensembl_gene_id,
columns=c("ENSEMBL","ENTREZID"),
keytype="ENSEMBL"
) -> ensembl_to_entrez'select()' returned 1:many mapping between keys and columnsensembl_to_entrez <- dplyr::rename(ensembl_to_entrez, "ensembl_gene_id" = ENSEMBL)
head(ensembl_to_entrez) ensembl_gene_id ENTREZID
1 ENSG00000000003 7105
2 ENSG00000000005 64102
3 ENSG00000000419 8813
4 ENSG00000000457 57147
5 ENSG00000000460 55732
6 ENSG00000000938 2268Number of NAs.
table(is.na(ensembl_to_entrez$ENTREZID))
FALSE TRUE
28722 10968 Hallmark gene sets
Use {msigdb}.
if (!require("msigdb", quietly = TRUE))
BiocManager::install("msigdb")
if (!require("ExperimentHub", quietly = TRUE))
BiocManager::install("ExperimentHub")
Attaching package: 'AnnotationHub'The following object is masked from 'package:Biobase':
cacheif (!require("GSEABase", quietly = TRUE))
BiocManager::install("GSEABase")
Attaching package: 'graph'The following object is masked from 'package:XML':
addNodesuppressPackageStartupMessages(library(msigdb))
suppressPackageStartupMessages(library(ExperimentHub))
suppressPackageStartupMessages(library(GSEABase))Query an ExperimentHub object.
eh <- ExperimentHub(ask = FALSE)
AnnotationHub::query(x = eh, pattern = 'msigdb')ExperimentHub with 49 records
# snapshotDate(): 2024-10-24
# $dataprovider: Broad Institute, Emory University, EBI
# $species: Homo sapiens, Mus musculus
# $rdataclass: GSEABase::GeneSetCollection, list, data.frame
# additional mcols(): taxonomyid, genome, description,
# coordinate_1_based, maintainer, rdatadateadded, preparerclass, tags,
# rdatapath, sourceurl, sourcetype
# retrieve records with, e.g., 'object[["EH5421"]]'
title
EH5421 | msigdb.v7.2.hs.SYM
EH5422 | msigdb.v7.2.hs.EZID
EH5423 | msigdb.v7.2.mm.SYM
EH5424 | msigdb.v7.2.mm.EZID
EH6727 | MSigDB C8 MANNO MIDBRAIN
... ...
EH8296 | msigdb.v7.5.1.hs.SYM
EH8297 | msigdb.v7.5.1.mm.EZID
EH8298 | msigdb.v7.5.1.mm.idf
EH8299 | msigdb.v7.5.1.mm.SYM
EH8300 | imex_hsmm_0722 Latest version.
AnnotationHub::query(x = eh, pattern = 'msigdb.*hs.EZID') |>
tail(1) -> msigdb_hs_latest
msigdb_hs_latestExperimentHub with 1 record
# snapshotDate(): 2024-10-24
# names(): EH8294
# package(): msigdb
# $dataprovider: Broad Institute
# $species: Homo sapiens
# $rdataclass: GSEABase::GeneSetCollection
# $rdatadateadded: 2023-07-03
# $title: msigdb.v7.5.1.hs.EZID
# $description: Gene expression signatures (Homo sapiens) from the Molecular...
# $taxonomyid: 9606
# $genome: NA
# $sourcetype: XML
# $sourceurl: https://data.broadinstitute.org/gsea-msigdb/msigdb/release/7.5...
# $sourcesize: NA
# $tags: c("Homo_sapiens_Data", "Mus_musculus_Data")
# retrieve record with 'object[["EH8294"]]' Download.
msigdb_hs_ezid <- eh[[names(msigdb_hs_latest)]]see ?msigdb and browseVignettes('msigdb') for documentationloading from cacheCollections.
table(sapply(lapply(msigdb_hs_ezid, collectionType), bcCategory))
c1 c2 c3 c4 c5 c6 c7 c8 h
299 6180 3726 858 28005 189 5219 700 50 Create gene lists from the Hallmark collection.
wanted <- sapply(lapply(msigdb_hs_ezid, collectionType), bcCategory) == "h"
hallmark_gs <- msigdb_hs_ezid[wanted]
hallmark_gs_list <- lapply(hallmark_gs, geneIds)
class(hallmark_gs_list)[1] "list"names(hallmark_gs_list) <- names(hallmark_gs)
head(hallmark_gs_list)$HALLMARK_TNFA_SIGNALING_VIA_NFKB
[1] "3726" "2920" "467" "4792" "7128" "5743" "2919" "8870"
[9] "9308" "6364" "2921" "23764" "4791" "7127" "1839" "1316"
[17] "330" "5329" "7538" "3383" "3725" "1960" "3553" "597"
[25] "23645" "80149" "6648" "4929" "3552" "5971" "7185" "7832"
[33] "1843" "1326" "2114" "2152" "6385" "1958" "3569" "7124"
[41] "23135" "4790" "3976" "5806" "8061" "3164" "182" "6351"
[49] "2643" "6347" "1827" "1844" "10938" "9592" "5966" "8837"
[57] "8767" "4794" "8013" "22822" "51278" "8744" "2669" "1647"
[65] "3627" "10769" "8553" "1959" "9021" "11182" "5734" "1847"
[73] "5055" "4783" "5054" "10221" "25976" "5970" "329" "6372"
[81] "9516" "7130" "960" "3624" "5328" "4609" "3604" "6446"
[89] "10318" "10135" "2355" "10957" "3398" "969" "3575" "1942"
[97] "7262" "5209" "6352" "79693" "3460" "8878" "10950" "4616"
[105] "8942" "50486" "694" "4170" "7422" "5606" "1026" "3491"
[113] "10010" "3433" "3606" "7280" "3659" "2353" "4973" "388"
[121] "374" "4814" "65986" "8613" "9314" "6373" "6303" "1435"
[129] "1880" "56937" "5791" "7097" "57007" "7071" "4082" "3914"
[137] "1051" "9322" "2150" "687" "3949" "7050" "127544" "55332"
[145] "2683" "11080" "1437" "5142" "8303" "5341" "6776" "23258"
[153] "595" "23586" "8877" "941" "25816" "57018" "2526" "9034"
[161] "80176" "8848" "9334" "150094" "23529" "4780" "2354" "5187"
[169] "10725" "490" "3593" "3572" "9120" "19" "3280" "604"
[177] "8660" "6515" "1052" "51561" "4088" "6890" "9242" "64135"
[185] "3601" "79155" "602" "24145" "24147" "1906" "10209" "650"
[193] "1846" "10611" "23308" "9945" "10365" "3371" "5271" "4084"
$HALLMARK_HYPOXIA
[1] "5230" "5163" "2632" "5211" "226" "2026" "5236" "10397"
[9] "3099" "230" "2821" "4601" "6513" "5033" "133" "8974"
[17] "2023" "5214" "205" "26355" "5209" "7422" "665" "7167"
[25] "30001" "55818" "901" "3939" "2997" "2597" "8553" "51129"
[33] "3725" "5054" "4015" "2645" "8497" "23764" "54541" "6515"
[41] "3486" "4783" "2353" "3516" "3098" "10370" "3669" "2584"
[49] "26118" "5837" "6781" "23036" "694" "123" "1466" "7436"
[57] "23210" "2131" "2152" "5165" "55139" "7360" "229" "8614"
[65] "54206" "2027" "10957" "3162" "5228" "26330" "9435" "55076"
[73] "63827" "467" "857" "272" "2719" "3340" "8660" "8819"
[81] "2548" "6385" "8987" "8870" "5313" "3484" "5329" "112464"
[89] "8839" "9215" "25819" "6275" "58528" "7538" "1956" "1907"
[97] "3423" "1026" "6095" "1843" "4282" "5507" "10570" "11015"
[105] "1837" "136" "9957" "284119" "2908" "1316" "2239" "3491"
[113] "7128" "771" "3073" "633" "23645" "55276" "5292" "25824"
[121] "55577" "1027" "680" "8277" "4493" "538" "4502" "9672"
[129] "25976" "5317" "302" "5224" "1649" "5578" "2542" "7852"
[137] "1944" "1356" "8609" "1490" "9469" "7163" "56925" "124872"
[145] "10891" "596" "2651" "3036" "54800" "949" "6576" "6383"
[153] "839" "7428" "2309" "5155" "126792" "6518" "8406" "1942"
[161] "2745" "57007" "5066" "7045" "1634" "6478" "51316" "2203"
[169] "8459" "5260" "4627" "1028" "9380" "5105" "3623" "3309"
[177] "8509" "23327" "7162" "7511" "3569" "6533" "4214" "3948"
[185] "9590" "26136" "3798" "3906" "1289" "2817" "3069" "10994"
[193] "1463" "7052" "2113" "3219" "8991" "2355" "6820" "7043"
$HALLMARK_CHOLESTEROL_HOMEOSTASIS
[1] "2224" "1595" "3422" "2222" "1717" "6713" "3157" "50814"
[9] "4047" "4597" "3949" "7108" "230" "10682" "6319" "10654"
[17] "4598" "4023" "6309" "9415" "3156" "51478" "312" "6721"
[25] "5833" "55902" "467" "127" "23474" "1891" "875" "2990"
[33] "2194" "3958" "22809" "308" "94241" "1119" "2946" "39"
[41] "552" "5359" "1191" "54206" "57761" "58191" "51330" "71"
[49] "182" "5641" "26270" "493869" "10957" "118429" "114569" "928"
[57] "5468" "2731" "6811" "134429" "1499" "27346" "116496" "5165"
[65] "5329" "7869" "2770" "20" "6311" "4783" "214" "2171"
[73] "6282" "132864"
$HALLMARK_MITOTIC_SPINDLE
[1] "9181" "23332" "3832" "9493" "57679" "382" "4650" "4627"
[9] "10426" "9793" "29127" "57580" "50650" "4926" "6711" "11004"
[17] "3799" "7272" "324" "11190" "5048" "10435" "9371" "55704"
[25] "56992" "332" "116840" "4763" "7248" "996" "11064" "114791"
[33] "24137" "22919" "55785" "675" "5347" "5921" "4751" "8936"
[41] "7153" "7204" "9826" "10300" "9055" "54443" "55755" "9126"
[49] "10844" "9700" "55201" "201176" "9732" "29901" "3619" "394"
[57] "2934" "10276" "10128" "23637" "2317" "64411" "121512" "29"
[65] "55835" "4690" "1063" "9585" "10163" "4628" "1062" "9266"
[73] "4281" "3831" "57787" "127829" "9702" "8409" "393" "23580"
[81] "163786" "9113" "4983" "8976" "4296" "6654" "25" "7074"
[89] "23095" "6453" "134549" "8440" "9787" "613" "10048" "2037"
[97] "10801" "11104" "51174" "22974" "3797" "357" "85378" "6709"
[105] "23022" "23647" "9735" "84376" "25777" "58526" "1739" "2316"
[113] "79658" "8476" "23365" "4082" "51199" "5108" "10928" "7430"
[121] "85464" "983" "22930" "10160" "11346" "54509" "1894" "2035"
[129] "51735" "3835" "84333" "6780" "396" "6790" "26271" "51203"
[137] "5829" "9564" "23607" "11214" "10013" "22994" "3996" "23192"
[145] "5116" "7840" "11133" "667" "22920" "151987" "9411" "9462"
[153] "9133" "80119" "5922" "4739" "8243" "81" "5311" "7461"
[161] "998" "10403" "9874" "9344" "6904" "832" "1794" "2017"
[169] "10051" "10565" "7277" "4001" "10006" "6093" "55125" "699"
[177] "50628" "64857" "253260" "10018" "1778" "6624" "8874" "140735"
[185] "4643" "274" "4853" "5981" "10611" "89941" "8470" "11135"
[193] "7414" "6249" "23012" "7531" "9771" "55722" "1453"
$HALLMARK_WNT_BETA_CATENIN_SIGNALING
[1] "4609" "1499" "3714" "4851" "28514" "8313" "5664" "8321" "4855"
[10] "51176" "8312" "85407" "81029" "8454" "182" "9794" "2648" "2770"
[19] "7475" "5727" "9612" "27121" "3066" "22943" "6932" "7471" "8650"
[28] "6868" "1856" "5467" "23385" "10014" "894" "10023" "1454" "3516"
[37] "8325" "7157" "6502" "23493" "23462" "79885"
$HALLMARK_TGF_BETA_SIGNALING
[1] "7046" "4092" "7040" "64750" "57154" "659" "6498" "6497" "90"
[10] "56937" "9612" "5054" "3726" "4086" "4091" "23645" "7050" "5045"
[19] "4088" "2280" "6885" "657" "1499" "28996" "7071" "650" "2022"
[28] "324" "5494" "331" "999" "3397" "7044" "1028" "51592" "11031"
[37] "7082" "6574" "1025" "3399" "9241" "51742" "3460" "3398" "5499"
[46] "6711" "25937" "8412" "7057" "2339" "3065" "7323" "4053" "387" fgsea analysis
Rank metric distribution.
plot(hist(edger_res$rank_metric, breaks = 50))
Create ranks vector.
edger_res |>
dplyr::inner_join(y = ensembl_to_entrez, by = "ensembl_gene_id", relationship = "many-to-many") |>
dplyr::filter(!is.na(ENTREZID)) |>
dplyr::group_by(ENTREZID) |>
dplyr::mutate(ambiguous = ifelse(dplyr::n()>1, TRUE, FALSE)) |>
dplyr::filter(!ambiguous) -> res
my_ranks <- res$rank_metric
my_names <- as.character(res$ENTREZID)
names(my_ranks) <- my_names
plot(hist(my_ranks, breaks = 50))
The fgsea() function runs the pre-ranked gene set
enrichment analysis.
set.seed(1984)
fgseaRes <- fgsea(
pathways = hallmark_gs_list,
stats = my_ranks,
minSize=15,
maxSize=500,
nPermSimple = 200000
)Warning in preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam, : There are ties in the preranked stats (2.44% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.Warning in fgseaMultilevel(pathways = pathways, stats = stats, minSize =
minSize, : There were 36 pathways for which P-values were not calculated
properly due to unbalanced (positive and negative) gene-level statistic values.
For such pathways pval, padj, NES, log2err are set to NA. You can try to
increase the value of the argument nPermSimple (for example set it nPermSimple
= 2000000)Top 6 enriched pathways.
head(fgseaRes[order(pval), ]) pathway pval padj log2err
<char> <num> <num> <num>
1: HALLMARK_NOTCH_SIGNALING 3.559469e-05 0.0004983257 5.573322e-01
2: HALLMARK_ANGIOGENESIS 1.063389e-03 0.0074437263 4.550599e-01
3: HALLMARK_SPERMATOGENESIS 3.347883e-01 1.0000000000 4.547319e-03
4: HALLMARK_PANCREAS_BETA_CELLS 4.714677e-01 1.0000000000 3.415898e-03
5: HALLMARK_KRAS_SIGNALING_DN 8.693607e-01 1.0000000000 1.250537e-03
6: HALLMARK_APICAL_SURFACE 9.995800e-01 1.0000000000 6.844827e-05
ES NES size leadingEdge
<num> <num> <int> <list>
1: 0.6388199 2.4231911 31 2648, 11....
2: 0.4445579 1.7490871 35 6696, 70....
3: -0.5487978 -1.0257935 135 56903, 8....
4: -0.5614137 -1.0146014 40 5080, 47....
5: -0.5031756 -0.9463782 194 56154, 2....
6: -0.2987335 -0.5417546 44 116085, ....Plot the most significantly enriched pathway.
plotEnrichment(
hallmark_gs_list[[head(fgseaRes[order(pval), ], 1)$pathway]],
my_ranks
) +
ggplot2::labs(title=head(fgseaRes[order(pval), ], 1)$pathway)
sessionInfo()R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.5 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.20.so; LAPACK version 3.10.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] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] GSEABase_1.68.0 graph_1.84.1 annotate_1.84.0
[4] XML_3.99-0.17 ExperimentHub_2.14.0 AnnotationHub_3.14.0
[7] BiocFileCache_2.14.0 dbplyr_2.5.0 msigdb_1.14.0
[10] org.Hs.eg.db_3.20.0 AnnotationDbi_1.68.0 IRanges_2.40.1
[13] S4Vectors_0.44.0 Biobase_2.66.0 BiocGenerics_0.52.0
[16] fgsea_1.32.2 BiocManager_1.30.25 workflowr_1.7.1
loaded via a namespace (and not attached):
[1] DBI_1.2.3 rlang_1.1.4 magrittr_2.0.3
[4] git2r_0.35.0 compiler_4.4.1 RSQLite_2.3.9
[7] getPass_0.2-4 png_0.1-8 callr_3.7.6
[10] vctrs_0.6.5 stringr_1.5.1 pkgconfig_2.0.3
[13] crayon_1.5.3 fastmap_1.2.0 XVector_0.46.0
[16] labeling_0.4.3 utf8_1.2.4 promises_1.3.0
[19] rmarkdown_2.28 tzdb_0.4.0 UCSC.utils_1.2.0
[22] ps_1.8.1 purrr_1.0.2 bit_4.5.0
[25] xfun_0.48 zlibbioc_1.52.0 cachem_1.1.0
[28] GenomeInfoDb_1.42.1 jsonlite_1.8.9 blob_1.2.4
[31] highr_0.11 later_1.3.2 BiocParallel_1.40.0
[34] parallel_4.4.1 R6_2.5.1 bslib_0.8.0
[37] stringi_1.8.4 jquerylib_0.1.4 Rcpp_1.0.13
[40] knitr_1.48 readr_2.1.5 httpuv_1.6.15
[43] Matrix_1.7-0 tidyselect_1.2.1 rstudioapi_0.17.1
[46] yaml_2.3.10 codetools_0.2-20 curl_6.0.1
[49] processx_3.8.4 lattice_0.22-6 tibble_3.2.1
[52] withr_3.0.2 KEGGREST_1.46.0 evaluate_1.0.1
[55] Biostrings_2.74.1 pillar_1.9.0 filelock_1.0.3
[58] whisker_0.4.1 generics_0.1.3 vroom_1.6.5
[61] rprojroot_2.0.4 BiocVersion_3.20.0 hms_1.1.3
[64] ggplot2_3.5.1 munsell_0.5.1 scales_1.3.0
[67] xtable_1.8-4 glue_1.8.0 tools_4.4.1
[70] data.table_1.16.2 fs_1.6.4 fastmatch_1.1-4
[73] cowplot_1.1.3 grid_4.4.1 colorspace_2.1-1
[76] GenomeInfoDbData_1.2.13 cli_3.6.3 rappdirs_0.3.3
[79] fansi_1.0.6 dplyr_1.1.4 gtable_0.3.6
[82] sass_0.4.9 digest_0.6.37 farver_2.1.2
[85] memoise_2.0.1 htmltools_0.5.8.1 lifecycle_1.0.4
[88] httr_1.4.7 mime_0.12 bit64_4.5.2