Understanding the evolutionary history of a gene family may provide some insight into the mechanisms of its activity. Here we will look at two characteristics for a list of eukaryotic genes of interest:
- Groups that are co-located on a chromosome.
- Groups evolutionarily related through duplication (paralogs).
By expanding on work done in an earlier post, we will identify co-located and duplicated genes using information from the Ensembl genome browser.
Using our previous method of screen scraping with Beautiful Soup, we extend our Ensembl REST-like class to include two new functions. Location is represented as a list of chromosome name, start and end, and parsed from the Gene Summary page:
[sourcecode language="python"]
def location(self, organism, gene_id):
with self._get_open_handle("Gene", "Summary",
organism, gene_id) as in_handle:
soup = BeautifulSoup(in_handle)
loc_tab = soup.find("dd", id="tab_location")
link = loc_tab.find("a")
path, attrs = urllib2.splitattr(link["href"])
for attr in attrs:
if attr.find("r=") == 0:
key, val = attr.split("=")
chrom, location = val.split(":")
start, end = location.split("-")
return chrom, int(start), int(end)
raise ValueError("Did not find location: %s" % link)
[/sourcecode]
Similarly, we can retrieve the details on paralogs from the Ensembl Comparative pages for the genes. With these two functions, a list of the protein IDs of interest, and a dictionary containing species and Ensembl ID references, we collect all of the location and duplication information into python dictionaries. At the same time we maintain a backwards mapping of Ensembl IDs to our original ID:
[sourcecode language="python"]
ensembl_retriever = EnsemblComparaRest(cache_dir)
loc_info = dict()
dup_info = dict()
ensembl_to_uniprot = dict()
for cur_id in all_ids:
cur_rec = db[cur_id]
cur_ensembl_org = cur_rec["org_scientific_name"].replace(" ", "_")
for ensembl_id in cur_rec.get("db_refs_ensembl", []):
paralogs = ensembl_retriever.paralogs(cur_ensembl_org, ensembl_id)
chromosome, start, end = ensembl_retriever.location(
cur_ensembl_org, ensembl_id)
dup_info[ensembl_id] = paralogs
loc_info[ensembl_id] = (cur_rec["org_scientific_name"],
chromosome, start, end)
ensembl_to_uniprot[ensembl_id] = cur_id
[/sourcecode]
Now we want to flatten the paralog dictionary into a list of groups associated by duplication. This is done by using python sets; any groups with shared genes are combined and the resulting unique list is returned:
[sourcecode language="python"]
def examine_paralogs(dup_info, ensembl_to_uniprot):
cur_groups = []
all_base = dup_info.keys()
for base_id, dup_ids in dup_info.items():
overlap = set(dup_ids) & set(all_base)
if len(overlap) > 0:
new_group = set([ensembl_to_uniprot[x] for x in overlap |
set([base_id])])
is_unique = True
for exist_i, exist_group in enumerate(cur_groups):
if len(new_group & exist_group) > 0:
update_group = new_group & exist_group
cur_groups[exist_i] = update_group
is_unique = False
break
if is_unique:
cur_groups.append(new_group)
return [list(g) for g in cur_groups]
[/sourcecode]
We use two functions to similarly group together genes by location. The first function calculates the actual distance given the previously retrieved location information list of (organism, chromosome, start, end). Two items in different organisms or on different chromosomes are placed maximally far apart since they can't be co-located.
[sourcecode language="python"]
def location_distance(loc_one, loc_two):
if loc_one[:2] != loc_two[:2] or loc_one == loc_two:
return sys.maxint
else:
return max(abs(loc_one[3] - loc_two[2]),
abs(loc_two[3] - loc_one[2]))
[/sourcecode]
The next function creates a dictionary of genes co-located based on our threshold, and then uses the examine_paralogs flattening function we defined above to generate groups from this dictionary. The variable distance_thresh defines the distance in which two genes are considered co-located. For the example script, 1 megabase is used, but this can be adjusted according to your own personal definition of close.
[sourcecode language="python"]
def examine_location(loc_info, ensembl_to_uniprot, distance_thresh):
loc_close = collections.defaultdict(lambda: [])
for uniprot_id, loc_one in loc_info.items():
for cmp_id, loc_two in loc_info.items():
if location_distance(loc_one, loc_two) <= distance_thresh:
loc_close[uniprot_id].append(cmp_id)
return examine_paralogs(loc_close, ensembl_to_uniprot)
[/sourcecode]
The full script takes a file input with each line being a uniprot ID and gene name, and also requires a UniProt shelve database like we developed earlier. This shelve database provides the base ID to ensembl ID mappings and organism information, which can be parsed from UniProt XML files.
The resulting co-location and duplication groups are sets of genes which may share an intriguing evolutionary history. The paralogs have been computationally determined to be evolutionarily related based on sequence similarity and the composition of their gene family tree. The co-located genes may be co-selected during evolution, be the result of localized gene duplication obscured by extensive sequence change, or simply be close together based on chance. Further examination of your gene family in light of this information can help determine which of these hypotheses to favor.
comments powered by Disqus