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 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:

  3. Limit the GOEA to biological process, molecular function, and/or cellular compartment

  4. with uncorrected pvalues < 0.05 to the screen: (default)

  5. with uncorrected pvalues < 0.05 to the screen, grouped:
    –sections=…

  6. with uncorrected pvalues < 0.05 to files:
    –outfile=goea.xlsx

  7. with Benjamini-Hochberg pvalues < 0.05 (plus bonferroni, sidak, holm) –pval_field=fdr_bh

  8. with Benjamini-Hochberg pvalues < 0.05 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
    –method=fdr_bh –outfile=goea_fdr_bh_grpd.xlsx –sections…

  10. regardless of pvalue (ALL GO terms) –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

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