This is a long post; I am trying to find out why R's GO.db gives different results than Python's go_db https://github.com/endrebak/go_db.
go_db seems to give different results than the R package GO.db.
go_db computes ancestors by looking upwards in the tree; so that C is an ancestor of A if A is a C, A is part of C or C has part A or any such relation exists transitively.
Since go_db is a dinky script I threw together quickly and GO.db is the seemingly canonical Gene Ontology package I am guessing I am in the wrong, but if anyone can explain why, I would appreciate it. It is probably my usage of GO.db that is incorrect, I'd wager.
To run the scripts below you need to install wide_diapeR and go_db with
pip install -U go_db wide_diaper
and have the R package GO.db installed
from __future__ import print_function
from widediaper import R
from go_db import get_offspring
r = R()
r.load_library("GO.db")
r("cc_offspring = as.data.frame(GOCCOFFSPRING)")
r_cc_offspring = r.get("cc_offspring")
r_cc_offspring.columns = ["Offspring", "Parent"]
py_cc_offspring = get_offspring("CC")
py_go_0001533 = py_cc_offspring.loc[py_cc_offspring.Offspring == "GO:0001533"]
r_go_0001533 = r_cc_offspring.loc[r_cc_offspring.Offspring == "GO:0001533"]
# Parents in GO.db, but not go_db
set(r_go_0001533.Parent) - set(py_go_0001533.Parent)
# {'GO:0005622',
# 'GO:0005856',
# 'GO:0043226',
# 'GO:0043228',
# 'GO:0043229',
# 'GO:0043232',
# 'GO:0044424'}
# Parents in go_db, but not GO.db
set(py_go_0001533.Parent) - set(r_go_0001533.Parent)
# {'GO:0005886',
# 'GO:0012505',
# 'GO:0016020',
# 'GO:0030313',
# 'GO:0031975',
# 'GO:0071944',
# 'GO:0097653'}
As you can see above, there are great differences in what nodes the two packages consider to be ancestors of "GO:0001533".
Let us look at the go.obo file and see if we can make sense of this. The script below prints out all the ancestors of a node (including the node itself):
# download obo file from http://purl.obolibrary.org/obo/go.obo
obo_file = '/Users/labsenter/.go_db/go.obo' # Use the path to whereever you are storing your .obo file
obo_records = "".join(open(obo_file).readlines()).split("\n\n[Term]")[1:-1]
py_ancestors = py_cc_offspring.loc[py_cc_offspring.Offspring == "GO:0001533"].Parent
py_ancestors.append("GO:0001533")
py_ancestor_search_terms = ["id: " + ancestor for ancestor in py_ancestors]
r_ancestors = list(py_cc_offspring.loc[py_cc_offspring.Offspring == "GO:0001533"].Parent)
r_ancestors.append("GO:0001533")
r_ancestor_search_terms = ["id: " + ancestor for ancestor in r_ancestors]
def _print_records(search_terms):
for record in obo_records:
for search_term in search_terms:
if search_term in record:
for line in record.split("\n"):
for term in ["id:", "is_a:", "relationship:", "intersection_of:"]:
if term in line:
print(line)
print()
print("### Py ancestors\n")
_print_records(py_ancestor_search_terms)
print("### R ancestors\n")
_print_records(r_ancestor_search_terms)
The script above gives the following results:
### Py ancestors
id: GO:0001533
is_a: GO:0005886 ! plasma membrane
id: GO:0005575
alt_id: GO:0008372
id: GO:0005623
is_a: GO:0005575 ! cellular_component
id: GO:0005886
alt_id: GO:0005904
is_a: GO:0016020 ! membrane
is_a: GO:0044464 ! cell part
relationship: part_of GO:0071944 ! cell periphery
id: GO:0012505
is_a: GO:0044464 ! cell part
relationship: has_part GO:0005773 ! vacuole
relationship: has_part GO:0005886 ! plasma membrane
id: GO:0016020
is_a: GO:0005575 ! cellular_component
id: GO:0030313
is_a: GO:0031975 ! envelope
relationship: has_part GO:0005886 ! plasma membrane
id: GO:0031975
is_a: GO:0044464 ! cell part
id: GO:0044464
is_a: GO:0005575 ! cellular_component
intersection_of: GO:0005575 ! cellular_component
intersection_of: part_of GO:0005623 ! cell
relationship: part_of GO:0005623 ! cell
id: GO:0071944
is_a: GO:0044464 ! cell part
id: GO:0097653
is_a: GO:0044464 ! cell part
relationship: has_part GO:0005622 ! intracellular
relationship: has_part GO:0005886 ! plasma membrane
### R ancestors
id: GO:0001533
is_a: GO:0005886 ! plasma membrane
id: GO:0005575
alt_id: GO:0008372
id: GO:0005623
is_a: GO:0005575 ! cellular_component
id: GO:0005886
alt_id: GO:0005904
is_a: GO:0016020 ! membrane
is_a: GO:0044464 ! cell part
relationship: part_of GO:0071944 ! cell periphery
id: GO:0012505
is_a: GO:0044464 ! cell part
relationship: has_part GO:0005773 ! vacuole
relationship: has_part GO:0005886 ! plasma membrane
id: GO:0016020
is_a: GO:0005575 ! cellular_component
id: GO:0030313
is_a: GO:0031975 ! envelope
relationship: has_part GO:0005886 ! plasma membrane
id: GO:0031975
is_a: GO:0044464 ! cell part
id: GO:0044464
is_a: GO:0005575 ! cellular_component
intersection_of: GO:0005575 ! cellular_component
intersection_of: part_of GO:0005623 ! cell
relationship: part_of GO:0005623 ! cell
id: GO:0071944
is_a: GO:0044464 ! cell part
id: GO:0097653
is_a: GO:0044464 ! cell part
relationship: has_part GO:0005622 ! intracellular
relationship: has_part GO:0005886 ! plasma membrane
Since you can see that GO:0001533 is a GO:0005886 which again is a GO:0016020, GO:0016020 should be considered an ancestor of GO:0001533. Unfortunately, that is not the case in GO.db. If anyone knows why, please explain!
