Package 'goat'

Description

cluster significant genesets from test_genesets() by geneset similarity (separately for each 'geneset source')

Usage

cluster_genesets(x, genelist, hclust_method = "ward.D2")

Arguments

Value

a list with elements genesets (param x), similarity, hc_row, hc_col

Description

naively darken a color by mixing in black

Usage

darken_color(color, frac = 0.1)

Arguments

Value

adjusted value for input color

Description

while the Bioconductor respository is extensive, contains data for many species and is a part of a larger infrastructure, it might contain outdated GO data when the user is not using the latest R version. If users are on an R version that is a few years old, so will the GO data from Bioconductor be.

As an alternative, we store gene2go data from NCBI (for Human genes only!) at the GOAT GitHub repository. This function allows for a convenient way to download this data and then parse the genesets.

Alternatively you can browse the file in the data branch of the GOAT GitHub repository and download these files manually, then load them via the GOAT R function load_genesets_go_fromfile().

To view all available data you can open this URL in a browser; https://github.com/ftwkoopmans/goat/tree/data

New data is automatically added biannually. The first available version is 2024-01-01, the next 2024-06-01, then 2025-01-01, and so on.

Usage

download_genesets_goatrepo(
  output_dir,
  type = "GO",
  version = "2024-01-01",
  ignore_cache = FALSE
)

Arguments

Value

result from respective geneset parser function. e.g. if parameter type was set to"GO" (default), this function returns the result of load_genesets_go_fromfile(). These data returned by this function is typically used as input for filter_genesets(), c.f. full example at documentation for test_genesets()

Examples

# note: this example will download 2 files of approx 10MB in total

# store the downloaded files in the following directory. Here, the temporary file
# directory is used. Alternatively, consider storing this data in a more permanent location.
# e.g. output_dir="~/data/go" on unix systems or output_dir="C:/data/go" on Windows
output_dir = tempdir()

# download data files with GO annotations, version 2024-01-01 (default parameter)
# these are then parsed with the load_genesets_go_fromfile() function
# if the files are already available at output_dir, the download step is skipped
genesets_asis = download_genesets_goatrepo(output_dir)

### for a basic example on how to use the data obtain here,
### refer to the example included at function documentation of: test_genesets()

Description

Downloads OMICs-based datasets that were used in the GOAT manuscript from the GOAT GitHub page. This file is cached in the output directory and only needs to be downloaded once. Multiple datasets are included and their names include the respective PubMed identifiers (PMID).

If you encounter technical difficulties, try to;

1. download the file by copy/pasting this URL into your browser: https://github.com/ftwkoopmans/goat/raw/main/analyses/goat_manuscript_datasets.rda

2. load the data in R using the following 2 lines of code, here assuming you stored the downloaded file at C:/data/goat_manuscript_datasets.rda

load("C:/data/goat_manuscript_datasets.rda")

genelist = goat_manuscript_datasets.rda[["Wingo 2020:mass-spec:PMID32424284"]]

Usage

download_goat_manuscript_data(output_dir, ignore_cache = FALSE)

Arguments

Value

a list of genelist data tables. The names of the list represent the datasets, values in the list are data tables that can be used as a "genelist" in the GOAT R package

Description

filter a geneset table; intersect with an array of genes-of-interest then apply cutoffs on min/max genes per geneset

Usage

filter_genesets(
  genesets,
  genelist,
  min_overlap = 10L,
  max_overlap = 1500L,
  max_overlap_fraction = 0.5,
  min_signif = NA,
  max_size = NA,
  dedupe = FALSE
)

Arguments

Value

the input genesets filtered for the subset of rows that match user's filter parameters

Description

https://stackoverflow.com/a/8197703

Usage

gg_color_hue(n)

Arguments

Value

a color code (string)

Description

note that this file lacks parent/child relations, so we only learn 'direct annotations'

Usage

go_gene2go(f, taxid_filter = 9606)

Arguments

Value

a tibble with columns; source, source_version, id, name, genes, ngenes

Description

note that we remove links between GO terms that are across GO domains (e.g. no CC to MF relations) The only supported relations are those that match this regex; "^(is_a:|relationship: part_of|relationship: regulates|relationship: positively_regulates|relationship: negatively_regulates)"

Usage

go_obo(f, rename_namespace = TRUE, remove_obsolete = TRUE)

Arguments

Value

tibble with ontology terms and their relations

Description

ASCII logo for this package

Usage

goat_logo()

Value

package logo as a string

Description

there parameters are used by goat to efficiently perform geneset testing without bootstrapping

Usage

goat_nulldistributions

Format

goat_nulldistributions

a data.frame with precomputed GOAT null distribution parameters

Description

Print package version and logo to console

Usage

goat_print_version()

Value

prints to console without returning a value

Description

simple wrapper around utils::packageVersion()

Usage

goat_version()

Value

package version as a string

Description

download link: https://www.genenames.org/download/statistics-and-files/ table: "Complete dataset download links" –>> "Complete HGNC approved dataset" –>> download the "TXT" table filename is typically something like hgnc_complete_set.txt URL; https://ftp.ebi.ac.uk/pub/databases/genenames/hgnc/tsv/hgnc_complete_set.txt

Usage

hgnc_idmap_table(filename)

Arguments

Value

a long-format table with columns; hgnc_id, hgnc_symbol, type, value

alternatively;

table: "Total Approved Symbols" –>> "TXT" / "text file in TSV format" filename is typically something like non_alt_loci_set.txt

Description

naively lighten a color by mixing in white

Usage

lighten_color(color, frac = 0.1)

Arguments

Value

adjusted value for input color

Description

parse genesets in GMT format where gene identifiers are numeric Entrez gene IDs

Usage

load_genesets_gmtfile(filename, label)

Arguments

Value

tibble with columns; source (character), source_version (character), id (character), name (character), genes (list), ngenes (int)

Example data;

URL: https://www.gsea-msigdb.org/gsea/msigdb/human/collections.jsp#C5 download this data: KEGG subset of curated pathways –>> NCBI (Entrez) Gene IDs filename should be something like "c2.cp.kegg.v2023.1.Hs.entrez.gmt"

Examples

# TODO: update the filename to your downloaded file
  f = "C:/DATA/c2.cp.kegg.v2023.1.Hs.entrez.gmt"
  if(file.exists(f)) genesets_asis = load_genesets_gmtfile(f, label = "KEGG")

Description

Download and import genesets from the GO database using the Bioconductor infrastructure. Use the goat::load_genesets_go_fromfile function for more fine-grained control over the GO database version that you use; it allows you to import NCBI gene2go files

Usage

load_genesets_go_bioconductor(include_child_annotations = TRUE)

Arguments

Details

Note that org.Hs.eg.db pulls data semi-annually from NCBI gene2go, but the GO database version returned by this function is tied to the version of the org.Hs.eg.db on your computer (this is controlled by the Bioconductor infrastructure).

The actual GO database version that is retrieved is returned by this function in the source_version column.

Value

table with columns; source (character), source_version (character), id (character), name (character), genes (list), ngenes (int)

Description

This function is used to load Gene Ontology (GO) genesets from files that you manually downloaded from the links below. This enables the use of the latest data from GO (in contrast, Bioconductor GO data may lag behind current data considerably). To construct genesets from available raw data, download the "gene2go" file (the gene annotations) from below NCBI link and download the GO OBO (ontology terms and relations to respective parent/child terms) from below geneontology.org link. Provide the full path to the downloaded file to this function. Both "gzipped" and "uncompressed" files are supported.

We encourage you to rename the files after your downloaded them such that the date of download in incorporated; this ensures you can always keep track of the GO database version that was used! For example, rename the downloaded "gene2go.gz" file to "gene2go_2024-01-31.gz".

Download link for gene2go file; https://ftp.ncbi.nih.gov/gene/DATA/gene2go.gz

Download link for gene ontology OBO file; http://current.geneontology.org/ontology/go.obo

Usage

load_genesets_go_fromfile(
  file_gene2go,
  file_goobo,
  include_child_annotations = TRUE
)

Arguments

Value

table with columns; source (character), source_version (character), id (character), name (character), genes (list), ngenes (int)

Examples

# TODO: update the filenames to your downloaded files
  file_gene2go = "C:/DATA/gene2go_2024-01-01.gz"
  file_goobo = "C:/DATA/go_2024-01-01.obo"
  if(file.exists(file_gene2go) && file.exists(file_goobo)) {
    genesets_asis = load_genesets_go_fromfile(file_gene2go, file_goobo)
  }

Description

Workflow;

• obtain the input file from; https://www.syngoportal.org

• click "bulk download SynGO release ..." for SynGO release of interest

• unzip

• call this function with the full file path to the 'syngo_ontologies.xlsx' file

Usage

load_genesets_syngo(filename, gene_database = "entrez")

Arguments

Value

table with columns; source (character), source_version (character), id (character), name (character), genes (list), ngenes (int)

Examples

# TODO: update the filename to your downloaded file
  f = "C:/DATA/SynGO_bulk_download_release_20210225/syngo_ontologies.xlsx"
  if(file.exists(f)) genesets_asis = load_genesets_syngo(f)

Description

-log10 transform a vector of p-values, replacing zeros with some limit/threshold

Usage

minlog10_fixzero(x, limit = 2.22e-16)

Arguments

Value

input parameter x transformed to -log10

Examples

pval = c(0, 10^-6, 0.001, 0.01, 1, NA, -Inf, Inf, NaN)
  cbind(
    input = pval,
    # default; replace zeros with typical R machine precision for doubles
    minlog10_default = minlog10_fixzero(pval),
    # alternatively, replace zero with lowest non-zero pvalue in input
    minlog10_limit_from_data = minlog10_fixzero(pval, limit = NA)
  )

Description

Adjust p-values for all genesets, grouped by 'source' then adjust for the number of 'sources'

Usage

padjust_genesets(
  genesets,
  method = "BH",
  cutoff = 0.01,
  correct_sources = TRUE
)

Arguments

Value

updated genesets table

Description

This can be convenient to prepare the significant/test/foreground set for classical ORA, e.g. test_genesets() with parameter method = "fisherexact". Note that the GOAT geneset enrichment algorithm does not use data in the 'signif' column of the input genelist.

Usage

partition_genes(
  genes,
  col,
  decreasing = FALSE,
  use_abs = FALSE,
  cutoff = NULL,
  fraction = NULL,
  topn = NULL
)

Arguments

Value

input table genes with results in the "signif" column

Examples

# note: this example will download 1 files of approx 4MB

# store the downloaded files in the following directory. Here, the temporary file
# directory is used. Alternatively, consider storing this data in a more permanent location.
# e.g. output_dir="~/data/goat" on unix systems or output_dir="C:/data/goat" on Windows
output_dir = tempdir()

# Download an example gene list, i.e. one of the datasets analyzed in the GOAT manuscript.
datasets = download_goat_manuscript_data(output_dir)
genelist = datasets$`Wingo 2020:mass-spec:PMID32424284`

# example 1: significant hits
genelist = partition_genes(genelist, col="pvalue_adjust", decreasing=FALSE, cutoff=0.01)
cat(sum(genelist$signif), "/", nrow(genelist), "are signif\n")

# example 2: abs(effectsize) >= 5
genelist = partition_genes(genelist, col="effectsize", decreasing=TRUE, use_abs=TRUE, cutoff=5)
cat(sum(genelist$signif), "/", nrow(genelist), "are signif\n")

# example 3: top 10% 'best' p-values
genelist = partition_genes(genelist, col="pvalue", decreasing=FALSE, fraction = 0.1)
cat(sum(genelist$signif), "/", nrow(genelist), "are signif\n")

Description

plot the geneset similarity matrix as a heatmap

Usage

plot_heatmap(
  x,
  output_dir,
  colors = grDevices::hcl.colors(100, "Viridis", rev = FALSE),
  fontsize = 10
)

Arguments

Value

does not return a value, plots are printed to device or files depending on output_dir parameter

Examples

# note; this example downloads data when first run, and typically takes ~60seconds

# store the downloaded files in the following directory. Here, the temporary file
# directory is used. Alternatively, consider storing this data in a more permanent location.
# e.g. output_dir="~/data/goat" on unix systems or output_dir="C:/data/goat" on Windows
output_dir = tempdir()

## first run the default example from test_genesets() to obtain geneset results
datasets = download_goat_manuscript_data(output_dir)
genelist = datasets$`Wingo 2020:mass-spec:PMID32424284`
genesets_asis = download_genesets_goatrepo(output_dir)
genesets_filtered = filter_genesets(genesets_asis, genelist)
result = test_genesets(genesets_filtered, genelist, method = "goat",
  score_type = "effectsize", padj_method = "bonferroni", padj_cutoff = 0.05)

# prior to running this function, cluster the genesets
clusters = cluster_genesets(result, genelist)

# use the plot heatmap function and try various color palettes
plot_heatmap(clusters, output_dir, colors = hcl.colors(100, "Viridis", rev = FALSE))
plot_heatmap(clusters, output_dir, colors = hcl.colors(100, "Inferno", rev = FALSE))
plot_heatmap(clusters, output_dir, colors = hcl.colors(100, "Lajolla", rev = TRUE))
plot_heatmap(clusters, output_dir, colors = hcl.colors(100, "Mako", rev = FALSE))
plot_heatmap(clusters, output_dir, colors = hcl.colors(100, "Turku", rev = TRUE))
plot_heatmap(clusters, output_dir, colors = hcl.colors(100, "Grays", rev = TRUE))

Description

Lollipop chart or barplot visualization of geneset enrichment testing results

Usage

plot_lollipop(
  x,
  output_dir,
  only_reduced = FALSE,
  plot_type = "lollipop",
  show_pvalue = FALSE,
  score_xaxis = "minlogp",
  score_color = ifelse(is.data.frame(x) && "score_type" %in% colnames(x) &&
    is.character(x$score_type) && any(x$score_type %in% c("effectsize_up",
    "effectsize_down")), "updown", "minlogp"),
  score_color_limits = "source",
  score_color_updown = c("#E57373", "#5E7ABC"),
  max_ngenes = NA,
  topn = NA,
  padj_cutoff = NA
)

Arguments

Value

if output_dir is NA, a list of ggplot2 objects. Otherwise, write plots to file and do not return any value

Examples

# note; this example downloads data when first run, and typically takes ~60seconds

# store the downloaded files in the following directory. Here, the temporary file
# directory is used. Alternatively, consider storing this data in a more permanent location.
# e.g. output_dir="~/data/go" on unix systems or output_dir="C:/data/go" on Windows
output_dir = tempdir()

## first run the default example from test_genesets() to obtain geneset results
datasets = download_goat_manuscript_data(output_dir)
genelist = datasets$`Wingo 2020:mass-spec:PMID32424284`
genesets_asis = download_genesets_goatrepo(output_dir)
genesets_filtered = filter_genesets(genesets_asis, genelist)
result = test_genesets(genesets_filtered, genelist, method = "goat",
  score_type = "effectsize", padj_method = "bonferroni", padj_cutoff = 0.05)

# generate lollipop charts for each GO domain (CC/BP/MF), with geneset -log10
# adjusted p-value on the x-axis and color-coding by geneset up/down-regulation
plot_lollipop(
  result, output_dir, plot_type = "lollipop", topn = 50,
  score_xaxis = "minlogp", score_color = "updown"
)

# alternatively, as a barplot
plot_lollipop(
  result, output_dir, plot_type = "barplot", topn = 50,
  score_xaxis = "minlogp", score_color = "updown"
)

# alternatively, color-code genesets by enrichment of significant genes using
# parameter `score_color="oddsratio"` . See further `score_geneset_oddsratio`
# function documentation for definition/computation of this score.
plot_lollipop(
  result, output_dir, plot_type = "lollipop", topn = 50,
  score_xaxis = "minlogp", score_color = "oddsratio"
)

Description

plot geneset distance matrix as a network

Usage

plot_network(
  clusters,
  src,
  show_clusters = TRUE,
  show_text = FALSE,
  topn_edges = 5,
  clr_default = "#29b6f6",
  node_color_palette = goat::gg_color_hue
)

Arguments

Value

a ggplot2 object

Description

For each provided geneset, a volcano plot of all genelist log2fc and p-values with respective geneset constituents highlighted

Usage

plot_volcano(
  x,
  genelist,
  plot = TRUE,
  topn_labels = 0,
  color_default = "#B0B0B0",
  color_highlight = "#ef5350",
  color_label = "#000000",
  pointsize = 2,
  pointalpha = 0.75,
  labelsize = 7
)

Arguments

Value

if plot==FALSE, a list of ggplot2 objects. Otherwise, does not return any value

Examples

# note; this example downloads data when first run, and typically takes ~60seconds

# store the downloaded files in the following directory. Here, the temporary file
# directory is used. Alternatively, consider storing this data in a more permanent location.
# e.g. output_dir="~/data/goat" on unix systems or output_dir="C:/data/goat" on Windows
output_dir = tempdir()

## first run the default example from test_genesets() to obtain geneset results
datasets = download_goat_manuscript_data(output_dir)
genelist = datasets$`Wingo 2020:mass-spec:PMID32424284`
genesets_asis = download_genesets_goatrepo(output_dir)
genesets_filtered = filter_genesets(genesets_asis, genelist)
result = test_genesets(genesets_filtered, genelist, method = "goat",
  score_type = "effectsize", padj_method = "bonferroni", padj_cutoff = 0.05)

## example 1; select top10 GO CC terms from the geneset testing results
result_subset = result |> filter(source == "GO_CC") |> arrange(pvalue) |> head(n = 10)
pdf(paste0(output_dir, "/volcano_CC_top10.pdf"), width = 4, height = 4)
plot_volcano(result_subset, genelist)
dev.off()

## example 2;, select small genesets that are significant and have
## near-exclusive enrichment in either up up/down-regulated genes
# first, add geneset directionality scores to our results
result = score_geneset_directionality(result, genelist)
# next, subset the geneset result table
result_subset = result |>
  filter(signif & ngenes <= 50 & abs(score_directionality_rank) > 0.6) |>
  arrange(pvalue_adjust)
# finally, create plots. Note that the genelist contains a column 'symbol'
# which we use here to print labels for the topN genes per plotted geneset
pdf(paste0(output_dir, "/volcano_signif_ngenes50_directionality06.pdf"), width = 4, height = 4)
plot_volcano(result_subset, genelist, topn_labels = 10)
dev.off()

Description

rank^2 can yield huge values, e.g. a large genelist of N=50000 genes would imply max gene score = 50000^2 = 2.5e+09. Thus we rescale these scores between 0~1000 so downstream applications (e.g. sum of 10000 gene scores) don't explode to huge numbers.

Usage

rankscore(x, sort1, sort2, sort3, colname)

Arguments

Value

data.frame x with added column colname, containing gene scores between 0 and 1000

Examples

x = data.frame(gene=c(1, 2, 3, 4, 5),
                 pvalue=c(0.01, 1, 1, 0.1, 1),
                 effectsize=c(-2, 0.25, 0.5, 1, 0.25))
  print(x, row.names = FALSE)
  print(rankscore(x, sort1 = -1*x$pvalue, sort2 = abs(x$effectsize),
        sort3 = x$gene, colname="score") |>
    arrange(desc(score)), row.names = FALSE)

Description

Gene score array, from low to high scores

Usage

rankscore_fixed_order(n)

Arguments

Value

result of (1:n)^2 / (n^2 / 1000)

Description

Analyses are performed independently per 'source' of genesets. The result of this function is the geneset table with a newly appended column 'signif_and_reduced'

Usage

reduce_genesets(
  clusters,
  simscore_threshold = 0.9,
  universe_fraction = 0.25,
  signifgenes_fraction = 0.9
)

Arguments

Value

the genesets table from the clusters parameter, with results in column "signif_and_reduced"

Description

Works for any filtered geneset table generated by this package, e.g. results from any of these functions; filter_genesets(), test_genesets() or simplify_genesets().

The genelist table is required to 1) lookup gene symbols and 2) determine the ordering of genes within-geneset. The latter makes it such that the output table shows the most important genes (in context of user's data) first.

All output columns that list gene identifiers or symbols are trimmed to 25000 characters. Depending on your input gene identifier type and dataset filtering settings (e.g. max geneset size), this might lead to some truncation. For example, 15 characters per ensembl gene ID (if these are used in provided genesets instead of default NCBI Entrez) and 1 char as delimiter, only the top ~1500 ensembl gene identifiers are included. If your geneset filtering allows for larger genesets and you want to check if any gene identifiers are lost upon saving, keep an eye out for the trailing ellipsis ('...'). Because genes in each geneset of the output table are sorted by order of importance as indicated in the genelist table, this shouldn't be an issue in typical use-cases where the output table is to check the topN most important/significant genes of a geneset of interest.

Usage

save_genesets(x, genelist, filename, arrange_genes = TRUE)

Arguments

Value

if filename is NA, returns the validated and formatted geneset result table. Otherwise, writes the table to file and does not return a value

Description

Importantly, the scope/utility of this score is limited to help users post-hoc filter for genesets that contain mostly up/down-regulated genes. However, this might not coincide with the geneset pvalues / significance. For example, genesets may exclusively contain genes with a positive effectsize but at the same time these can all be minor effects/values and thus the geneset is not significant or less significant than other genesets with the exact same "directionality score". For example, genesets may contain both up- and down-regulated genes but still be significant when testing with GOAT and using score_type='effectsize'

The scores computed with this function can help in post-hoc interpretation of GOAT results to further narrow down all significant genesets to a subset with strong directionality. For example, after test_genesets() we can filter the results for A) significant genesets and B) that contain at most N genes and C) that are near-exclusively up/down-regulated. Bringing this all together (also useful for other types of geneset testing, like ORA, score_type="pvalue", etc); result = test_genesets(genesets_filtered, genelist, method = "goat", score_type = "effectsize", padj_method = "bonferroni", padj_cutoff = 0.05) result = score_geneset_directionality(result, genelist) result |> filter(signif & ngenes <= 100 & abs(score_directionality_rank) >= 0.95)

Usage

score_geneset_directionality(genesets, genelist)

Arguments

Value

input genesets with results in 3 columns; score_directionality_rank is the weighted gene score where gene values are the sign of their effectsize and weights are linearly proportional to their inverse ranks in the input genelist. score_directionality_rank2 is similar, but now using rank^2 gene weights to boost the influence of most-important genes in the input genelist. score_directionality_value uses the absolute gene effectsizes as gene weights Note that the latter is least robust as it depends on the scaling of input data!

score definitions;

geneset directionality score = weighted mean of respective genes, where gene weights are 1 minus their relative rank in up/down-regulation (depending on negative/positive effectsize) and the value for each gene is -1 or 1 depending on up/down-regulation (sign of effectsize).

pseudocode;

1. gene values and weights A) gene weight between 0 and 1 for the subset of upregulated genes / positive effectsizes;

• r = for the subset of genes with effectsize > 0, compute gene rank (1 = highest effectsize, N = smallest effectsize that is greater than zero)

• weight = 1 - r/max(r) B) analogous to (A) for the subset of genes with negative effectsize C) result per gene: value = sign of its effectsize, weight = 0 if effectsize 0, otherwise respective weights from (A) or (B)

1. geneset score_directionality = weighted mean over values/weights of respective genes

Description

gs_signif = number of significant genes in geneset G that intersect with user's genelist (i.e. foreground genes in G) gs_all = number of genes in geneset G that intersect with user's genelist (i.e. foreground+background genes in G) k_signif = total number of significant genes in user's genelist k_all = total number of genes in user's genelist

gs_signif/gs_all = ratio of foreground genes in geneset G k_signif/k_all = ratio of overall foreground genes (i.e. expected value for a random geneset)

oddsratio = (gs_signif/gs_all) / (k_signif/k_all)

Usage

score_geneset_oddsratio(genesets, genelist)

Arguments

Value

input genesets with results in column "score_oddsratio"

Description

replacement for stringr::trunc() so we don't need a package dependency for just 1 function (our code was adapter therefrom)

Usage

string_trunc_right(string, width, trim_left = FALSE)

Arguments

Value

truncated variant of input string

Description

Map the the symbol column in a table to HGNC human gene IDs by matching official gene symbols and synonyms

Usage

symbol_to_entrez(x, hgnc)

Arguments

Value

entrez gene IDs are returned in the "gene" column of table x. Additionally, columns "entrez_id", "hgnc_id" and "hgnc_symbol"

Examples

# TODO: update the filename to your downloaded file
  # download instructions in the documentation of `hgnc_idmap_table()`
  f = "C:/DATA/HGNC/hgnc_complete_set.txt"

  if(file.exists(f)) {
    df = data.frame(symbol = c("vamp2", "STXBP1", "UNC18", NA, "PSD95", "NOT-A-GENE"))
    hgnc = hgnc_idmap_table(f)
    df = symbol_to_entrez(df, hgnc)
    print(df)
  }

Description

Perform geneset enrichment testing using any supported method

Usage

test_genesets(
  genesets,
  genelist,
  method,
  padj_method = "BH",
  padj_sources = TRUE,
  padj_cutoff = 0.01,
  padj_min_signifgenes = 0L,
  ...
)

Arguments

Details

After application of the enrichment testing algorithm (e.g. GOAT, ORA or GSEA), multiple testing correction is applied to obtain adjusted p-values using padjust_genesets. That function will first apply the specified pvalue adjustment procedure in the padj_method parameter within each 'source' in the genesets table. Second, it applies Bonferroni adjustment to all p-values according to the number of different geneset sources that were tested (or set padj_sources = FALSE to disable).

For example, if the input is a genesets table that contains GO_CC, GO_BP and GO_MF genesets, first multiple testing correction is applied within each source (e.g. using FDR if so desired) and afterwards a Bonferroni correction is applied based on 3 repeated analyses.

Note that this is more rigorous than typical GO tools; hypothetically, one could split all GO_CC pathways into 1000 different databases/'sources' and then run enrichment testing. Consequently, the multiple testing burden is reduced if one doesn't adjust p-values for the number of 'sources' as we do here.

Value

the input genesets, with results stored in columns 'pvalue', 'pvalue_adjust' and 'signif'

Examples

#' # note; this example downloads data when first run, and typically takes ~60seconds

## Basic example for a complete GOAT workflow
# Downloads test data to your computer and stores it at current working directory
# Refer to the GitHub documentation for elaborate documentation and a worked example

# store the downloaded files in the following directory. Here, the temporary file
# directory is used. Alternatively, consider storing this data in a more permanent location.
# e.g. output_dir="~/data/go" on unix systems or output_dir="C:/data/go" on Windows
output_dir = tempdir()

# download an example gene list
datasets = download_goat_manuscript_data(output_dir)
genelist = datasets$`Wingo 2020:mass-spec:PMID32424284`

# download GO genesets
genesets_asis = download_genesets_goatrepo(output_dir)

# filter genesets for sufficient overlap with the genelist, then apply GOAT
genesets_filtered = filter_genesets(genesets_asis, genelist)
result = test_genesets(genesets_filtered, genelist, method = "goat",
  score_type = "effectsize", padj_method = "bonferroni", padj_cutoff = 0.05)

# print first 10 rows of the result table
print(result |> select(source, name, ngenes, pvalue_adjust) |> utils::head(n=10))

Description

In most cases, it's more convenient to call the more generic test_genesets function which also deals with multiple-testing correction (per geneset source)

It is assumed that the genesets and genelist parameters are in sync, i.e. genesets provided here is the result of the filter_genesets() function (using the same genelist table)

Same as hypergeometric for non-EASE, but slower because stats::fisher.test isn't vectorized

Only genesets with at least 1 significant gene are subjected to statistical testing (e.g. NA is returned for genesets without significant genes)

Usage

test_genesets_fisherexact(
  genesets,
  genelist,
  require_nsignif = 1L,
  use_ease = FALSE
)

Arguments

Value

input genesets table with results in the "pvalue" column

See Also

test_genesets

Examples

# note; this example downloads data when first run, and typically takes ~60seconds

# store the downloaded files in the following directory. Here, the temporary file
# directory is used. Alternatively, consider storing this data in a more permanent location.
# e.g. output_dir="~/data/goat" on unix systems or output_dir="C:/data/goat" on Windows
output_dir = tempdir()

## first run the default example from test_genesets() to obtain input data
datasets = download_goat_manuscript_data(output_dir)
genelist = datasets$`Wingo 2020:mass-spec:PMID32424284`
genesets_asis = download_genesets_goatrepo(output_dir)
genesets_filtered = filter_genesets(genesets_asis, genelist)

## example: same results between Fisher-exact and hypergeometric tests
result_hg = test_genesets_hypergeometric(genesets_filtered, genelist, require_nsignif = 3L)
result_fe = test_genesets_fisherexact(genesets_filtered, genelist, require_nsignif = 3L)
all.equal(result_hg$pvalue, result_fe$pvalue)

Description

In typical use-cases, one applies test_genesets() instead with parameter method="goat" , which in turn will use test_genesets_goat_precomputed (and not this function).

Usage

test_genesets_goat_bootstrap(
  genesets,
  genelist,
  score_type,
  niter = 5e+05,
  verbose = FALSE
)

Arguments

Value

see test_genesets_goat_precomputed

Description

In typical use-cases, one applies test_genesets() instead with parameter method="goat" , which in turn will use test_genesets_goat_precomputed (and not this function).

Usage

test_genesets_goat_fitfunction(
  genesets,
  genelist,
  score_type,
  niter = 5e+05,
  verbose = FALSE
)

Arguments

Details

Optionally, use this function to ignore precomputed/bundled null distribution estimates and perform new permutations and function fitting (e.g. because you want to test the effect of huge niter parameter, but beware of RAM requirements)

Value

see test_genesets_goat_precomputed

Description

In most cases, it's more convenient to call the more generic test_genesets function which also applies multiple-testing correction (per geneset source) to the geneset p-values computed by this function.

This is the canonical geneset test function for GOAT that uses precomputed null distributions that are bundled with the GOAT package

Usage

test_genesets_goat_precomputed(genesets, genelist, score_type)

Arguments

Value

input genesets table with results in the "pvalue", "score_type" columns. "zscore" column: A standardized z-score is computed from geneset p-values + effectsize direction (up/down) if tested. Importantly, we here return standardized z-scores because the GOAT geneset score (mean of gene scores) is relative to the respective geneset-size-matched null distributions (a skewed normal)! In contrast, the standardized z-scores are comparable between genesets (as are the pvalues obviously).

Only if either (or both) the effectsize-up/down was tested, the direction of regulation has been tested (effectsize_abs and pvalue score types are agnositic to up/down regulation). So when score_type was set to any of effectsize/effectsize_down/effectsize_up, the z-scores are negative values in case the "score_type" output column is "effectsize_down".

See Also

test_genesets

Examples

# note; this example downloads data when first run, and typically takes ~60seconds

# store the downloaded files in the following directory. Here, the temporary file
# directory is used. Alternatively, consider storing this data in a more permanent location.
# e.g. output_dir="~/data/goat" on unix systems or output_dir="C:/data/goat" on Windows
output_dir = tempdir()

## first run the default example from test_genesets() to obtain input data
datasets = download_goat_manuscript_data(output_dir)
genelist = datasets$`Wingo 2020:mass-spec:PMID32424284`
genesets_asis = download_genesets_goatrepo(output_dir)
genesets_filtered = filter_genesets(genesets_asis, genelist)

### we here compare GOAT with precomputed null distributions against
### a GOAT function that performs bootstrapping to compute null distributions on-demand

# apply goat with precomputed null (default) and goat with on-demand bootstrapping
result_precomputed = test_genesets(genesets_filtered, genelist, method = "goat",
  score_type = "effectsize", padj_method = "bonferroni", padj_cutoff = 0.05) |>
  # undo sorting by p-value @ test_genesets(), instead sort by stable IDs
  arrange(source, id)
result_bootstrapped = test_genesets(genesets_filtered, genelist, method = "goat_bootstrap",
  score_type = "effectsize", padj_method = "bonferroni", padj_cutoff = 0.05, verbose = TRUE) |>
  arrange(source, id)

# tables should align
stopifnot(result_precomputed$id == result_bootstrapped$id)
# no missing values
stopifnot(is.finite(result_precomputed$pvalue) &
          is.finite(is.finite(result_bootstrapped$pvalue)))

# compare results
plot(result_precomputed$pvalue, result_bootstrapped$pvalue)
abline(0, 1, col=2)

plot(minlog10_fixzero(result_precomputed$pvalue),
     minlog10_fixzero(result_bootstrapped$pvalue))
abline(0, 1, col=2)

summary(minlog10_fixzero(result_precomputed$pvalue) -
        minlog10_fixzero(result_bootstrapped$pvalue))

Description

In most cases, it's more convenient to call the more generic test_genesets function which also deals with multiple-testing correction (per geneset source)

Usage

test_genesets_gsea(
  genesets,
  genelist,
  score_type = NULL,
  parallel_threads = 1L,
  gseaParam = 1,
  nPermSimple = 50000,
  gsea_genelist_col = NULL,
  gsea_scoretype = NULL,
  random_seed = 123
)

Arguments

Value

input genesets table with results in the "pvalue", "score_type" and "gsea_nes" columns

See Also

test_genesets

Description

In most cases, it's more convenient to call the more generic test_genesets function which also deals with multiple-testing correction (per geneset source)

It is assumed that the genesets and genelist parameters are in sync, i.e. genesets provided here is the result of the filter_genesets() function (using the same genelist table)

Only genesets with at least 1 significant gene are subjected to statistical testing (e.g. NA is returned for genesets without significant genes)

Usage

test_genesets_hypergeometric(genesets, genelist, require_nsignif = 1L)

Arguments

Value

input genesets table with results in the "pvalue" column

See Also

test_genesets

Description

refer to the goat::treemap_plot() function for a complete example

Usage

treemap_data(
  geneset_ids,
  genesets,
  genesets_test_result,
  simplify = "none",
  toplevel_max_ngenes = NA
)

Arguments

Value

data structure needed for treemap_plot()

Description

simple wrapper around the treemap R package. To customize this plot, copy/paste its code and tweak parameters as desired

Usage

treemap_plot(x, label_group = FALSE)

Arguments

Value

a ggplot2 object constructed by treemap::treemap()

Examples

# note; this example downloads data when first run, and typically takes ~60seconds

# store the downloaded files in the following directory. Here, the temporary file
# directory is used. Alternatively, consider storing this data in a more permanent location.
# e.g. output_dir="~/data/goat" on unix systems or output_dir="C:/data/goat" on Windows
output_dir = tempdir()

## first run the default example from test_genesets() to obtain geneset results
datasets = download_goat_manuscript_data(output_dir)
genelist = datasets$`Wingo 2020:mass-spec:PMID32424284`
genesets_asis = download_genesets_goatrepo(output_dir)
genesets_filtered = filter_genesets(genesets_asis, genelist)
result = test_genesets(genesets_filtered, genelist, method = "goat",
  score_type = "effectsize", padj_method = "bonferroni", padj_cutoff = 0.05)

# subset GO CC results
x = result |> filter(signif & source == "GO_CC")
tm = treemap_data(
  geneset_ids = x$id,
  genesets = genesets_filtered,
  genesets_test_result = x,
  simplify = "leaf_only" # options: none/leaf_only/prune_singletons/pvalue
)
treemap_plot(tm$treemap_plotdata)