## ----knitr, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( comment = "#>", error = FALSE, warning = FALSE, eval = TRUE, message = FALSE ) ## ----install, eval = FALSE---------------------------------------------------- # if (!requireNamespace("BiocManager", quietly = TRUE)) { # install.packages("BiocManager")} # # BiocManager::install("EMMA") ## ----setup-------------------------------------------------------------------- library("EMMA") ## ----loadlibraries------------------------------------------------------------ library("EMMA") library("macrophage") library("DESeq2") library("org.Hs.eg.db") library("clusterProfiler") library("mosdef") library("topGO") library("GO.db") ## ----differential_expression-------------------------------------------------- # load data data(gse, package = "macrophage") # set up design dds_macrophage <- DESeqDataSet(gse, design = ~ line + condition) # preprocess rownames(dds_macrophage) <- substr(rownames(dds_macrophage), 1, 15) keep <- rowSums(counts(dds_macrophage) >= 10) >= 6 dds_macrophage <- dds_macrophage[keep, ] # set seed for reproducibility set.seed(2711) # sample randomly for 2k genes selected_genes <- sample(rownames(dds_macrophage), 2000) dds_macrophage <- dds_macrophage[selected_genes, ] # run DESeq dds_macrophage <- DESeq(dds_macrophage) # get de res for 1st contrast IFNg_vs_naive <- results(dds_macrophage, contrast = c("condition", "IFNg", "naive"), lfcThreshold = 1, alpha = 0.05) IFNg_vs_naive <- lfcShrink(dds_macrophage, coef = "condition_IFNg_vs_naive", res = IFNg_vs_naive, type = "apeglm") IFNg_vs_naive$SYMBOL <- rowData(dds_macrophage)$SYMBOL # get de res for 2st contrast SL1344_vs_naive <- results(dds_macrophage, contrast = c("condition", "SL1344", "naive"), lfcThreshold = 1, alpha = 0.05) SL1344_vs_naive <- lfcShrink(dds_macrophage, coef = "condition_SL1344_vs_naive", res = SL1344_vs_naive, type = "apeglm") SL1344_vs_naive$SYMBOL <- rowData(dds_macrophage)$SYMBOL # sort by adjusted p value de_list <- list( IFNg_vs_naive = IFNg_vs_naive, SL1344_vs_naive = SL1344_vs_naive ) de_list <- lapply(de_list, function(df) { df <- df[order(df$padj), ] df <- df[!is.na(df$padj) & df$padj <= 0.05, ] df }) # set background gene list gene_universe <- rownames(dds_macrophage) ## ----withoutemma-------------------------------------------------------------- de_res <- de_list$IFNg_vs_naive fea_res <- enrichGO(gene = rownames(de_res), universe = gene_universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.1) ## ----EMMA_run_1--------------------------------------------------------------- # perform FEA, but with EMMA! fea_res <- enrichGO(gene = rownames(de_res), universe = gene_universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.1) |> EMMA_run() # simply pipe your call to EMMA_run() # check res fea_res ## ----EMMA_run_2--------------------------------------------------------------- # you can also pass the function name and its namespace # e.g. `clusterProfiler::enrichGO(...)` fea_res_nobg <- EMMA_run( clusterProfiler::enrichGO( # no universe set gene = rownames(de_res), keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.1, readable = TRUE ) ) ## ----EMMA_show---------------------------------------------------------------- EMMA_show(fea_res) ## ----EMMA_get_record---------------------------------------------------------- emma_record <- EMMA_get_record(fea_res) # get all the record emma_record ## ----EMMA_method-------------------------------------------------------------- # get the method record emma_record$method ## ----argument_form------------------------------------------------------------ fea_res_no_param <- enrichGO(gene = rownames(de_res), universe = gene_universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.1, readable = TRUE) |> EMMA_run(args_form = "unevaluated") # when we don't want the values stored # else set to evaluated (default) # check EMMA_get_record(fea_res_no_param) ## ----EMMA_explain, results = 'asis', message=TRUE----------------------------- EMMA_explain(fea_res, get_citation = TRUE) ## ----mosdef_eg---------------------------------------------------------------- mosdef_fea_res <- mosdef::run_goseq(de_genes = rownames(de_res), bg_genes = gene_universe, mapping = "org.Hs.eg.db", id = "ensGene", genome = "hg19") |> EMMA_run(store_session_info = FALSE, args_form = "unevaluated") # quick inspection EMMA_get_record(mosdef_fea_res) ## ----custom_eg---------------------------------------------------------------- # a custom function (not from a package) my_custom_function <- function(gene, universe = NULL, ontology = "BP", id_type = "ENTREZID", org_db_name = "org.Hs.eg.db", organism = "hsapiens") { res1 <- mosdef::run_topGO(de_genes = gene, bg_genes = gene_universe, ontology = ontology, gene_id = id_type, mapping = org_db_name, add_gene_to_terms = TRUE) res2 <- gprofiler2::gost(query = gene, organism = organism, custom_bg = gene_universe) return(list(topGO_res = res1, gost_res = res2 )) } # run analysis with EMMA frankenstein_fea <- my_custom_function( gene = rownames(de_res), universe = gene_universe, ontology = "BP", id_type = "ENSEMBL", org_db_name = "org.Hs.eg.db", organism = "hsapiens" ) |> EMMA_run(store_session_info = FALSE, args_form = "unevaluated") # quick inspection EMMA_get_record(frankenstein_fea) ## ----add_custom_metadata------------------------------------------------------ frankenstein_fea2 <- EMMA_add_custom_metadata(res = frankenstein_fea, extra = list( wrapped_function_topGO = "runTest", notes = "any other meaningful info")) EMMA_get_record(frankenstein_fea2)$extra ## ----emma_and_dde------------------------------------------------------------- dde <- DeeDeeExperiment::DeeDeeExperiment(sce = dds_macrophage, de_results = IFNg_vs_naive, enrich_results = list( IFNg_vs_naive = fea_res_no_param) ) fea <- DeeDeeExperiment::getFEA(dde, format = "original") EMMA_get_record(fea) ## ----EMMA_freeze-------------------------------------------------------------- if (requireNamespace("renv", quietly = TRUE)) { project_path <- tempdir() dir.create(project_path) EMMA_freeze(project = project_path, file = "analysis.lock", pkgs = loadedNamespaces(), prompt = FALSE, force = TRUE) } ## ----emma_with_lapply--------------------------------------------------------- fea_res_list <- lapply(de_list, function(dea){ #### perform fea with EMMA_run #### enrichGO(gene = rownames(dea), universe = gene_universe, keyType = "ENSEMBL", OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.1, readable = TRUE) |> EMMA_run() }) #### print the metadata summary for each contrast #### invisible(lapply(names(fea_res_list), function(nm) { cat("\n###", nm) EMMA_show(fea_res_list[[nm]]) })) #### add more info to the metadata #### fea_res_list <- setNames(lapply(names(fea_res_list), function(nm) { EMMA_add_custom_metadata(fea_res_list[[nm]], extra = list(contrast_name = nm)) }), names(fea_res_list)) #### get the raw metadata #### record_list <- list() record_list <- setNames(lapply(names(fea_res_list), function(mn){ EMMA_get_record(fea_res_list[[mn]]) }), names(fea_res_list)) record_list$SL1344_vs_naive$extra record_list$SL1344_vs_naive$method$package_name ## ----sessioinfo--------------------------------------------------------------- sessionInfo()