Let me back up a bit. We're into grid computing. We have a boat-load of data from various plant genomes that we like to go to town on. On our campus we happen to have a 2,000 core compute farm (known as BlueHelix) and we use Sun Grid Engine (SGE) for job management. One nice feature of SGE is the ability to submit array jobs. An array job is a convenient SGE semantic for grouping tasks that are very similar. You submit an array job once and SGE spawns individual compute tasks that inherit the environment and run parameters of the job. Tasks are identified by their index number. A simple array job submission looks like this:

$ qsub -N job_name -t 1-10 ./script.sh

The job is named job_name and is composed of 10 tasks. While the job is queued or running, it's very easy to track and manipulate that job using its chosen name. For example, to find the status of the array job, run:

$ qstat -j job_name

SGE also lets you specify command-line parameters on special lines within a submission script that are only interpreted at job runtime. So you can submit an array job by running

$ qsub ./script.sh

where the first few lines of ./script.sh might look like this:

#!/bin/sh
 
#$ -N job_name
#$ -t 1-10

Each spawned task is assigned a unique environment variable, $GE_TASK_ID, that identifies its index in the array. It is this variable that distinguishes one task from the next (it can be used programmatically inside ./script.sh). This presents a slight complication because not all jobs can be naturally indexed by number. How can an array job be used to process a batch of arbitrarily-named files, for example? It is the submitter's responsibility to dereference the task ID to a meaningful value that can be processed by a given task.

Today I had a specific task: Dump all the maize chromosomes, soft-masking repeats, from an Ensembl database. Maize has 10 chromosomes, and an additional 11th chromosome we call UNKNOWN, where a small number of unanchored BAC clones are arbitrarily concatenated. I've written a handy script in the past that uses the EnsEMBL API to dump any specified region, with the option to soft- or hard-mask repeats. Dumping the ~2 Gb maize genome serially takes a good number of hours. With repeat-masking (the genome is more than 80% repetitive), it also requires a lot of memory. Naturally, we can use parallelization to speed up the dump. The simplest approach (not necessarily the fastest, or the best use of a 2,000 node cluster) is to divide the genome by chromosome. The job array would submit 11 tasks. Of course, the tasks aren't normalized, as the chromosome lengths vary rather wildly (The largest chromosome, 1, is more than twice the size of the smallest real chromosome, 10).

Within the submission script, we can use the $GE_TASK_ID as is to pass into the dump script the chromosome being dumped by the task. This is all good for 10 out of 11 chromosomes. We just need to handle the special case of the unanchored chromosome. The submission script looks like this:

#!/bin/sh
 
#$ -N job_name
#$ -t 1-11
 
chr=$GE_TASK_ID
[ $chr -eq 11 ] && chr="UNKNOWN"
 
./dump.pl --region $chr

Not very generalizable but illustrates how the task ID can be used to dereference meaningful values (a good general approach is to have task IDs be line numbers in a control file where each line contains the meaningful value).

A problem I ran into is that, using the default submission parameters, 5 out of 11 chromosome tasks failed because they ran out of memory (ironically, the largest chromosome was dumped successfully, probably because it ran on a high-memory node). I wanted to resubmit just those failed tasks without rerunning the 6 completed tasks. But as I said in the beginning, I didn't find a straightforward way to resubmit those 5 jobs. The -t parameter does not accept a delimited list of tasks, unfortunately (it does let you do other nifty things, such as specifying one range, e.g., 3-8).

The solution I came up with took advantage of another cool SGE feature that allows you to pass in a script through standard input. This lets you write a wrapper script that does some preprocessing before printing (typically via a heredoc) a shell script to standard output. Submitting then looks like this:

$ ./submit.sh | qsub

For this illustration, I had the following chromosomes fail to dump: 2, 4, 6, 7, and 8. To generally specify specific tasks to run, I used a bash array variable that held the failed chromosomes, then used the array length to pass in the array job parameters. The final script looks something like this:

#!/bin/sh
 
chrs=(2 4 6 7 8)
 
cat <<SCRIPT
#!/bin/sh
 
#$ -N dump_chromosomes
#$ -t 1-${#chrs[*]}
 
chrs=(`echo "${chrs[*]}"`)
 
chr=\${chrs[\$SGE_TASK_ID-1]}
 
./dump.pl --region \$chr
SCRIPT

This isn’t that pretty, but it gets the job done, no pun intended. Some sigils ($) had to be escaped to avoid interpolation. That’s because we have to distinguish between variables that are interpreted as the wrapper is executed and variable names that remain intact so they can be interpreted later when the task is submitted.

The first declaration creates an array of chromosomes that need to be re-dumped (this can easily have been the command-line argument list, which would allow you to run the job as ./submit.sh 2 4 6 7 8). The script then prints another script where the job array parameter is reconstructed using the length of the array (${#chrs[*]}). A slight inconvenience is that we had to also reconstruct the same array within the generated script so that the task can use it to dereference its corresponding chromosome.

This approach no longer treats the $GE_TASK_ID as incidental chromosome names, but rather dereferences meaningful values out of some array ($chrs, in this example). This is one illustration of how job array task IDs can be dereference in a general way. Submitting jobs by a generated script also allows you, as a sanity check, to see the script that would be submitted (by not piping it into qsub).

So what have we learned today? Besides the fact that the maize genome is huge and mostly repetitive, there are various sleek ways to submit jobs to a compute farm using SGE. The SGE qsub command is very flexible and provides powerful ways of submitting — and resubmitting — parallel jobs.

The post essentially argues that bioinformaticians are — to mix metaphors — too quick to jump on the cloud computing band wagon. It often makes more sense — and cents — to analyze genomic data locally (or privately). His cost-benefit analysis is spot on, but applies to the state of the art. The state of bioinformatics is, however, changing very rapidly and the scale (pun very much intended) is likely to tip.

Cloud and High-Performance Computing

Cloud Computing has become a bit of an overused catch-phrase, especially in bioinformatics. In the sense that it is a service, cloud computing is extremely essential for the life sciences. It allows researchers to integrate data sets and infer latent relationships among them. Reference data sets like the human genome are made ubiquitously available on “the cloud,” and scientists are then able to query those data sets seamlessly against newly-procured experimental data.

The term “cloud computing” has been liberally applied to high-performance computing (HPC). Biological data is complex and requires fairly sophisticated algorithms to interrogate it in a reasonable timeframe. Fortunately, many bioinformatics use-cases, such as read-mapping, are embarrassingly parallel and hence scalable. The problem can be rapidly chopped into fragments, each analyzed in parallel, and then combined in the end. Other problems whose solution rests on complex interdependent data structures, such as de novo genome assembly, cannot be as easily decomposed.

Cloud computing has also become synonymous with the MapReduce framework and Hadoop, its primary open-source implementation. This is somewhat unfortunate because Hadoop provides other benefits that make it suitable for high-performance distributed computing. CloudBurst and Crossbow, both developed in the same University of Maryland group, are the first published bioinformatics applications to use Hadoop. They smartly and effectively demonstrate the utility of MapReduce and cloud computing in bioinformatics. They leverage important Hadoop features such as task management, configurability, rapid prototyping, and fault-tolerance. But they have certain limitations. After reading the papers, I remain unconvinced that they are more efficient or scale better than an embarrassingly-parallel pipeline. In Crossbow, read-mapping is relegated to the mapping step while SNP-calling is relegated to the reduce step, but I’m not sure that this is the best fit of the MapReduce paradigm to the problem. Also, while the correctness of the pipelines is well-demonstrated for the data sets used, I doubt that their component applications (e.g., Bowtie, the unfortunately-named SOAPsnp) are universally applicable. Making the read-mapping and SNP-calling configurable would be a step in the right direction.

The point is that CloudBurst and Crossbow are heralded as HPC applications (even though their performance is similar to other existing pipelines) when in fact their primary appeal is to cloud computing (in the strict definition), namely in that they are Software-as-a-Service (SaaS) or Platform-as-a-Service (PaaS). In the HPC sense, Dooling is correct to point out that analyzing data locally makes more economic sense in the near term. In the cloud computing sense, it is hard to measure the benefits, since many of them are intangible (for example, distributed access to biological data and analytics).

IT Resources

According to Dooling's numbers, "unless you are just sequencing a few genomes, you are probably better off purchasing a (possibly single node) cluster." What his numbers don't take into account is the overhead of running a (possibly single node) cluster. While the fixed cost of purchasing computer equipment might be manageable, especially compared to chemical reagents, the operational costs of running a data center are substantial. Computer equipment needs to be continually serviced, be it for software, security, or kernel patches, or for unscheduled maintenance. In addition, energy costs for running a data center are high and expected to increase in the near future.

Data Storage and Localization

Another major consideration for analyzing data, be it locally or on the cloud, is its shuttling and storage costs. An Illumina sequencing run can easily generate hundreds of gigabytes of compressed data. That data comprises both sequenced reads and quality information. The Crossbow paper reported that transferring 103 Gb of sequence data to the Amazon EC2 cluster for processing cost $28, which is about a third of the cost to process it ($85). It also took 1 hour and 15 minutes for the transfer. That effectively increases the wall time for a Crossbow job, especially when compared to a locally-run instance (and that could be an optimistic estimate, since many of the smaller research groups — poised to be beneficiaries of cloud computing services like Crossbow — don't have access to local servers with significant upload bandwidth).

On the other hand, the cost of marshaling data to and from the Amazon EC2 cluster as well as of long-term storage on the Amazon Simple Storage Server (S3) already reflects the total data center costs that Amazon has to bear. Long-term storage of data in a local data center is expensive, especially when thinking about backup requirements and data availability. Data storage costs should not be separate from computer equipment purchases or the aforementioned IT resource considerations.

Sequencing Scale-up

Next-gen — or rather now-gen — sequencing platforms like 454, Illumina, ABI SOLiD, and Helicos, are pushing the bioinformatics and biodata envelope. More and more emphasis is placed on data wrangling and management. And these sequencers are becoming more affordable (e.g., improved protocols require less reagents) as the data they produce is getting larger (e.g., longer high-quality reads). But more than improved technology, the niche for these sequencers is growing rapidly. It didn't take long to realize that these sequencers can be used for cheaper genotyping over microarray-based approaches. But now with RNA-seq and ChIP-seq applications, it's becoming apparent that sequencing is not just for genome projects anymore.

In the meantime, the next generation (Gen. 3) of sequencing platforms, like Pacific Biosciences' SMRT technology or the IBM nanopore sequencer are already on the horizon, promising higher and faster throughput. They are likely to bring yet another disruptive shift in DNA interrogation.

As Dooling aptly puts it, "the more sequencing data people get, the more they want." As sequencing costs get lower and yields get higher, their niche in the scientific marketplace will grow. More and more research initiatives will rely on sequencing, at the very least for post-hoc validation. And as the data grows at an exponential rate, small research groups will discover when it's too late that they aren't equipped to handle the deluge. Cloud computing, and especially cloud storage, may in the very near future provide an option for long-term persistent storage of raw biological data. There will be some incurred costs (data upload and storage fees), but those costs will be miniscule compared to the local IT infrastructure that will be required to manage the data.

As things stand, the role of cloud computing in bioinformatics is fairly limited and undefined. There isn’t a complete suite of analytic utilities available for non-developers to manage and analyze their data on the cloud. As proposed by Dooling, it might make sense for a researcher to process all their data on a cheap commodity computer they can order on NewEgg.com.

But as things stand, change is the only constant. Bioinformatics is very much in flux, as the life sciences are quickly becoming a data science. In the future, small research projects will not only rely on high-throughput sequencing, but on bringing multiple systems-wide data sets to bear. Researchers will not be able to hire or train a developer as an afterthought. IT infrastructure is an absolute necessity at the inception of the project. The question is, where does it make most sense to manage IT resources, locally or remotely? The answer to that question will depend on an institution’s existing infrastructure. Large IT departments with a mature data center might be able to handle the capacity. But for small groups it make more sense to outsource the IT infrastructure to a growing, scalable data center known as the cloud.

Applications like Crossbow illustrate that cloud computing can play a very important role in the life sciences. Over time robust cloud-based analytics and visualization tools will be made available not only to programmers but directly to scientists. Ubiquitous access to both data and tools will promote both good science and collaboration. That is a primary promise of cloud computing. A secondary promise is that high-throughput analytics will allow users to interrogate the data more rapidly. MapReduce frameworks obviate the need to develop robust pipelines and let users concentrate on specific queries. As more MapReducible bioinformatics tools such as Crossbow are brought online, the scientific paradigm is likely to shift from hypothesis-driven to hypothesis-generating.

There is a limit in the mysql client and server for the maximum size for network transfer. This is one of the main reasons for chunking. We are moving away from running analyses on the seq level coord system much anymore. Having features on toplevel coords makes access much faster, so generating everything on that coord system makes sense. Doing this makes the underlying chunking fairly irrelevant. The exception would be if you wanted to map to some third coord system through the chunks, but the reason chunks are usually made is because there is no real assembly information (an agp), just the final assembled sequences.

The example that I am working on right now is the Arabidopsis lyrata genome, which has a maximum scaffold length of 33 Mbase. So what is the ideal chunk length (if any) in this case?

Having looked at the Ensembl schema (the dna.sequence field is LONGTEXT), the MySQL manual, and the configuration of the public Ensembl database, it seems that the maximum realistic length of a sequence is 16 Mbase. Although I could increase my mysqld’s max_allowed_packet to encompass my 33 Mbases I feel this would be unwise (and incompatible with most peoples clients), so I’m stuch with chunking for this genome.

One problem with chunking is that it introduces artificial sequence breaks in downstream analyses that operate on the seqlevel coord_system (input_type_id=CONTIG, in ensembl-pipeline parlance). I cannot seem to find a list of which analyses use which input types, but I assume that algorithms which have been designed to run on sequences such as clones will use seqlevel; examples include RepeatMasker (max sequence length 4 Mbase), GenScan (max 5 Mbase) and so on.

I have always blindly used a chunk size of 100 kbase in the past, but I’m now thinking that 1 Mbase may be more appropriate. We’ll have to wait and see…

In this release, we bring you 26625 gene annotations from TAIR, Gramene, SGN and MaizeGDB, 8558 QTL annotations from Gramene, 9832 germplasm associations from SGN, MaizeGDB and NASC.

The following number of annotations were added for the first time: 16016 gene annotations from TAIR, 3 gene annotations from Gramene, 2 gene annotations from SGN and 3928 germplasm annotations from MaizeGDB.

For genes curated by TAIR and SGN, you may also find links to their Gene Ontology (GO) pages through PO browser.

The new ontology files and the database dump are available for download.

To submit plant ontology term requests, we encourage researchers to use SourceForge PO tracker.

The Plant Ontology Consortium
web: http:www.plantontology.org
e-mail: po-dev at plantontology.org

The project is funded by National Science Foundation, USA, (Grant No. DBI-0703908)

Now, most Perl programmers know that it’s very useful to pass in a hash (or hashref!) of parameters for functions that are very very long. It’s useful and keeps things tidy.

Compare, which would you rather read:

function(0, 'true', 'apple', undef, undef, undef, undef, undef, 'porcupine');

or

function(
    'reptile_count'  => 0,
    'in_forest'      => 'true',
    'preferred_fruit'=> 'apple',
    'alternative'    => undef,
    'nomenclature'   => undef,
    'predator'       => undef,
    'food_source'    => undef,
    'known_zoos'     => undef,
    'animal'         => 'porcupine',
);

Even better, all of those undefs are superfluous, since the key’s value is already undef. We don’t need to pass them in. Well, unless of course there’s a default value that we instead want to explicitly set to undef, but more on that in a second.

function(
    'reptile_count'   => 0,
    'in_forest'       => 'true',
    'preferred_fruit' => 'apple',
    'animal'          => 'porcupine',
);

And it’s much tidier. Just by looking at the function, you can probably guess that the animal you’re dealing with is a porcupine, he’s in the forest, he likes apples, and there are no reptiles about. Good luck figuring that out from the original version of the function with all the unnamed parameters. Very tidy, very neat, very compact.

But then, almost inevitably, I see functions that have huge blocks of code at the top of them to copy out the hash arguments and set defaults.

sub function {
    my %args = @_;
 
    my $reptile_count   = $args{'reptile_count'}         || 0;
    my $in_forest       = $args{'in_forest'}             || 'false';
    my $preferred_fruit = $args{'preferred_fruit'};
    my $alternative     = $args{'alternative'};
    my $nomenclature    = $args{'nomenclature'};
    my $predator        = $args{'predator'}              || 'humans';
    my $food_source     = $args{'food_source'};
    my $known_zoos      = $args{'known_zoos'};
    my $animal          = $args{'animal'};
    .
    .
    .
 
    do_something($animal);
}

And it just gets worse as more and more arguments are added. The block grows gigantically and it’s just a huge messy wart. Line after line after line of useless setup information.

Plus, there’s a subtle issue which sometimes pops up. Take a look at the example function, for the predator variable. If no value is passed in, then it uses the default of “humans”. But what if you actually wanted to set it to undef? Say the porcupine actually has no predators (does it?). This function won’t even allow it. Admittedly, that may be a good thing. Or it may not. It can be a subtle bug to track down.

Secondly, tell me quickly, at a glance, which variables have defaults and what are they? You have to scan the whole block to see them. Admittedly, proper columnation of your defaults on the right side will help, but still.

My preference is to just leave all of the passed in arguments in a hash and be done with it. No need to copy out into separate variables. Sure, you have a smidgen bit more typing later in your function, but you also have the added advantage of knowing exactly when you’re operating with an argument to your function as opposed to a variable you created yourself.

Compare my version:

sub function {
    my %args = @_;
    .
    .
    .
    do_something($args{'animal'});
}

Done! Much simpler. “But wait!”, you protest. “You’re no longer setting the defaults!” Easily done. Make them part of your hash definition:

sub function {
    my %args = (
        'reptile_count' => 0,
        'in_forest'     => 'false',
        'predator'      => 'humans',
        @_
    );
    .
    .
    .
    do_something($args{'animal'});
}

Ta da! You end up with much cleaner, simpler code. At a glance, you can see at the top of your function just which arguments have defaults set for them.

As an added bonus, you can now explicitly set “predator” to undef and that’s what your function will get. Here, you’re interpreting a list as a hash, so whichever key definition occurs last wins. Early on, it sets ‘predator’ to ‘humans’ (the default), but then inside of @_ we may have passed in ‘predator’ => undef, so it happily uses that instead.

Okay, but maybe you wanted to ignore keys specifically set to undef. We can still easily do this, we’re just going to define 2 separate hashes.

Finally:

sub function {
    my %defaults = (
        'reptile_count' => 0,
        'in_forest'     => 'false',
        'predator'      => 'humans',
    );
 
    my %args = @_;
 
    #set passed args to defaults, if no value passed in
    %args = map {$_, $args{$_} || $defaults{$_}} keys %args;
    .
    .
    .
    do_something($args{'animal'});
}

In my opinion, the resulting code is substantially tidier. You have a small block of defaults at the top, and the rest of the function specifically identifies which variables were arguments to begin with. No big, messy, superfluous wart to wade through.

Since then, I’ve been doing lots of things. In no particular order:

• Working fairly fruitlessly on bring our new blast.gramene.org machine online
• Answering user email, fixing CMap and Mart bugs
• Hacking a bit on SQL::Translator to help create Gramene schema ER diagrams
• Thinking a lot about how we could improve Gramene’s look and feel. This was the subject of a presentation to the group, one I’ll revise for the upcoming Gramene retreat.
• Reading a book on project managment and considering how I could improve. Apparently I should be having regular one-on-one meetings with people — revolutionary!
• Writing the 2008-9 Gramene grant summary. This has been a nice exercise for me as I go over the Gramene experimental plan and consider what we promised to do, what we’ve done, what we’ve got left.
• Polishing a brief paper on CMap for the journal Bioinformatics
• I created a test video of a Gramene tutorial I drafted and solicited opinions. It’s decent, but needs much improvement. I’ll tackle this again soon, breaking the 10+ minute intro into many smaller, 2-3 minute intros on each section. I will increase the font size and decrease the screen resolution to try to make the text more readable and will also use Omnidazzle, a free mouse highlighter, to make it easier to see my pointer. Sheldon asked me what software I used to record. I purchased (for $99) Screenflow, which I found relatively easy to use. I should also point out that, though I don’t really like the sound of my own voice (few do, I know), I believe I sound relatively decent because I used a decent microphone run through an inexpensive pre-amp

For several months I’ve really been thinking about ways to radically update CMap. My recent work on the CMap paper and answering some user email has made me want to do something even more, but then this week I really got in my head to try to incorporate Circos output into CMap. Looking into the Circos code, I think it will take some refactoring since all the logic is in a single Perl script. This will need to be move out into a module to make it easier to call it from other code. Since Shiran also wants to create Circos views from Ensembl, perhaps we could work together to create an API for Circos. We’ve each contacted the author, Martin Krzywinski at the BC Cancer Agency (Sheldon’s old employer), and we’re waiting to hear how we can move forward on contributing to that. It would be nice to see perhaps a Google Code project setup with SVN, Wiki, etc. I could be really excited to contribute to such a project. I will also add that I was very happy to see that Martin apparently found the SQLFairy on his own and thought that Circos would be another cool way to visualize schemas, too, so he created Schemaball. It says it’s limited to MySQL, so I wouldn’t mind turning that around and creating a Circos producer module for SQL::Translator so that this would be available for any database syntax. First, though, I need Circos to have an API, and that will require refactoring the code.

Anyway, here’s just one idea how Circos could radically improve CMap views. Take a typical, boring table from the correspondence matrix showing how two map sets compare to each other and create this view:

Pretty neat, I think. I could see making Circos an output option for the matrix and for map comparisons. At first I got excited thinking I might find a way to use Circos for all CMap’s graphical output since it is cool and allows multi-way comparisons as opposed to CMap’s limitation that maps can only be compared to those on their immediate right or left, but then I realized that the standard CMap view does have value such as when you want to align contigs to a reference map for order and orientation. Also, Shiran seems to think that performance with Circos might be a tad slow for dynamic queries, so it would be extra eye-candy that users would wait on. Lastly, there is the issue of making the image clickable with image maps. Circos doesn’t create anything like that, and, in fact, I seems like a very tricky think to do given that the regions are very curvy.

Along with creating new, exciting visual output, I’d love to completely rewrite the user interface for CMap. It’s time for it to have real Web 2.0 goodness. I’ve wondered if I should write it in Java as an applet or a separate download, but my colleague on the Gramene project, Terry, feels that this has been a hindrance to people accepting his GDPC tool on Gramene. While researching other ideas, I happened across the Chrome Experiments website and was thoroughly intrigued with what you can do nowadays with Javascript + HTML5 + Canvas + SVG (Jim pointed out all the technologies at work there). I’m thinking that’s gotta be the way to go. Whatever direction I choose, what I’d like is a much more dynamic comparative map viewer, from the form to the image itself. I would love to be able to update just a part of the image, e.g., if the user wanted to hide all the ESTs on a map, I could just redraw that one map rather than the whole image, but this might also entail redrawing the correspondence lines to neighboring maps, so I’m wondering if it would be possible to create the image in layers such as background, map, correspondences, legend, etc.

Another reason I would like to overhaul CMap is the performance issue. I’ve learned a lot since I originally designed the schema, and I’d like a do-over on several parts where I think I could make things much faster (particularly in how I store the feature aliases). Performance is not terrible anymore on the Gramene site since we considerably lowered the number of features we put on maps, but I would really love to see it be able to handle at least 1 million features per map — at present we try to keep it to under 50-100K.

At the heart of CMap is the idea that two things have a correspondence, they are related in some manner be it a sequence similarity or the same reagent being used in two experiments or whatnot. Modeling this efficiently has never been easy. Here is the core of the schema that shows how it’s currently modeled:

You can see in the feature/feature_correspondence/correspondence_evidence/correspondence_lookup area that things get a little circular, and that’s usually a bad thing. My first pass at this (oh, so many, many years ago, hacking on my little white iBook while riding the LIRR — my, the things I remember) had just the correspondence table, but then I realized that I had to have a table that showed the relationship going in both directions, hence the lookup table that says both “A->B” and “B->A” such that there are two lookup records for each correspondence. Ben came along and decided that, if the data was going to have to be denormalized to handle web queries, it might as well have mapping information in it in order to decrease the number of joins required, so he added the start/stop/type/map fields. Each correspondence can have one or more evidences.

To use these tables, I have some set of features in a region of a given reference map, and I want to find all the other maps that have features with correspondences optionally including or excluding some set of evidence criteria. I want the species and name of the map sets, the map names, and the number of correspondences. In the many times I’ve thought about how I could simplify or denormalize this schema further, I’ve never really come up with a better solution. I can’t see a way to pre-calculate anything. Every time a user crops a map or excludes some feature types, I have to go back and requery, joining millions of rows across several tables to figure out the choices I can present in the comparative maps menu. If nothing else, I’ve thought about binning the features such that I could perhaps at least get rid of the start/stop range queries, which can be a killer for performance. Another option I’ve considered is getting rid of the evidence table. I know this removes some flexibility (e.g., only show me correspondences based on BLAST hits with a score greater than “X”), but maybe it’s not worth it — or maybe I move it directly into the correspondence record itself, and each evidence gives you a new correspondence record. No, that would make it hard to count unique correspondences. Anyway, you see how complicated this can get.

I’ve done some looking about on modeling pairwise relationships but haven’t gotten far. I’d love to see how Facebook models “friends” — I know they have several million users, each of which can have hundreds of friends, so how do they efficiently store and query this? If anyone has any ideas for a radically improved version of how to show that a thing of one type in a particular location is related to something else, I’d love to hear about it. Something that is as simple as possible, but no simpler would be perfect.

That is all.

Gramene is going to be rolling out an autocomplete feature to offer suggestions to users sometime in the near future. You can sample the wonderfully suggestive goodness here until we go live.

The state of javascript in the world today continuously amazes me. I started writing javascript in the fairly early days. I go back to Netscape 2 with it and got seriously cranking developing things around 1997. I largely abandoned it in early 1999 when I got rather irate at the current implementations as I tried to write a clone of Pac-Man in it. I reached the point where Pac-Man would run around the maze and eat dots, but as I started to implement the ghosts I ran into a stumbling block – when the timers that controlled the ghosts were global objects, it ran fine. When they were attributes of the ghost object (so the ghost could move itself), then performance ground to a crawl. I got so annoyed at the arbitrariness of it all that I shelved javascript and didn’t look at it for a few years.

When I finally picked it up again, I was amazed to see that the current state of the world was actually all of the magical things that I’d been promised many years before. A consistent, easy to understand Document Object Model. The ability to independently update parts of the page. Dynamic pages that actually work. I was hooked again.

So. To Gramene. I had this bug kicking around for months but never took the time to delve into it. Conceptually it’s very easy – monitor the text boxes the user is typing in. Whenever they let go of the keyboard key, fire off an event to a handler. That’ll look at what’s typed in, and bundle it off to ship off to the server. The server will take that information and look it up…somehow…on the server side to get a list of suggestions. Those suggestions go back to the client and pop up in a little box below the text. Simple!

Mostly. At first I wanted to use Ken’s quick search routines and database, but that seemed like it might be slow for our purposes. And it didn’t have quite the information I wanted to see either. So him and I brainstormed on it and came up with a vocabulary list concept. We set up lists of vocabulary words for individual parts of the site – after all, if you’re typing in a text box on a Gene search page, you probably don’t want to see QTL information popping up. Ken was then nice enough to go populate these vocabulary lists for me. For those keeping score at home, we currently have 3,649,306 distinct vocabulary terms spread across 54 different vocabulary lists.

A few wrapper scripts and that part was done. The AJAX to hand off to the server side was straightforward. I’d written some AJAX routines for myself a few years back, and we integrated those into Gramene as part of the Help system. I needed to tweak a few options to make it operate properly in this case (for example – we only want to fire off one AJAX request at a time so as not to get suggestions back from the server out of sync), but it was basically done. The pop up box with the suggestions was just a DIV tag that was absolutely positioned slightly below the text box in question. We had a javascript function for that too. As a side note, finding the absolute position of a relative element on the page is still quite the nuisance, but I guess I can’t have everything I want yet.

I added in options for the user to disable suggestions, and then we were done! Time to blog it! This is when Ken fouled up the whole works by saying, “I’m having trouble autocompleting with more than one box in a form.” My response? “Crap.”

I had ever so carefully engineered and designed the thing and never really taken that into account. Because I wanted to keep the javascript simple to use, I set it up so that additional options for the script were set as hidden form elements. Which vocabulary list to use, for example. This worked fine when there was only one autocomplete box in the form. But once you had more than one, it blew up. After all, you might (as Ken did) want to define two separate search boxes which each use their own distinct vocabulary list. How would the server side script know which vocabulary item to use?

Further, how would it know which input box to look at?

I came up with 3 ways to solve the problem – the bad way, the good enough way (with some caveats), and the flexible way.

I skipped over the bad way, since it was, well, bad. I could tie the client side javascript and forms to the server side implementation. So that the functions the programmer calls in the browser do some particular parsing and analysis to submit the proper information to the server for the script to function. This is fairly easy to do, but introduced a tight coupling to the code.

Up until this point, the javascript and the server side were completely independent. It didn’t really care what the script was on the server side, it just handed off its parameters and got something back, which it stuck in a box. It could use the default implementation, or any other one. I felt that this flexibility was quite valuable – it allows us to easily deploy the code to our other probjects and it would allow us to easily open source it for anyone to use it, regardless of integration with Gramene. Maybe you have your own fast word lookup on your server. You shouldn’t need to change your code much (other than formatting the output) to use the results.

The good enough way was to extend the javascript to accept a lot of additional parameters. This way, Ken could hardwire in a bunch of values to allow the client side code to remain independent of the server. I set up a bunch of default values in the javascript that handle the usual case. You can even still pass in your values with hidden form elements (and are encouraged to!). But you can also send in all the options you need to your javascript. Voila! We have multiple forms!

I only consider that approach to be “good enough” because it has an irritating nuisance – what if the additional parameters we want to pass in are user configurable? Right now, it’s hardwired that one box searches the genes vocabulary list, and the other searches the QTL vocabulary list, for example. But it’s conceivable that we may want the user to choose which part of the site they’re searching, and then use that to populate the list. So no hardwiring values.

At that point, we need to analyze particular form elements and update arguments and hand around values. Or build new forms on the fly with carefully crafted definitions of what should be in them. Or any number of other nonsense like that. Or maybe something even more clever that I didn’t think of. Regardless, I opted to defer for another day until we actually needed the functionality.

As a result, after all that rambling about coupling and options and parsing, for most cases, it’s terribly simple to use.

<input
   type = "text"
   name = "search_for"
   id = "search_for"
   onkeyup = 'suggest(this.id)'
   onblur = 'hide_suggestions()'
   onfocus = 'show_suggestions(this.id)' />

And you’re off to the races. I always like it when software can be powerfully complicated on the backend, and completely hidden from the end developer so that most of the time, it’s trivially simple to use.

I will concentrate on getting a ‘minimal’ database populated, i.e. the sequences, genome assembly and genes (without additional annotation).

In my experience, the exact approach for each species differs, mainly due to the different formats and semantics of the source data. The following description will probably need tweaking for a particular species.

Basic input data;

• The sequences that are used in the assembly, preferably the unassembled contigs, but the assembled chromosome or supercontig sequences can be used if needed.
• A description of the assembly; how the contigs are mapped to form larger regions. This can be multiple levels, e.g. contigs-to-clone-to-supercontig-to-chromosome for instance. If you are working from large, assembled sequences you may need to ‘fake’ an assembly to chop the sequences into pieces that can be efficiently stored in the database.
• A file that describes the gene structures; the locations of exons and CDSes, and how these are grouped into transcripts and genes. It is also useful to have cDNA and peptide sequences for the genes for validation.

Preparation

You need to create an empty ensembl-schema database. The SQL to generate the schema is provide by a file in the ensembl project repository: table.sql

Loading the contig sequences

Only use this method if you have short sequences and an assembly file. If you have long sequences (e.g. whole chromosomes), see the Loading chromosome sequences and faking the assembly section.

Ensembl provides a script for loading sequences from a FASTA file into the database:
load_seq_region.pl. A typical invocation of this script is:

$ perl load_seq_region.pl -dbhost host -dbuser user -dbname my_db -dbpass **** \
  -coord_system_name contig -rank 4 -sequence_level -fasta_file sequence.fa

Loading the assembly

Once you have the sequences loaded, you need to load the identifiers of other regions that appear in the assembly. Again using the load_seq_region.pl script, these identifiers can be loaded either from a FASTA file (the sequences will be ignored);

$ perl load_seq_region.pl -dbhost host -dbuser user -dbname my_db -dbpass **** \
  -coord_system_name clone -rank 3 -fasta_file clone.fa

Or from an agp file (assembled regions);

$ perl load_seq_region.pl -dbhost host -dbuser user -dbname my_db -dbpass **** \
-coord_system_name chromosome -rank 1 -agp_file genome.agp

Once you have all the pieces of the assembly in the database (typically by running load_seq_region.pl a number of times), you need to load the assembly itself, i.e. the order/orientation/overlap between the assembly pieces. If you have agp file(s), ensembl provides a script to loat it: load_agp.pl. A typical invocation would be;

$ perl load_agp.pl -dbhost host -dbuser user -dbname my_db -dbpass **** \
  -assembled_name chromosome -assembled_version NCBI34 \
  -component_name contig -agp_file genome.agp

You will need to run this script for each level of assembly that you have.

Loading chromosome sequences and faking the assembly

The fastest way to get an assembly into Ensembl is to take the fully assembled chromosome sequences (or other large assembled component) for a species, an load this straight into Ensembl. As it is inefficient to work with huge sequences in single MySQL records, the approach is to chunk into smaller pieces and fake an assembly. A script that does this is available from Gramene: load_assembly_from_fasta.pl. Typical invocation would be:

$ perl load_assembly_from_fasta.pl \
  --ensembl_registry file=/path/to/ensembl.registry --species=Homo_sapiens \
  --coord_system=chromosome -assembly_version=v1 --chunk_size=100000

Loading the genes

Gene structures are usually contained in a gff-derived file format. The semantics of gff itself are a fairly lose, so it is better to use a more specific format, such as gtf or gff3. For a starting point, I direct you to a recent Gramene script used to load grape genes into the database; load_genes_from_grape_gff.pl. Typical invocation:

$ perl load_genes_from_grape_gff.pl \
  --ensembl_registry file=/path/to/ensembl.registry --species=grape \
  --logic_name=genes_irgsp /path/to/grape_genes.gff

Footnote

This is a very basic treatment on the subject of loading data into Ensembl. There are many optimisations and display issues that could have been addressed. As time allows I would like to extend these notes into a longer article on loading into Ensembl. If you have any suggestions, please post!

Even though Gramene is releasing only twice a year, we’re allowing interim updates of the site and important upgrades when they come along. Recently we recommissioned our old web server, filetta, as a public-access (and read-only) MySQL server with all of our Ensembl dbs, the markers db, and a few other dbs. I was able to use that to help a grad student at UT answer some questions about our genes db — how easy to tell him to connect directly to a running copy that he could query. (And, of course, I pointed him to the “mysqldump” so he could set it up locally for himself.)

Another release I made just this week was to reinstate the OMAP stacked maps in CMap. This was a challenge for me because Ben and Bonnie always did this in the past, so I had to follow their docs which ended up being pretty out-dated since so many things have changed since they both left. I hope that it will be easier this next release.

Since users often ask about our schemas, I finally finished up a “Build” action that I’d started long ago. If, from the Gramene root directory, you “./Build schemas” and “./Build schema_diagrams,” the Build.PL script will do a “mysqldump” of all our defined “” (in “gramene.conf”) into a “schemas” directory and then use “sqlt-graph” (from my SQL::Translator project) to create a schema diagram. It’s quite useful documentation to see how all the tables in the schemas are related.

One other nice improvement we’ve made recently to Gramene is the use of an actual blog for our news items rather than updating our index page all the time. I installed WordPress on bivouac and then run a cron job on the top of the hour to grab the latest three posts from the RSS feed, format them into HTML (taking care to respect Unicode characters after the problems with Pankaj’s post in Hindi), and put that into a file that gets included via SSI on the home page. It works very nicely, and I’m happy to see a much smoother system for publishing and archiving news items.

As we approach the data freeze in Gramene, I’m more determined than in previous builds to make it mean something. It seems data keeps flying well after the freeze, but we’re going to have to make this stick this time if there is to be time to build CMap and the Marts in a reasonable time. We’re only 13% through 156 items in our “data” build, so I’m wondering how much gets jettisoned or done by then. I’m slowly getting more comfortable with my status as a “decider” and enforcer, so I’ll just have to bring down the hammer when the time comes. Watch out!

So in a nutshell, aside from the technical aspects of things, what have I learned thus far?  First of all, programming is really fun!  As a biologist working with deep sequencing data, it is also really essential…at least a basic knowledge anyway.  Even knowing only what I know now would have made my research as a graduate student that much easier.  There are only so many lines you can populate in an excel spreadsheet before the computer crashes.  Especially with some of these Solexa datasets…such analyses would be virtually impossible without the proper codes.  Also, very importantly, I learned that programmers love Starbucks.  Nuff said…I’m in good company :)