# Gene Ontolog Enrichment Analyses (GOEA) from the command-line ## GOEA command-line examples 1. [with **pvalues < 0.05 corrected using Benjamini-Hochberg multiple test correction**](#1-enrichment-analysis-using-various-annotation-formats) using annotation in a: * Original id-to-GOs text format * GAF file * GPAD file * NCBI's gene2go file 2. **Exclude or include annotations by evidence code** Use with GAF, GPAD, or NCBI's gene2go files: * [**Exclude annotations inferred from Electronic Annotation (IEA)**](#2a-exclude-annotations-inferred-from-electronic-annotation-iea) * [**Include only annotations inferred from experimental evidence**](#2b-include-only-annotations-inferred-from-experimental-evidence) * [Get **list of evidence codes**](#2c-get-list-of-evidence-codes) 3. [**Limit the GOEA to _biological process_, _molecular function_, and/or _cellular compartment_**](#3-limit-the-goea-to-biological-process-molecular-function-andor-cellular-compartment) 4. [with **uncorrected pvalues < 0.05** to the screen](#4-print-go-terms-with-uncorrected-pvalues--005-to-the-screen-default): (default) 5. [with **uncorrected pvalues < 0.05** to the screen, **grouped**](#5-print-go-terms-with-uncorrected-pvalues--005-to-the-screen-grouped): **--sections=...** 6. [with **uncorrected pvalues < 0.05** to files](#6-print-go-terms-with-uncorrected-pvalues--005-to-an-xlsx-file): **--outfile=goea.xlsx** 7. [with **Benjamini-Hochberg pvalues < 0.05** (plus bonferroni, sidak, holm)](#7-print-go-terms-with-benjamini-hochberg-pvalues--005) **--pval_field=fdr_bh** 8. [with **Benjamini-Hochberg pvalues < 0.05** to an **xlsx** file](#8-print-go-terms-with-benjamini-hochberg-only-pvalues--005-to-an-xlsx-file) **--method=fdr_bh --outfile=goea_fdr_bh_flat.xlsx** 9. [with **Benjamini-Hochberg pvalues < 0.05** to an **xlsx** file, **grouped**](#9-print-go-terms-with-benjamini-hochberg-only-pvalues--005-to-an-xlsx-file-grouped) **--method=fdr_bh --outfile=goea_fdr_bh_grpd.xlsx --sections...** 10. [regardless of pvalue **(ALL GO terms)**](#10-print-all-go-terms-regardless-of-pvalue) **--pval=-1** ## Example Details CMD: goatools find_enrichment data/study data/population data/association --outfile=goea.xlsx,goea.tsv --pval_field=fdr_bh ### 1) Enrichment Analysis using various annotation formats Arguments to print results where Benjamini-Hochberg values are less than 0.05: ``` --pval=0.05 # print pvalue's < 0.05 --method=fdr_bh # Use Benjamini-Hochberg multiple test correction on uncorrected p-values --pval_field=fdr_bh # print fdr_bh values < 0.05 (rather than uncorrected pvalues) --outfile=FILENAME.xlsx # Write to an Excel spreadsheet ``` ### 1A) Original id-to-GOs text format ``` goatools find_enrichment data/study.txt data/population.txt data/association.txt --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_id2gos.xlsx ``` ### 1B) GAF file ``` goatools find_enrichment ids_stu.txt ids_pop.txt goa_human.gaf --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_gaf.xlsx ``` ### 1B) GPAD file ``` goatools find_enrichment ids_stu.txt ids_pop.txt goa_human.gpad --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_gpad.xlsx ``` ### 1D) NCBI's gene2go file #### Human is the default: ``` goatools find_enrichment ids_stu_gene2go_9606.txt ids_pop_gene2go_9606.txt gene2go --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_gene2go_9606.xlsx ``` #### Specify mouse using the **--taxid** argument: ``` goatools find_enrichment ids_stu_gene2go_10090.txt ids_pop_gene2go_10090.txt gene2go --taxid=10090 --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_gene2go_10090.xlsx ``` ### 2) Exclude or include annotations by evidence code #### 2a) Exclude annotations inferred from Electronic Annotation (IEA) ``` EXCLUDE all IEA annotations: --ev_exc=IEA # Evidence codes with IEA are excluded ``` ``` goatools find_enrichment ids_stu_gpad.txt ids_pop_gpad.txt goa_human.gpad --ev_exc=IEA --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_gpad.xlsx ``` #### 2b) Include only annotations inferred from experimental evidence ``` INCLUDE evidence codes by group of by code: 1) --ev_inc=Experimental # By group name 1) --ev_inc=Experimental,Similarity # By group names 2) --ev_inc=EXP,IDA,IPI,IMP,IGI,IEP # List all the Experimental codes ``` **Specify the experimental evidence codes by their group name:** ``` goatools find_enrichment ids_stu_gpad.txt ids_pop_gpad.txt goa_human.gpad --ev_inc=Experimental --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_gpad.xlsx ``` **Specify the experimental evidence codes explicitly:** ``` goatools find_enrichment ids_stu_gpad.txt ids_pop_gpad.txt goa_human.gpad --ev_inc=EXP,IDA,IPI,IMP,IGI,IEP --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_gpad.xlsx ``` #### 2c) Get list of evidence codes The argument, **--ev_help**, will cause the Evidence codes to be printed: ``` goatools find_enrichment --ev_help ``` ##### Evidence Code Help You can use either group names, like Experimental, or codes on the commandline: ``` EVIDENCE CODES: Experimental: EXP Inferred from Experiment IDA Inferred from Direct Assay IPI Inferred from Physical Interaction IMP Inferred from Mutant Phenotype IGI Inferred from Genetic Interaction IEP Inferred from Expression Pattern Similarity: ISS Inferred from Sequence or structural Similarity ISO Inferred from Sequence Orthology ISA Inferred from Sequence Alignment ISM Inferred from Sequence Model used in manual assertion IGC Inferred from Genomic Context IBA Inferred from Biological aspect of Ancestor IBD Inferred from Biological aspect of Descendant IKR Inferred from phylogenetic determination of loss of key residues (manual assertion) IRD Inferred from Rapid Divergence from ancestral sequence (manual assertion) IMR Phylogenetic determination of loss of key residues in manual assertion Combinatorial: RCA Inferred from Reviewed Computational Analysis High_Throughput: HTP Inferred from High Throughput Experimental HDA Inferred from High Throughput Direct Assay HMP Inferred from High Throughput Mutant Phenotype HGI Inferred from High Throughput Genetic Interaction HEP Inferred from High Throughput Expression Pattern Author: TAS Traceable Author Statement used in manual assertion NAS Non-traceable Author Statement used in manual assertion Curatorial: IC Inferred by Curator No biological data: ND No biological Data available Automatic: IEA Inferred from Electronic Annotation ``` ### 3) Limit the GOEA to _biological process_, _molecular function_, and/or _cellular compartment_ Use the --ns option to limit the GOEA to any of: | NS | Namespace | |----|--------------| | BP | Biological Process | | MF | Molecular Function | | CC | Cellular Component | Use the two-letter abbreviation as the value in the --ns option: ``` Namespace examples: --ns=BP # Only run one GOEA on the Biological Process branch --ns=BP,MF # Run two GOEAs: One on Biological Process and one on MOlecular Function --ns=CC # Only run one GOEA on the Molecular Function branch ``` #### Without --ns, all three branches will be analyzed: ``` goatools find_enrichment ids_stu_gene2go_9606.txt ids_pop_gene2go_9606.txt gene2go --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_gene2go_9606.xlsx ``` #### With **--ns=MF**, only the _**molecular function**_ branch will be analyzed: ``` goatools find_enrichment ids_stu_gene2go_9606.txt ids_pop_gene2go_9606.txt gene2go --ns=MF --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_gene2go_9606.xlsx ``` #### With **--ns=BP,MF**, run GOEAs on both the _**biological process**_ and the _**molecular function**_ branches: ``` goatools find_enrichment ids_stu_gene2go_9606.txt ids_pop_gene2go_9606.txt gene2go --ns=BP,MF --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_gene2go_9606.xlsx ``` #### With **--ns=MF** and **--ev_inc=IPI**: Run a GOEA on the _**molecular function**_ branch for just the evidence code, **IPI**. **IPI** includes evidence codes which are inferred from Physical Interaction ``` goatools find_enrichment ids_stu_gene2go_9606.txt ids_pop_gene2go_9606.txt gene2go --ns=MF --ev_inc=IPI --pval=0.05 --method=fdr_bh --pval_field=fdr_bh --outfile=results_gene2go_9606.xlsx ``` ### 4) Print GO terms with uncorrected pvalues < 0.05 to the screen (default) The default is to print all GO terms with uncorrected P-values < 0.05 (default) ``` goatools find_enrichment data/study.txt data/population.txt data/association.txt ``` |GO|NS|ep|name|study|pop|p_uncorrected|depth|# study|p_bonferroni|p_sidak|p_holm|p_fdr_bh|study_items |--|--|--|----|-----|---|-------------|-----|-----------|------------|-------|------|--------|----------- |GO:0036211|BP|e|protein modification process|33/276|1725/33239|8.070e-06|5|33|0.0491|0.04789|0.04910|0.00982|AT1G13580... |GO:0006464|BP|e|cellular protein modification process|33/276|1725/33239|8.070e-06|6|33|0.0491|0.04789|0.04910|0.00982|AT1G13580... ### 5) Print GO terms with uncorrected pvalues < 0.05 to the screen, grouped optional attribute: **--sections=goatools.test_data.sections.data2018_07_find_enrichment** Acceptable values for **--sections** are: | File type | Example value |---------------|------------------------------------------------------------ | Python module | goatools.test_data.sections.data2018_07_find_enrichment | text file | sections.txt | Python file | goatools/test_data/sections/data2018_07_find_enrichment.py ``` goatools find_enrichment data/study.txt data/population.txt data/association.txt --sections=goatools.test_data.sections.data2018_07_find_enrichment protein GO:0036211 BP e 9.83e-03 1320 5.14 D05 33/276 1725/33239 protein modification process ... GO:0006464 BP e 9.83e-03 1318 5.14 D06 33/276 1725/33239 cellular protein modification process ... ... phosph GO:0016301 MF e 4.64e-02 355 4.87 D04 25/276 1310/33239 kinase activity ... ... ``` ### 6) Print GO terms with uncorrected pvalues < 0.05 to an xlsx file optional attribute: **--outfile=goea_uncorr.xlsx** ``` $ goatools find_enrichment data/study.txt data/population.txt data/association.txt --outfile=goea_uncorr.xlsx ... 253 items WROTE: goea.xlsx ``` ### 7) Print GO terms with Benjamini-Hochberg pvalues < 0.05 optional attribute: **--pval_field=fdr_bh** optional attribute: **----outfile=goea_fdr_bh.xlsx,goea_fdr_bh.tsv** ``` $ goatools find_enrichment data/study.txt data/population.txt data/association.txt --outfile=goea_fdr_bh.xlsx,goea_fdr_bh.tsv --pval_field=fdr_bh 17 items WROTE: goea.xlsx 17 items WROTE: goea.tsv ``` ### 8) Print GO terms with Benjamini-Hochberg (only) pvalues < 0.05 to an xlsx file optional attribute: **--method=fdr_bh** optional attribute: **--outfile=goea_fdr_bh_flat.xlsx** ``` $ goatools find_enrichment data/study.txt data/population.txt data/association.txt --outfile=goea_fdr_bh_flat.xlsx --method=fdr_bh 17 items WROTE: goea_fdr_bh_flat.xlsx ``` ### 9) Print GO terms with Benjamini-Hochberg (only) pvalues < 0.05 to an xlsx file, grouped optional attribute: **--method=fdr_bh** optional attribute: **--outfile=goea_fdr_bh_grpd.xlsx** optional attribute: **--sections=goatools.test_data.sections.data2018_07_find_enrichment** ``` $ goatools find_enrichment data/study.txt data/population.txt data/association.txt --outfile=goea_fdr_bh_grpd.xlsx --method=fdr_bh --sections=goatools.test_data.sections.data2018_07_find_enrichment 17 items WROTE: goea_fdr_bh_grpd.xlsx ``` ### 10) Print all GO terms, regardless of pvalue optional attribute: **--pval=-1** ``` $ goatools find_enrichment data/study.txt data/population.txt data/association.txt --outfile=goea_all.xlsx,goea_all.tsv --pval=-1 6088 items WROTE: goea.xlsx 6088 items WROTE: goea.tsv ``` ### 11) Decorrelate the GO graph with the *elim* algorithm optional attribute: **--algorithm=elim** By default every GO term is scored independently (`--algorithm=classic`). Because annotations are propagated up the GO DAG, a broad parent term inherits the genes of its more specific children, and so tends to look enriched whenever a child is -- producing long result lists of overlapping, redundant terms. The **elim** algorithm ([Alexa 2006](https://doi.org/10.1093/bioinformatics/btl140)) addresses this. It walks the DAG bottom-up, and whenever a term is significant it removes that term's genes from *all* of the term's ancestors before they are scored. A parent therefore has to be enriched for reasons of its own, not merely by inheritance. ``` $ goatools find_enrichment data/study.txt data/population.txt data/association.txt --algorithm=elim ``` Notes: * `--algorithm` is unrelated to `--method`, which chooses the multiple-testing correction. * `--elim_cutoff` (default **0.01**) is the p-value at or below which a term's genes are eliminated from its ancestors. This matches topGO. The original paper divides the cutoff by the number of scored terms; enable that with `--elim_bonferroni`. * elim requires a one-sided Fisher test (`--alternative=greater`), which it selects for you; asking for another alternative is an error. Under a two-sided test a significantly *depleted* term would eliminate its genes too, inverting the algorithm. This also means elim p-values are not directly comparable to `classic` ones, which are two-sided by default. * elim scores each term conditioned on its neighbours, so its p-values are not independent. topGO's documentation suggests reading them as already accounting for multiple testing; goatools warns if you apply a correction on top of them. Threshold the raw p-value (`--pval_field=uncorrected`) to follow topGO's practice -- stacking a correction is markedly more conservative. * elim requires propagated counts, so it is incompatible with `--no_propagate_counts`. Given the same GO DAG, propagated annotations and gene universe, this implementation reproduces topGO's `elim` scores to floating-point precision. Results will not match a stock topGO run out of the box: goatools propagates over `is_a` only unless you pass `--relationships`, while topGO's `GO.db` graph also uses `part_of` and `regulates`. Copyright (C) 2010-present, DV Klopfenstein, Haibao Tang, et al. All rights reserved.