GeneProf Web Services

If you're running a Unix-like operating system, you're probably familiar with the concept of 'piping' the output of one commandline program into another (cp. this Wikipedia article on Unix pipelines). You can use wget to retrieve GeneProf data as a stream like this:

Here's an example of how to retrieve a FASTQ sequence file, which is then filtered for sequences containing the nucleotides ACTG (in order) and written to a file called filterd-seqs.fq:

You can combine the outputs of several web services to achieve more advanced results. For instance, let's first look up the GeneProf ID of the mouse gene with the Ensembl ID ENSMUSG00000024406 (that's the transcription factor (TF) Pou5f1), and then look for all genes that are putatively regulated by this TF across all ChIP-seq datasets available in the GeneProf database and then check for genes from the "Sox" family:

Hint: Of course, you're not limited at all to using standard Unix commands! There's nothing holding you back from making use of great bioinformatics commandline tools such as, for example, Biopieces, FASTX-Toolkit or BEDTools.

Another use of GeneProf's web services is for loading data into R. Many types of data can be exported directly as a binary RData objects -- which can be easily loaded into an existing R session like this:

Let's put this into a function (called loadGeneProfData) for even easier use:

We can use this to write a function that combines two web service calls to look up the GeneProf Gene ID for a gene symbol and then to retrieve gene expression data for this gene:

Let's combine these methods to retrieve some expression data and generate a few plots:

In the example above, we first get the GeneProf ID's of the mouse genes with the common symbols Sox2 and Pou5f1 (with the getGeneProfID function) and then use this ID to query expression data (in RPKM format, by default) via the Get Gene Expression Values for a Gene servlet described above. We could now do anything we like with these values, for example, we haven chosen here to plot two histograms of the RPKM values for each gene and a scatterplot comparing them.

In addition to the expression values, the same web service makes it possible to retrieve additional annotation data, e.g. for the cell type of each observation. Let's get the additional annotations for one of the genes (the annotations for the other would be the same, so no need to get them twice) and then use these annotations to plot an annotated scatter plot for a selection of cell types (the dots for each cell type will have a different colour in this plot):

The complete R script for this example is available here.

GeneProf stores quite a lot of genomic data, e.g. from alignments of short reads or enriched regions, such as peaks detected in ChIP-seq experiments. It is possible to load this data directly into external applications that understand the commonly used WIG and/or BED file formats such as genome browsers, like the UCSC genome browser, the Integrative Genomics Viewer, and others. The key web services required for this task are Genomic Data as BED Files (BED) and Genomic Data as WIG Files (WIG).

We'll try it with the UCSC Genome Browser first. On the genome browser homepage, we need to click the Add/Manage Custom Tracks button to start up the interface for uploading custom tracks. There's a big text box which allows users to enter the URLs under which custom tracks can be found and this is where the URLs for the relevant GeneProf web service calls are meant to go. Let's load in alignments of ChIP-seq datasets (Pol2 ChiP-seq and Input DNA control from this experiment based on data from Sultan et al. (2008)) as WIG tracks to see the coverage of reads in these datasets across the genome as well as the coordinates of ChIP-seq peaks that were identified in the analysis of these tracks. Let's focus on chromosome 3 only, for faster loading times. Copy and paste the URLs below into the UCSC's text box:

Just click the Submit button and that's it. Your data should be imported into the UCSC browser rather quickly. You can then continue using the genome browser as you are used to, e.g. have a look at the region chr3:37,253,759-37,439,292 with a nice peak at the TSS of GOLGA4.

Fancy another browser? The same URL patterns should work with IGV (see screenshots below) and probably pretty much every other genome browser. Note that, unfortunately, none of the browsers we've tried seem to work very well with URLs that contain query parameters (that's the argument behind the question mark (?), so further configuration of the tracks will not be possible unless you download the tracks first and upload them afterwards.

Galaxy is an open-source, web-based platform for genomic data analysis that has attracted a vibrant user-base over the recent years. GeneProf data, especially genomic data, can easily be loaded into Galaxy for further analysis.

As a simple example, we'll load two sets of ChIP-seq peaks for the factors Stat3 and Klf4 into Galaxy using a call to the Genomic Data as BED Files (BED) service. We'll then use Galaxy's built-in tools to intersect those regions to discover regions of potential interaction between the two factors.

On the Galaxy homepage, we start by opening the Upload File tool (found under Get Data). We'll upload two sets of ChIP-seq peaks, both contained in the GeneProf dataset gpDS_11_12_125_2 using web service calls to the URLs:

Simply copy & paste both URLs into the text box marked "URL / Text". The files we're retrieving are from the mouse genome and come in BED format, so let's select bed where it says "File Format" and Mouse July 2007 (NCBI37/mm9) where it says "Genome". In the end, click the Execute button to start the upload.

This should only take a moment (depending on how busy the servers are just now) and then two new datasets should be available in the Galaxy workspace for you to play with. For instance, we can now use the Intersect tool (under Operate on Genomic Intervals) to overlap both sets of peaks. In the configuration of this tool, just pick any one of the two datasets under "of:" and the other one under "that intersect:" and the execute the tool. After a few minutes you should be left with a new dataset of all overlapping peaks.

That's just one example, but in the same way you can upload lots of other genomic or sequence data from GeneProf directly into Galaxy enabling you to exploit both tools to their best potential!

General-purpose pipeline execution engines have many prospective uses and have enjoyed immense popularity in many scientific fields. One of the most well-known solutions in the life sciences is probably Taverna. Using GeneProf's web services it is possible to wire in GeneProf data with other tools into arbitrary pipelines, let's look at a simple example..

The key to using GeneProf data in Taverna is the REST Service template offered by Taverna. You can use this template to connect to the GeneProf web services listed above, filling in wild cards ({mywildcard}) with user inputs as required.

In the example workflow shown above, a concatenation of two web service calls is employed to retrieve gene expression data for an arbitrary gene. We first use the web service Get the GeneProf ID of a Gene to get the internal GeneProf gene ID for a gene based on an external identifier or name. We configure the Taverna module with the URL pattern of the web service leaving space for three wildcards (step 1):

Please also change the 'accept' type to text/plain, because this is the type of data returned by that service. The wild cards are to be filled with user-defined inputs, so we next add three user input modules and connect them up to the web service box (step 2).

The GeneProf web service returns a plain text file with one GeneProf ID per line. We need to parse these into a list of Strings in order to use them futher on. The local Taverna service Split string into string list by regular expression will do the job (step 3): Connect the responseBody output of the REST service module to the string input of the regular expression splitter and define string constant to feed into the regex input. The string constant should read \n, i.e. the regular expression is to split the input along all "new line" characters.

Finally, we can add another REST service module to address the Get Gene Expression Values for a Gene service. Again, we leave space for two wildcards:

The value of the {id} wild card comes from the parsed output of the regular expression string splitter module and the {id} wild card receives the input from the same user input parameter as the first module (step 4). The 'accept' type of this module should be application/xml since we're requesting XML data this time.

In the end, we define the workflow outputs of the new Taverna pipeline and we can now rename the modules used and annotate inputs and outputs to make the workflow a little nicer.

You can download the complete example Taverna workflow here.

Other web sites across the globe can dynamically retrieve live GeneProf data and include it in their web pages using AJAX and JSONP cross-domain request, for example via jQuery (also see this StackOverflow post).

As an example, we'll implement a simple HTML page that displays a search form for genes. Using jQuery, the form will request matching genes from the GeneProf web services (Search Genes) and list them on the page along with their GeneProf IDs, a description and a plot for its gene expression values in a selection of cell types.

Click here to show the complete example page in an <iframe>!

Instead of listing the entire source code for the page, let's just look at its most important component, the JSONP request (the complete source code is available here):

As you'll notice the url attribute points to the URL of the GeneProf web service and dataType is set to jsonp (not json!). There are two important callback functions: success and error, the implementation of which defines in which way the response data is being dealt with.

In our case, we parse out a part of the response (gene names, etc.) and create some HTML code accordingly, and then fire up another AJAX call to another web service to retrieve some expression data for each gene. We then use the d3.js library to render an SVG bar plot showing the average expression in all available observations for a selection of cell types (this requires a compatbile browser -- pretty much anything but IE8 or earlier should do!).

Check out the complete source code of the example page here and have a good read through the jQuery documentation and you'll find you can easily integrate GeneProf data into your own web site, too!

GeneProf web services can be easily integrated with many popular scripting languages, e.g. Python or Perl. Most languages provide libraries for parsing XML files, so that's usually the kind of output you'd want to retrieve from the GeneProf web services.

Let's look at a simple example in Perl: We'll retrieve a bunch of gene expression values (RPKM) for an arbitrary gene and calculate the average gene expression per cell type. To make this as easy as possible we'll be using the REST::Client library.

Complete example script: Perl example script, Required Perl modules: XML::LibXML and REST::Client (along with their dependencies)

So, let's write some Perl scripts..

First, we want to set up a new REST web service client and configure it with the root URL of the GeneProf web services. We then use this client to perform the search and retrieve the data in XML format from the Get Gene Expression Values for a Gene service (we'll arbitrarily pick the human gene with the ID #2981 here as an example, but you may, of course, substitute other values as you please).

REST::Client provides a mechanism to automatically interpret the response as XML, so we can easily use XPath to get a list of all returned observations, which we can iterate one by one:

In order to calculate the average expression per cell type, we need to get (a) the expression value for each observation and (b) the cell type from the sample annotation of the observation. The expression value can be retrieved from the RPKM child node of the observation:

In order to get the Cell_Type sample annotation, we first get the sample child node and then retrieve the value from there:

Some observations might not have a cell type annotated, so we need to fill in missing values:

Assuming we had already defined as hash map groupAverages, we could now add up the totals and observation count for the current cell type (which we use later on to calculate the average):

After the main loop and the data retrieval work is done, we just need to print out the results:

You can download the complete Perl example script here.

In this advanced example, we'll be using the same techniques as in the simple example above to achieve a more complex task: We'll search GeneProf for all experiments with linked publications in Cell Stem Cell, retrieve metadata about these experiments, from which we'll find out what the main outputs of the analyses were and then we take all the mouse-specific gene-centric data from these datasets and merge them all together into one big table, which will be printed out as a tab-separated text file.

Complete example script: Perl example script, Required Perl modules: XML::LibXML and REST::Client (along with their dependencies)

Step 1: Search for experiments

So in the first step, we need to set up a new REST client and configure it with the root URL of the GeneProf web services. We then use this client to perform the search and retrieve the data in XML format from the Search Experiments service.

REST::Client provides a mechanism to automatically interpret the response as XML, so we can easily use XPath to get a list of all returned experiments, which we can iterate one by one:

Step 2: Find the main outputs of each experiment

Next, we need to get metadata about each experiment using the Metadata about a GeneProf Experiment service. We'll use the with-outputs=true parameter to include a list of all output datasets in the metadata.

To retrieve the metadata we just use another call with the REST client:

Using XPath we can now easily parse out all the outputs -- we can even directly filter out datasets with data type FEATURES based on the reference dataset pub_mm_ens58_ncbim37 (that's the mouse reference dataset!):

Step 3: Retrieve and merge the actual data

The last web service we'll be using is Data as XML. In the previous steps, we've only been retrieving metadata about experiments and datasets, whereas this service retrieves actual (sort of tabular) data, in this case gene expression values and regulatory data specific to genes.

Again, we'll use the previously set up REST client to retrieve the data:

The response will contain lots of data for genes organised in rows and columns. Add the beginning, the XML file contains additional meta information about the individual columns. We'll parse both separately:

In order to merge all data together, we'll build up to large hashes-of-hases, called allData and allColumns. Since all datasets are based on the same reference the internal GeneProf IDs can be used as a merge criterion for the data rows. At the same time, the column identifiers should be unique unless the columns are shared between datasets, e.g. C_NAME defines gene names and is identical for all datasets, so it can be safely overridden:

Step 4: Write out the results

Add this point, the two variables allData and allColumns should point to two large hashes-of-hashes containing all the data with retrieved merged into one, so all that's left to do is write out the data!

We'll write the data as tab-delimited plain text, so let's first write a header line:

And now the data:

Yo, that's it. You can download the complete Perl example script here.

Lastly, let's look at an example using an object-oriented programming language, Java. We'll write a little program that looks up genes by name and then calculates the average gene expression value (RPKM) for each of these genes per cell type -- based on all the RNA-seq data in GeneProf. In this example, we'll only use basic Java classes, but in practice you might want to give any of the REST client libraries out there a try.

Complete example code: Java Code, Compiled Example Program:: GeneProfWebServicesJavaClient

We'll not go through the entire code of the program here (the complete source code is available here), but just focus on the most important parts. So let's start by defining a new class that will act as our web service client:

We can now define some generic methods to retrieve data as XML or plain text from the GeneProf web services (the code's abbreviated a bit -- of course, you should always make sure to close streams and connections properly!):

Using those generic methods, we can define others that perform more specific operations using the GeneProf web services, like this function which looks up internal GeneProf gene IDs corresponding to a given gene name / symbol:

The remainder of the source code is a bit too bulky to post here, but please just refer to the (commented) source code available here.

Once you compile the code into a JAR file (or download it here), you can run the program like this: