Do enrichment analyses using the Human Phenotype Ontology (HPO)¶

HPO files located in goatools/notebooks/data/hpo

fin_study = 'data/hpo/genes.list'         # Study genes
fin_pop   = 'data/hpo/gobackground.list'  # Population genes
fin_obo   = 'data/hpo/hp.obo'             # DAG containing HPO terms
fin_anno  = 'data/hpo/hpo.annotation.tab' # Annotation of genes-to-HPO terms

1) Read list of study genes and population genes¶

from goatools.utils import read_geneset

study_ids = read_geneset(fin_study)
population_ids = read_geneset(fin_pop)

    2,252 READ: data/hpo/genes.list
   14,446 READ: data/hpo/gobackground.list

2) Load the human phenotype ontology DAG¶

from goatools.obo_parser import GODag

godag = GODag(fin_obo)

data/hpo/hp.obo: fmt(1.2) rel(hp/2021-02-28) 19,498 Terms

3) Load the annotations of genes to sets of HPO terms¶

from goatools.anno.idtogos_reader import IdToGosReader

annoobj = IdToGosReader(fin_anno, godag=godag)

id2gos = annoobj.get_id2gos()

HMS:0:00:00.594349 187,934 annotations READ: data/hpo/hpo.annotation.tab 
4531 IDs in loaded association branch, human_phenotype

4) Run enrichment analysis on HPO terms¶

from goatools.go_enrichment import GOEnrichmentStudy

goeaobj = GOEnrichmentStudy(
    population_ids,
    annoobj.get_id2gos(),
    godag,
    methods=['bonferroni', 'fdr_bh'],
    pvalcalc='fisher_scipy_stats')

4531 IDs in loaded association branch, human_phenotype

Load  Ontology Enrichment Analysis ...
Propagating term counts up: is_a
 26%  3,823 of 14,446 population items found in association

5) Run an enrichment analysis on HPO terms¶

results = goeaobj.run_study_nts(study_ids)

Runing  Ontology Analysis: current study set of 2252 IDs.
 29%    650 of  2,252 study items found in association
100%  2,252 of  2,252 study items found in population(14446)
Calculating 9,331 uncorrected p-values using fisher_scipy_stats
   9,331 terms are associated with  3,823 of 14,446 population items
   5,802 terms are associated with    650 of  2,252 study items
  METHOD bonferroni:
       2 GO terms found significant (< 0.05=alpha) (  2 enriched +   0 purified): local bonferroni
     489 study items associated with significant GO IDs (enriched)
       0 study items associated with significant GO IDs (purified)
  METHOD fdr_bh:
      12 GO terms found significant (< 0.05=alpha) ( 12 enriched +   0 purified): statsmodels fdr_bh
     498 study items associated with significant GO IDs (enriched)
       0 study items associated with significant GO IDs (purified)

6) Print the results¶

print('namespace       term_id  e/p pval_uncorr Benjamimi/Hochberg Bonferroni  study_ratio population_ratio')
print('--------------- -------- --- ----------- ------------------ ----------  ----------- ----------------')
pat = '{NS} {GO} {e}    {PVAL:8.2e}           {BH:8.2e}   {BONF:8.2e} {RS:>12} {RP:>12}'
for ntd in sorted(results, key=lambda nt: [nt.p_uncorrected, nt.GO]):
    if ntd.p_fdr_bh < 0.05:
        print(pat.format(
            NS=ntd.NS,
            GO=ntd.GO,
            e=ntd.enrichment,
            RS='{}/{}'.format(*ntd.ratio_in_study),
            RP='{}/{}'.format(*ntd.ratio_in_pop),
            PVAL=ntd.p_uncorrected,
            BONF=ntd.p_bonferroni,
            BH=ntd.p_fdr_bh))
print('e: enriched')
print('p: purified')

namespace       term_id  e/p pval_uncorr Benjamimi/Hochberg Bonferroni  study_ratio population_ratio
--------------- -------- --- ----------- ------------------ ----------  ----------- ----------------
human_phenotype HP:0000007 e    3.40e-09           3.18e-05   3.18e-05     459/2252   2322/14446
human_phenotype HP:0000080 e    4.22e-08           1.97e-04   3.94e-04     142/2252    590/14446
human_phenotype HP:0000556 e    7.80e-06           2.26e-02   7.28e-02      74/2252    288/14446
human_phenotype HP:0000512 e    1.21e-05           2.26e-02   1.13e-01      48/2252    166/14446
human_phenotype HP:0000572 e    1.21e-05           2.26e-02   1.13e-01      48/2252    166/14446
human_phenotype HP:0012868 e    2.69e-05           3.85e-02   2.51e-01      10/2252     16/14446
human_phenotype HP:0003138 e    3.65e-05           3.85e-02   3.40e-01       8/2252     11/14446
human_phenotype HP:0031970 e    3.65e-05           3.85e-02   3.40e-01       8/2252     11/14446
human_phenotype HP:0000510 e    4.07e-05           3.85e-02   3.80e-01      49/2252    178/14446
human_phenotype HP:0030453 e    4.12e-05           3.85e-02   3.85e-01      53/2252    197/14446
human_phenotype HP:0032558 e    5.07e-05           3.94e-02   4.73e-01       9/2252     14/14446
human_phenotype HP:0032560 e    5.07e-05           3.94e-02   4.73e-01       9/2252     14/14446
e: enriched
p: purified