Blue Collar Bioinformatics

Location and duplication information from Ensembl

with 2 comments

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:

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)

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:

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

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:

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 |
            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
            if is_unique:
    return [list(g) for g in cur_groups]

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.

def location_distance(loc_one, loc_two):
    if loc_one[:2] != loc_two[:2] or loc_one == loc_two:
        return sys.maxint
        return max(abs(loc_one[3] - loc_two[2]),
                   abs(loc_two[3] - loc_one[2]))

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.

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.

Written by Brad Chapman

January 31, 2009 at 3:11 pm

2 Responses

Subscribe to comments with RSS.

  1. Hi Brad

    Just came across your site and you might be interested in looking at what we are doing at Xyggy ( which building item-based search services. This page will tell you more about item-based search – There are also 5 demos to play around with to get the idea.

    Why have I written to you? We have long thought how useful item-based search would in the biotechnology industry but it is incredibly crowded with all sorts of people and companies doing all sorts of search things with data. Plus, there is BLAST. So, we’ve been reticent about investing time and resources into this sector. However, you may know different and would be good to communicate on the subject.

    We are a startup located in Mountain View, California and Cambridge, England.

    Best …


    dinesh vadhia

    February 4, 2009 at 10:18 pm

    • Dinesh;
      Thanks for the note. The Xyggy demos look great and your search approach makes a lot of sense. There are certainly a lot of ways to search for items related to a protein of interest, but not really one default way that dominates the space. You can look based on: sequence similarity (with BLAST, from your example), literature, proteins it interacts with, other genes it is expressed with, proteins it is evolutionarily related to; that’s off the top of my head. You also have the problem, like in movies or music, of having two different types of searches. One is when you have just encountered a protein, and would like to know all the obvious things you should look at that an expert would already know. The second is when you are already familiar with a protein and want to dig deeper and pull out novel similarities for research.

      This blog helps me play around with some of those ideas while I am thinking about my own work. In terms of Xyggy, I do think that scientists would be open to new search approaches that help organize and present related proteins or papers based on many criteria. There isn’t an obvious path to a great search, but rather a lot of room for innovation.


      Brad Chapman

      February 5, 2009 at 7:06 am

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: