Saaien Tist

Biovis/Visweek recap

2011-11-10T16:09:00.000+01:00

Finally time to write something about the biovis/visweek conference I attended about a week ago in Providence (RI)... And I must say: they'll see me again next year. (Hopefully @infosthetics will be able to join me then). Meanwhile, several blog posts are popping up discussing it (see here and here, for example).

This was the first time that biovis (aka the IEEE Symposium on Biological Data Visualization) was organized. It's similar to the 2-year old vizbi, but has an agenda that is more focused on research in visualization rather than application of visualization in the biological sciences. Really interesting talks, posters and people.

The biovis contest
This first installment of biovis included a data visualization contest, focusing on "specific biological problem domains, and based on realistic domain data and domain questions". The topic this year was on eQTL (expression quantitative trait loci) data, and I'm really happy that Ryo Sakai - now a PhD student in my lab - won the "biologists' favourite" award!! The biologist jury was impressed with the ease in which his visualizations of the eQTL data highlighted and confirmed prior knowledge, and how it suggested directions for further experiments. It was interesting to see that there was a huge variation in the submissions, going from just showing the raw data in an interesting way (which Sakai-san did) to advanced statistical and algorithmic munging of the input data and visualizing the end result (which the winner of the "dataviz professionals' favourite" award did). See how this relates to my previous post on humanizing bioinformatics?

Interesting talks - amazing (good & bad) talks
As this was the first time that I attended visweek, I was really looking forward to the high quality presentations/papers and posters. Overall, I got what I wanted. But there were some examples of papers and posters that I have major doubts about (taking into account that I have to be humble here in talking about people working in the field for far longer than I do).
One example that seemed pretty counterintuitive was a presentation by Basak Alper from Microsoft about a new set visualization technique that they baptized LineSets. The main issue that they want to solve is the visualization of intersections of >3 sets (up to 3 you'd just use Venn diagrams). Their approach is to connect the different elements from a set by a line; hence: linesets. However, I (and many others with me) felt that this approach has some very serious drawbacks. Most of all, it suggests that there is an implicit ordering of the elements, which there isn't. In the image below, for example, line sets were used to connect Italian restaurants (in orange) and Mexican restaurants (in purple). That's the only thing this visualization wants to do: tell me which of the restaurants are Italian and which are Mexican. But give this picture to 10 people, and every single one of them will think that the lines are actually paths or routes between these restaurants. Which they're not... The example below shows data that has specific positions on a map, but they demonstrate this approach on social networks as well.

LineSets

Another example comes from the biovis conference: TIALA or Time Series Alignment Analysis. Suppose you have the time-dependent gene expression values for a single gene, which you'd plot using a line plot. Now what would you do if you have that type of data for 100 genes? Would you put those plots into 3D? I know I wouldn't... And better still: would you then connect these plots so you end up with some sort of 3D-landscape? That's like connecting the tops of a barchart displaying categorical data with a line...

TIALA - Time Series Alignment Analysis

But of course there were plenty really good talks as well. Some of the talks I really enjoyed are those about HiTSEE (by Bertini et al) on the analysis of high-throughput screening experiments, EVEVis (Miller et al) on multi-scale visualization for dense evolutionary data, arc length-based aspect ratio selection (Talbot et al) which is an alternative to banking to 45 degrees, drawing road networks with focus regions (Haunert et al), and especially DICON which showed an amazing application of visual analysis of multidimensional clusters using healthcare data.

HiTSEE

EVEVis

Road networks with focus regions

DICON - interactive visual analysis of multidimensional clusters

Meeting interesting people
But of course this was very much about meeting interesting people as well. It was really nice to exchange ideas again with the biovis crowd (Nils Gehlenborg, Cydney, Tamara, Will Ray, ...), and I finally had the chance to have a chat with @filwd Enrico. All those discussions with Thorri from Icelandic DataMarket were both useful and fun (as was our day hanging out in town, chatting to the Occupy Providence woman (forgot her name, I'm afraid) and trying to find a good hat).
At the airport on my way back, as I was trying to find out how to get to Brussels (as our flights were cancelled due to the weather), a chap comes to me and introduces himself as someone from Belgium. From Leuven. From our very own faculty. So together with @infosthetics Andrew that now makes three of us :-)

Anyway: I'll definitely be back next year (have to play some more official role anyway) and already looking forward to it.

Humanizing Bioinformatics

2011-10-19T17:12:00.000+02:00

I was invited last week to give a talk at this year's meeting of the Graduate School Structure and Function of Biological Macromolecules, Bioinformatics and Modeling (SFMBBM). It ended up being a day with great talks, by some bright PhD students and postdocs. There were 2 keynotes (one by Prof Bert Poolman from Groningen (NL) and one by myself), and a panel discussion on what the future holds for people nearing the end of their PhDs.

My talk was titled "Humanizing Bioinformatics" and received quite well (at least some people still laughed at my jokes (if you can call them that); even at the end). I put the slides up on slideshare, but I thought I'd explain things here as well, because those slides will probably not convey the complete story.

Let's ruin the plot by mentioning it here: we need data visualization to counteract the alienation that's happening between bioinformaticians and bright data miners on the one hand, and the user/clinician/biologist on the other. We need to make bioinformatics human again.

Jim Gray from Microsoft wrote a very interesting book "The Fourth Paradigm - Data-intensive Scientific Discovery". Get it. Read it. He describes how the practice of research has changed over the centuries. In the First Paradigm, science was very much about describing things; the Second Paradigm (last couple of centuries) saw a more theoretical approach, with people like Keppler and Newton defining "laws" that described the universe around them. The last few decades saw the advent of computation in the research field, which allowed us to take a closer look at reality by simulating it (the Third Paradigm). But just recently - so Jim Gray says - we're moving into yet another fundamental way of doing science. We have moved into an age where there is just so much data generated that we don't know what to do with it. This Fourth Paradigm is that of data exploration. As I see it (but that's just one way of looking at it, and it doesn't want to say anything about what's "better" than what), this might be a definition for the difference between computational biology and bioinformatics: computational biology fits within the Third Paradigm, while bioinformatics fits in the Fourth.

Being able to automatically generate these huge amounts of data (e.g. in genome sequencing) does mean that biologists have to work with ever bigger datasets, using ever more advanced algorithms that use ever more complicated data structures. This is not about just some summary statistics anymore; it's support vector machine recursive feature elimination, manifold learning and adaptive cascade sharing trees and stuff. Result: biologist is at a loss. Remember Dr McCoy in Star Trek saying "Dammit Jim, I'm a doctor, not an electrician/cook/nuclear physicist" whenever the captain let him do stuff that is - well - not doctorly? (Great analogy found by Christophe Lambert). It's exactly the same for a clinician nowadays. In order to do a (his job: e.g. decide on a treatment plan for a cancer patient), he has to first do b (set up hardware that can handle the 100s of Gb of data) and c (devise some nifty data mining trickery to get his results). Neither of which he has the time or training for. "Dammit Jim, I'm a doctor, not a bioinformatician". Result: we're alienating the user. Data mining has become so complicated and advanced, that the clinician is at a complete loss. Heck, I'm working at a bioinformatics department and don't understand half of what they're talking about. So what can the clinician do? His only option is to trust some bioinformatician to come up with some results. But this is a blind trust: he has no way of assessing the results he gets back. This trust is even more blind than the one you give the guy who repairs your car.

As I see it, there are (at least) four issues.

What's the question?
Data generation used to be really geared towards proving or disproving a specific hypothesis. The researcher would have a question, formulate some hypothesis around it, and then generate data. Although that same data could already be used to answer other unanticipated questions as well, this really became an issue with easy, automated data generation; DNA sequencing being a prime example. You might ask yourself "does this or that gene have a mutation that lead to this disease?", but the data you generate (i.c. exome sequences) to answer this question can be used to answer hundreds of other questions as well. You just don't know what questions yet...
Statistical analysis and data mining are indispensable for (dis)proving hypothesis, but what if we don't know the hypothesis? As many others in the field, I believe that data visualization can give us some clues at what to investigate further.

Let's for example look at this example hive plot by Martin Krzewinski (for what B means: see the explanation at the hive plot website). Suppose you're given a list of genes in E.coli (or a list of functions in the linux operating system) and the network between those genes (or functions). Using clever visualization, we can define some interesting questions that we can look into using statistics or data mining. For example: why do we see so many workhorse genes in E.coli? Does this reflect reality, and what would that mean? Or does it mean that our input network is biased? What is so special about that very small number of workhorse functions in linux that have that high connectivity? These are questions that we need to be presented to us.

What parameters should I use?
Second issue: the outcome from most data mining/filtering algorithms depend tremendously on the right parameters. But it can be very difficult to actually find out what those parameters should be. Does there actually exist a "right" set of parameters for this or that algorithm? Also, tweaking some arguments just a little bit can have vast effects on the results, while you can change other parameters as much as you want, but it won't affect the outcome whatsoever.

Turnbull et al. Nature Genetics 2010

Can I trust this output?
Issue number 3: if I am a clinician/biologist and a bioinformatician hands me some results, how do I know if I can trust those results? Heck, being a bioinformatician myself and writing some program to filter putative SNPs, how do I know that my results are correct? Suppose there are 3 filters that I can apply consecutively, with different combinations of settings.

Looking at exome data, we main information that we can use for assessing the results of SNP filtering are the fact that you should end up with 20k-25k SNPs, and a transition/transversion ratio of 2.1 (if I remember correctly). But there's many different combinations of filters that can give these summary statistics. The state of the art (believe it or not) is to just run many different algorithms and filters independently, and then take the intersection of the results...

I can't wrap my head around this...
And finally, there's the issue of too much information. Not just the sheer amount, but of different data sources. It's actually not really too much information per se, but too much to keep into one head. Someone trying to decide on a treatment plan for a cancer patient, for example, will have to combine data from heterogeneous datasets, multiple abstraction levels and multiple sources. He'll have to look into patient and clinical data, family/population data, MR/CT/Xray scans, tissue samples, gene expression data and pathways. That's just too much. His cognitive capacities are fully engaged in trying to integrate all that information, rather than in answering the initial question.

Visualization... part of the solution
I'm not saying anything new here when I suggest that data visualization might be part of the solution to these problems. As current technologies and analysis methods have alienated the end-user from his own results, visualization can reach over and cross this gap. The rest of the presentation is basically about some basic principles in data visualization, which I'll not go further into here.

All in all, I think the presentation went quite well:

Visualize This (by Nathan Yau) arrived...

2011-09-01T10:27:00.003+02:00

Last Friday I received my long-anticipated copy of "Visualize This" by Nathan Yau. On its website it is described as a "practical guide on visualization and how to approach real-world data". You can guess what my weekend looked like :-)

Overall, I believe this book is a very good choice for people interested in getting started in data visualization. Not only does it provide the context in which to create visualizations (chapters 1, 2 and 9), it also handles different tools for creating them: R, protovis, flash.... Apart from chapter 3 that is dedicated entirely to that topic, different examples in the book were created using different tools, which gives people a good feel of what's possible in each and how "hard" or "easy" the coding itself is for the different options. Different chapters discuss different types of data that you could encounter: patterns over time, proportions, relationships, ...

There were some minor points in the book that I'd mention if they asked me to review it (but that's according to me, and I don't want to pretend to be an expert). First of all, it would have been nice if Nathan had gone a little bit deeper into theories behind what is seen as good visualization. In the first chapter ("Telling Stories with Data") he does mention Cleveland & McGill in a side-note, but I think that information (along with Gestalt laws, etc) definitely deserves one or two full paragraphs, if not half a chapter. I also don't completely agree with the use of a stacked barchart (about page 109). From my experience, they're worth less than the time it takes to create them. After all, it's impossible to compare any groups other than the one that is at the bottom (and therefore has a common "zero"-line). For example: look at the first picture below. This shows the number of "stupid things done" by women and men, stratified over 5 different groups (A-F). Although it is easy to compare total stupidity per group (group C is doing particularly bad), as well as that for men, we can't see which of the groups A, D or F scores the worst for women. And that's because they don't have a common origin. We could of course put the women next to the men, but then we'd loose the total numbers.

In the second plot, however, it is possible to compare women, men and totals. The bars for women are put next to those for men, but I've added a shaded larger bar at the back that shows the sum of the two. This plot was originally created in R using ggplot2, but I'm afraid I can't find back the reference that explained how to do this... Let me know if you can find it.

The contents of the book of course is not world-shattering. But that's not the point of the book. For people new to the field it's a great addition to their library (and I learned a thing or two myself as well). If you're interested in data visualization, go out and get it.

Visualizing the Tour de France

2011-07-26T14:48:00.000+02:00

UPDATE: I encountered a blog post by Martin Theus describing a very similar approach for looking at this same data (see here).

Disclaimer 1: This is a (very!) quick hack. No effort was put in it whatsoever regarding aesthetics, interactivity, scaling (e.g. in the barcharts), ... Just wanted to get a very broad view of what happened during the Tour de France (= biggest cycling event each year).
Disclaimer 2: I don't know anything about cycling. It was actually my wife who had to point out to me which riders could be interesting to highlight in the visualization. But that also meant that this could become interesting for me to learn something about the Tour.

Data was copied from the Tour de France website (e.g. for the 1st stage). Visualization was created in processing.

The parallel coordinate plot shows the standings of all riders over all 21 stages. No data was available for stage 2, because that was a team time-trial (so discard that one). At the top is the rider who came first, at the bottom who came last. Below the coordinate plot are little barcharts displaying the distribution in arrival time (in "number of seconds later than the winner") for all riders in that stage.

The highlighted riders are: Cavendish (red), Evans (orange), Gilbert (yellow), Andy Schleck (light blue) and Frank Schleck (dark blue).

So what was I able to learn from this?

Based on the barcharts you can guess which trips were in the mountains, and which weren't. You'd expect that the riders become much more separated in the mountains than on the flat. In the very last stage in Paris, for example, everyone seems to have arrived in one big group. Whereas for stages 12-14 the riders were much more spread. So my guess (and that's confirmed by checking this on the TourDeFrance website :-) is that those were mountain stages.
You can see clear groups of riders who behave the same. There is for example a clear group of riders who performed quite badly in stage 19 but much better in stage 20 (and bad in 21 again).
As the parallel coordinate plots were scaled according to the initial number of riders, we can clearly see how people left the Tour because the "bottom" of the later stages are empty.
We see that Cavendish (red) has very erratic performance. And it seems to co-incide with trips where the arrival times are spread out (= mountain trips?). This could mean that Cavendish is good on the flats, but bad in the mountains. Question to those who know something about cycling: is that true?
Philippe Gilbert started good (both on the flats and in the mountains), but became more erratic halfway through the Tour.

TenderNoise - visualizing noise levels

2011-07-13T20:11:00.001+02:00

A couple of days ago I bumped into this tweet by Benjamin Wiederkehr (@datavis): "Article: TenderNoise http://datavis.ch/q9pIxq" It describes a visualization by Stamen Design and others displaying noise levels at different intersections in San Francisco. They recorded these levels over a period of a few days in order to get an idea of auditory pollution. More information is here.

Although this particular visualization might be very useful for the people involved, I would like to explain some of the issues that I have with it, coming from a data-visualization-for-pattern-finding viewpoint.

I think there are many things that might be gleaned from this data which are not possible with the current visualization:

Is there a relationship between the noise patterns at different intersections? Based on the graphic at the bottom, we can conclude that on average noise level goes down during the night and up during daytime, but it would be nice if the visualization would give an indication of any aberrant patterns as well. Are there intersections that behave differently from others?
I don't see a real use for changing the graphic over time. I suspect that small multiples of area charts would work better to demonstrate the change over time (as e.g. the visual used here). Using the current approach it is very difficult to see how particular intersections change over time because (a) the display changes and you loose temporal context, and (b) the resolution is so hight that the blobs just flicker.
Concerning that flicker, it might be an option to bin the data in larger time blocks. For calculating the value in each block different approaches should be investigated, like the average value, the maximum, the minimum, or the most extreme value (be it maximum or minimum, based on comparison with the average).

It'd be interesting to get hold of these data and work on some alternatives (given the time...)

Why did I move into data visualization?

2011-06-30T16:28:00.000+02:00

Preamble: It's been very quiet on this blog since I left the Wellcome Trust Sanger Institute in the UK and took my position here at Leuven University in Belgium last October. Truth is that the type of work changed so profoundly that it takes a while to give it all a place in your head; let alone a blog. Until I remembered this morning why I started this blog in the first place: to help me order my thoughts in the first place. So it might have sped things up instead, actually...

Anyway... In this post I'd like to explain why I'm moving into the *data visualization* field. And it's not just because it's always nice to look at pretty pictures.

Statistics are great, but...

Ben Schneiderman, Professor at the University of Maryland, very eloquently stated that "The great fun of information visualization is that it gives you answers to questions you didn't know you had". I'd rephrase that a bit to "The great use of data visualization is that it gives you clues to questions you didn't know you had". To me, it's as much about finding the questions to ask as about finding the answers to those questions. Data visualization should not be used to "prove" things; that's what statistics is for. But the visualization can give you ideas on what statistical models to test. As do many others, I see a strong connection between statistics and data visualization. Taking a bit of a shortcut here, you could say that statistics is about proving what you expect, while visualization is about discovering what you didn't expect and refining those expectations.

From my own experience, I've seen that many (but not all!) statisticians look down upon data visualization with the argument that it can't proof anything. That's true. But their reaction then often becomes to throw away the baby with the bath water, instead of trying to see how both fields can benefit from each other. It's not always equally simple to convince people of the effectiveness of visualizations, but we're getting there...

Explain and explore

In the data visualization field, there is often the tension between explanation and exploration. The work I'll be doing here in Leuven will cover the whole spectrum. In the explanation corner, there is trying to make sense of complex data. For example helping cancer genetics researchers understand how tumours evolve (e.g. the phylogeny of cancer cells) or what the rearranged genome in those tumours looks like. This type of visualization sits downstream from the data analysis, after the data is churned. On the other hand, there are the exploration projects, where we focus on showing the raw(ish) data to help us decide on what type of analysis to perform, for example for investigating parameter-space for an algorithm. Of course many projects will fit somewhere in the middle...

The visualization model

Jarke van Wijk's paper "The Value of Visualization" (doi: 10.1109/VISUAL.2005.1532781) is a masterpiece in that it describes a comprehensive model of what visualization is and how we can quantify its effectiveness (cost). I'll just leave the picture here for you to contemplate over:

My inspiration

There are several people whose work I keep in mind when discussing what I want to do in my group; my sources of inspiration, so to speak. They have a very important thing in common: they don't take shortcuts in their work and are not afraid to really think about what their visualizations are intended to do.

These include:

Cydney Nielsen: ABySS-Explorer - a sequence assembly visualization tool
Miriah Meyer: Pathline - a tool for comparative functional genomics
Martin Krzywinski: Hive plots - rational network visualization

It's really exciting to work in this field; I'm looking forward to what the next few years will bring :-)

VizBi 2011 - looking back

2011-03-28T12:23:00.003+02:00

Has been a while (again) since my last post. It seems that the requirements on my time are just a little bit different from during my previous position... But I'd like to share a little bit about the VizBi conference that I attended 2 weeks ago.

This second installment of the VizBi conference was held at the Broad Institute in Cambridge, MA. This workshop focuses on visualizing biological data and harbours a new community of people with backgrounds ranging from bioinformatics and genomics to pure data visualization. Biological visualization is a very broad field, going all the way from abstract data visualization to enable finding patterns in data, to outreach for education, creating movies to explain things like DNA replication to the general public using Holywood-studio tools (for example molecularmovies and the work by Drew Berry).

Probably the best description of the relevance of this conference was given in the opening keynote by Eric Lander, professor at MIT, one of the people pulling the Human Genome Project and - entre autres - co-chair of the Council of Advisors on Science and Technology in the Obama administration. Some of his quotes:
- "If all data fitted nicely within one clear paradigm, we wouldn't need VizBi."
- "Things that are beautiful have huge communicative value"
- "Nowhere is visualization as important as in biology, because it's the bleeding edge right now, and has the messiest data and the messiest problems."

And my favourite: "We need to work on exchange rates; if one picture is worth only a thousand words, we're screwed."

I can highly recommend having a look at Tamara Munzner's keynote on visualization principles, available here.

Slides and videos of the talks will be available on the vizbi.org website.

Open Research Computation - a new journal from BioMedCentral

2010-12-13T10:28:00.003+01:00

As a colleague of mine said a couple of weeks ago: "if you don't publish it, it didn't happen". Scientific publications are the currency to advance a researcher's career. Looking for a new job? You better make sure your publication list is littered with first or second author papers in good (read: high impact factor) journals. Hoping to have your tenure track lead to tenure? Idem. Publish or perish.

Meanwhile, many bioinformaticians spend huge amounts of time developing software to make genetic or genomic research possible; research that just wouldn't happen if it was not for their custom-written tools, scripts and pipelines. Unfortunately, you often need the find function of your webbrowser or PDF reader to be able to pinpoint the lone bioinformatician in the author list. He's not the first or second author; he works in function of work by someone else.

Just like many others (see e.g. the alt-metrics manifesto), I feel quite strongly about the limitations of impact factors for judging a researcher's contribution to science. What about those papers that are published in "lower" journals but that are actually read more? What about blogs? What about your contributions on FriendFeed, seqanswers, quora, ...? How about all that software you wrote and that everyone can access on github? Don't you think these also have a teenie-weenie little value in the scientific discourse as well?

It will take time to change this skewed way of appreciating a researcher's work. We'll have to fight this from different angles, and subtly. We'll have to rethink what the term "publication" means, and how we can reference each other's work, including data that we made available in repositories and contributions to discussion forums. To start this process and give credit for the work of the computational scientist, BioMed Central launches a new journal today: Open Research Computation with Cameron Neylon at the helm as the editor-in-chief (and I'm grateful for having been asked to act as an editor). The aim of this journal is to provide a venue for programmers in research, who - as Cameron says it in his accompanying blog post - "either see themselves as software engineers or as researchers who code". The focus of the journal will be on making the tools available to the wider community of research programmers, with strict guidelines on code quality, code reusability and documentation.

I can only suggest you read Cameron's post and have a look at the www.openresearchcomputation.com website. The rest is up to you. Let's give the lone bioinformatician a venue to bring his work out in the public.

VCF, tab-delimited files and bioclojure

2010-08-12T18:08:00.005+02:00

A lot of the work I do involves extracting data from VCF files ("Variant Call Format"; see http://bit.ly/apUbi8). It's tab-delimited but not quite: some of the columns contains structured data rather than just a value, and the format of these columns might even be different for every single line.

An example line (with the header):

#CHROM POS   ID REF ALT QUAL   FILTER INFO                                            FORMAT   SAMPLE1
1      12345 .  A   G   249.00 0      MQ=23.66;DB;DP=89;MQ0=26;LowMQ=0.2921,0.2921,89 GT:DP:GQ 1/1:89:99.00

The INFO field is actually a list of tag/value pairs (except when it's just a tag), and the meaning of the data in the SAMPLE1 column is explained in the FORMAT column. Not only can different INFO tags be present on different lines, but the FORMAT can change line-by-line. Let it be clear that this is quite a bit of a pain to parse. What if I want to know the distribution of the depth at which each SNP is covered (i.e. the DP tag in the INFO field)?

So the first thing anyone does before working with such files: convert to "real" tab-delimited: the INFO field is spread over multiple columns and the same goes for the data in SAMPLE1. There is a vcf2tsv script available in vcftools, but it only extracts part of the data in the INFO field; not all tags in the INFO field will be represented in the tab-delimited output file.

I've tried to write a vcf2tsv method in clojure (see clojure.com and clojure.org). I'm now able to read/write VCF files without the tsv-intermediate using incanter (an R-like statistics environment using clojure), and also have a vcf2tsv command-line script that converts a VCF file to its tab-delimited counterpart, but this time including all the data that is in the file instead of a fixed selection of tags.

The vcf2tsv script convert the above sample data into:

CHROM POS   ID REF ALT QUAL   FILTER INFO-MQ INFO-DB INFO-MQ0 INFO-LowMQ       SAMPLE1-GT SAMPLE1-DP SAMPLE1-GQ
1     12345 .  A   G   249.00 0      23.66   1       26       0.2921,0.2921,89 1/1        89         99.00

So to generate that histogram of quality scores:

(use 'bioclojure)
(ns bioclojure)
(def snps (load-vcf "data.vcf"))
(with-data snps
  (view (histogram ($ :QUAL) :nbins 50)))

The result:

You can download the bioclojure library (which now only has this functionality :-) from github. That page also shows how to use the library. Feel free to fork and add new functionality; I only have time for very small incremental additions :-)

Postdoc position - Genomic variation discovery and visualization

2010-07-16T16:28:00.002+02:00

Just a short note...

Even though my position in Leuven only starts in October, I've already been involved in writing and defending a major grant. We've set up a consortium in Leuven (SymBioSys 2) consisting of 6 PIs "focusing on how individual genomic variation leads to disease through cascading effects across biological networks". This should be a good stepping stone to get my own lab running.

The grant concerns several workpackages (six in total), ranging from the analysis of the raw next-gen sequencing data, over the application of network algorithms for gene prioritization to the ultimate application in a couple of disease fields. I'll be in charge of the workpackage "Genomic variation discovery and visualization". Within this WP, we will apply and improve existing methods for the discovery, annotation and prioritization of SNPs and structural variation based on data available from the Leuven University hospital (HiSeq and 454). We'll develop a pipeline that can be used within and outside of the university. In addition - and I envision this to be the bigger part of the work - we'll work on visualization of these data. Several issues remain in this field: data can be too big to be readily visualized (e.g. locations of aberrantly mapped read pairs), or too complex (e.g. structural variation between individuals or family structures).

So I have a postdoc vacancy for someone to work on this project with a good understanding of genetics and next-generation sequencing. He/she should have good programming and statistics skill. Ideally, he/she will have visualization experience as well or be very interested to put his/her teeth into the subject. With visualization experience, I mean for statistical purposes (e.g. using R), but also more general (e.g. using Processing). Expect to be talking a lot about indexing methods, mapreduce, visual encoding, ...

The SymBioSys2 consortium consists of young scientists in the fields of bioinformatics, human genetics and cancer. As I'll only be starting my own group this fall the person doing the postdoc will have the chance to have a real impact on where my group will be going :-)

Encounter with incanter - about clojure, incanter and bioinformatics

2010-06-25T14:51:00.013+02:00

I have been a bit frustrated lately by the fact that for many of my analyses I have to write a ruby script to mangle my data first, then resort to R to add a statistic to each of the datapoints, go back to ruby to mangle the result, repeat, rinse, and finally make plots in R. Of course as a bioinformatician you're used to that and if necessary you write wrapper/pipeline scripts to handle this all for you if you know that this won't be the only time you have to do the analysis. But often you don't know. And breaking the analyses up into chunks just based on the tool you use will help you losing overview.

Being relatively new to R (less than 2 years), some of its idiosyncracies keep biting me. Just last week I had a dataframe to which I wanted to add a column which was the mean of two other columns. So from

1    9
8   16
2    4

I wanted to go to:

1    9    5
8   16   12
2    4    3

So you'd think (at least I would):

df$new_column <- mean(df$first_column, df$second_column)

and you'd be wrong. Just try it out... But it works for adding a column that is the sum of the other two, right? Don't get me wrong: many people are very good in R. Including the colleagues and ex-colleagues who I constantly ask for help. If you're reading this: you know who you are, and thank you.

There has been some talk on FriendFeed lately about the clojure programming language, which is a Lisp-inspired functional language that runs on the JVM. There is also incanter, a "clojure-based, R-like platform for statistical computing and graphics". Apart from the fact that functional languages are said to be good at working with huge data and concurrency, I was interested because, if incanter can do what I need from R, an incanter script is basically just a clojure script. So no switching between ruby/R/ruby/R/ruby/R; just do everything in the same place.

So I borrowed someone's "Programming Clojure" book, got onto some websites, downloaded incanter and gave it a spin this week. And my verdict: I'm impressed. It's a bit of a warp in your brain at first - being a functional language and all - but once you get the hang of it, it is easy to write. Of course I haven't done anything really difficult yet, but

You'll get used to the prefix notation ("+ 1 2") instead of the infix notation ("1 + 2"). At least I did. Maybe it helped that the calculator that I have been using for more than a decade - my trusted HP 42S - also doesn't use infix: "1 2 +".

So let's go for an example. Let's say I have a file with SNPs. For each SNP I have the position, a quality score and the reference and alternative alleles. Suppose I want to know what the ratio of transition over transversion is for a given set of quality score cutoffs. So the file looks like this:

1   123    29   A   G
1   245    93   C   G
1   832    51   C   T
1   1234   63   T   G
1   2345   75   A   C
1   8315    9   C   A
1   9213   59   T   G
...

For a list of quality score cutoffs (let's say [10,40,60,80,100]), I wanted to get the ratio of transitions to transversion. Output would look like this:

I'd be interested to know how people would do this in R. Please feel free to put solutions in the comments. I can only learn from them...

But how about doing this in clojure/incanter? Being new to the language, I first wrote a script that would add "ti" or "tv" as an extra column, and another one that would calculate ti/tv for those cutoff.

Step 1: Assigning "ti" or "tv" to SNPs

(use '(incanter core io charts))
(defn ti?
"Checks if two alleles constitute a transition"
[allele-1 allele-2]
(if (contains? #{"AG" "CT"} (apply str (sort [allele-1 allele-2]))) true false)
)

(defn tv?
"Checks if two alleles constitute a transversion"
[allele-1 allele-2]
(if (contains? #{"AC" "CG" "GT" "AT"} (apply str (sort [allele-1 allele-2]))) true false)
)

(defn ti-or-tv
"Returns 'ti' if two alleles are transition, otherwise returns 'tv'"
[allele-1 allele-2]
(if (ti? allele-1 allele-2) "ti"
 (if (tv? allele-1 allele-2) "tv" "other"))
)

(def data
(read-dataset "data_file.tsv" :header true :delim \tab)
)

(def titv-result (map #(ti-or-tv (:ref %) (:alt %)) (:rows data)))

(def data
(col-names
 (conj-cols data titv-result)
 [:chr :pos :qual :ref :alt :titv]))

(save data "data_file_with_titv.tsv")

This creates a file "data_file_with_titv.tsv" that has an additional column that either says "ti" or "tv".

Step 2: Calculate ti/tv using different quality score cutoffs

(use '(incanter core io charts))

(defn filter-qual
"Returns only rows of a dataset where quality score is better than cutoff"
[ds cutoff]
(sel ds :filter #(> (nth % 2) cutoff)))

(defn count-ti
"Count the number of rows that have 'ti' in the 6th column"
[ds]
(nrow (sel ds :filter #(= (nth % 5) "ti"))))

(defn count-tv
"Count the number of rows that have 'tv' in the 6th column"
[ds]
(nrow (sel ds :filter #(= (nth % 5) "tv"))))


(defn titv
"Calculate the ti/tv ratio in a given dataset"
[ds]
(float (/ (count-ti ds) (count-tv ds))))

(def filename
"data_file_with_titv.tsv")

(def data
(read-dataset filename :header true :delim \tab)
)

(apply max (sel data :cols :qual)) ; => max qual = 527.58

(def cutoffs (range 10 500 10))

(def titv-values (map titv (map #(filter-qual data %) cutoffs)))

(view (xy-plot cutoffs titv-values))

...and this nicely creates the graph that I want.

The code may seem strange, but once you get used to the prefix notation thing you'll be alright.

Of course there are downsides in using clojure/incanter as well. Being a young language, there is no application that can come close to R.app yet. There are plugins for eclipse, netbeans, vi and emacs, but I had trouble installing any of these. There is a one-click incanter Mac application, but it does not provide the same user-friendliness as R.app. My current (well: these few days) setup: I start an interactive clojure/incanter session

java -Xmx800m -cp /Library/Clojure/lib jline.ConsoleRunner clojure.main

and have a TextMate window open. I just write my script in TextMate and copy/paste to the command line. Not ideal, but workable.

I also don't like having to refer to column numbers in some occasions instead of always being able to refer to column names (see the "(nth % 2)" bits in the code). But this might be me: you can probably do that but I just don't know how yet.

Overall, I'd say that my first encounter with incanter has been a great success. The clojure language seems to be very consistent which greatly helps picking it up. You have to make a click in your head, but I'm sure any bioinformatician is easily capable of doing that. I would greatly recommend giving clojure and incanter a try. Get the book. Browse the web. Rewrite one of your small R scripts in incanter. Response on the incanter Google group is also very fast, including from the developer of incanter David Liebke (@liebke)

Further information can be found in the usual places: clojure.org, incanter.org and google. Also check out data-sorcery.org, the incanter blog.

Update: Tim Yates has created an R version of the incanter solution above. I'm impressed. There's no way I would have been able to write that... He's now in my "R guru" box. See here

Threads in ruby: probably not how to use them

2010-05-19T14:00:00.000+02:00

I should create an online labbook with code examples of how I do things. Keep going back to an example script I have to copy/paste the code for handling different threads in ruby. But I'll put it here for the moment :-)

Suppose I have a file with several millions of lines containing information on SNPs. And suppose I have a database that already contains data for those SNP. And suppose I want to update the entries in the database with the data from the input file.

Please note: this is a quick hack.

require 'rubygems'
require 'progressbar'

MAX_NR_OF_THREADS = 5
nr_of_lines = `wc -l input_file.tsv`.split(/ /)[0].to_i
pbar = ProgressBar.new('processing', nr_of_lines.to_f/MAX_NR_OF_THREADS)
File.open(input_file.tsv).each_slice(MAX_NR_OF_THREADS) do |slice|
  pbar.inc
  threads = Hash.new
  slice.each do |line|
    threads[line] = Thread.new do
      # do the actual line parsing, DB lookup and DB updates
    end
  end
  threads.values.each do |thread|
    thread.join
  end
end
pbar.finish

I know this is far from perfect:

I shouldn't need to create that array.
This way all concurrent threads wait for each other before the next slice is taken from the input file. If one of the 5 threads takes a really long time, the other ones will wait but could instead start parsing the next lines in the input file.

Don't think less of me for this code...

Tipping my toes in mongodb with ruby

2010-01-07T11:31:00.004+01:00

Read the excellent post by Neil Saunders on using ruby and mongodb to archive his posts on FriendFeed, prompting me to finally write down my own experiences with mongodb. So here goes...

Let's have a look at the pilot SNP data from the 1000genomes project. The data released in April 2009 contain lists of SNPs from a low-coverage sequencing effort in the CEU (European descent), YRI (African) and JPTCHB (Asian) populations. SNPs can be dowloaded from here; get the files called something.sites.2009_04.gz. The exercise that we'll be performing here, is to get an idea of how many SNPs are in common between those populations.

The input data
The input data contains chromosome, position, reference allele (based on reference sequence), alternative allele and allele frequency in that population. Some sample lines from CEU.sites.2009_04.gz:

1 223336 C G 0.017544
1 224176 C T 0.052632
1 224344 T A 0.824561
1 224419 C T 0.008772
1 224438 C T 0.412281
1 224472 T C 0.122807

What will it look like in a database?
If we would load this into a relational database, we could create a table with the following columns: chromosome, position, ceu_ref, ceu_alt, ceu_maf, yri_ref, yri_alt, yri_maf, jptchb_ref, jptchb_alt, jptchb_maf. However we would end up with a lot of NULLs in that table.

For example: there is a SNP at position 522,311 on chromosome one in the CEU and YRI populations, but not in the JPTCHB population.

chr  pos      ceu_maj   ceu_min   ceu_maf   yri_maj   yri_min   yri_maf   jptchb_maj  jptchb_min  jptchb_maf
1    223336   C         G         0.01754   G         C         0.47321   C           G            0.22034
1    522311   C         A         0.05263   C         A         0.33036   NULL        NULL         NULL

(In this example, I have already recoded ref/alt allele to major/minor allele. In some cases such as the SNP at 223,336 the major allele is not the same in each population.)

Actually, of the 21,742,359 SNPs, there are only 4,836,814 (about 22%) is present in each of these populations. Moreover, there are 13,957,866 (64%!) SNPs that are present in only one of the populations. Ergo: loads of NULLs.

This is where you can start thinking of using a document-oriented database for storing these SNP data: each document will be tailored to a specific SNP and will e.g. not refer to the JPTCHB population if it it not present in that population. Enter mongodb. Apparently the "mongo" comes from "humongouos". I hope it will live up to that name...

The above two SNPs could be represented in two json documents like this:

{
  "_id" : "1_223336",
  "chr" : "1",
  "pos" : 223336,
  "kg" : {
    "ceu" : {
      "major" : "C",
      "minor" : "G",
      "maf" : 0.017544
    },
    "yri" : {
      "maf" : 0.473214,
      "major" : "G",
      "minor" : "C"
    },
    "jptchb" : {
      "maf" : 0.220339,
      "major" : "C",
      "minor" : "G"
    }
  }
}

{
  "_id" : "1_522311",
  "chr" : "1",
  "pos" : 522311,
  "kg" : {
    "ceu" : {
      "major" : "C",
      "minor" : "A",
      "maf" : 0.05263
    },
    "yri" : {
      "maf" : 0.33036,
      "major" : "D",
      "minor" : "A"
    }
  }
}

Getting the data in the database
So how do we load this data into a mongo database? You'll obviously need a mongo daemon running, but I'll leave that as an exercise to the reader as it's very easy and well-explained on the mongodb website. My approach was to first import all CEU SNPs straight from a json-formatted file, and then parse the YRI and JPTCHB files

1. Loading the CEU SNPs
The mongoimport tool takes a json-formatted file and - guess what - imports it into a mongo database. Just reformat the lines in the original CEU.sites.2009_04 to json (see below) and run mongoimport --db test --collection snps --file ceu.json. Ceu.json looks like this:

{_id: "1_211", chr: "1", pos: 211, kg: {ceu: {major: "A", minor: "G", maf: 0.201754}}}
{_id: "1_216", chr: "1", pos: 216, kg: {ceu: {major: "T", minor: "A", maf: 0.035088}}}
{_id: "1_229", chr: "1", pos: 229, kg: {ceu: {major: "A", minor: "C", maf: 0.017544}}}

2. Loading the YRI and JPTCHB SNPs
Loading the YRI and JPTCHB SNPs is a bit trickier because we need to check if the SNP is already defined. If it's not we'll just add a new document, but if it is we have to update the existing document and add the population-specific details. For this I use a ruby script with the mongo gem (I already added an additional column to the input files containing a unique ID per SNP consisting of chromosome and position):


require 'rubygems'
require 'mongo'
require 'progressbar'

db = Mongo::Connection.new.db("test")
coll = db.collection('snps')

pbar = ProgressBar.new('processing', 13759844)
File.open('YRI.sites.2009_04.with_uid').each do |line|
  pbar.inc
  uid, chr, pos, allele1, allele2, allele2_freq = line.chomp.split(/\t/)
  allele2_freq = allele2_freq.to_f
  major, minor, maf = nil, nil, nil
  if allele2_freq.to_f <= 0.5
    major, minor, maf = allele1, allele2, allele2_freq
  else
    major, minor, maf = allele2, allele1, 1 - allele2_freq
  end
  
  snp = coll.find_one(:_id => uid)
  if snp.nil?
    snp = {:_id => uid, :chr => chr, :pos => pos.to_i, :kg => {:yri => {:major => major, :minor => minor, :maf => maf}}}
  else
    kg_data = snp['kg']
    kg_data['yri'] = {:major => major, :minor => minor, :maf => maf}
  end
  
  coll.save(snp)
end
pbar.finish

For each SNP, we first get the major and minor allele as well as the minor allele frequency based on the alternative allele frequency. Next, we search the entire SNP collection for the SNP ID (which consists of chromosome_position). If it doesn't exist, we create a completely new document, otherwise we fetch the "kg" data ("kg" stands for "1000genomes") and add a 'yri' key-value pair.

Do the same for JPTCHB and you should get a complete collection.

Summarizing the data
Going back to the question we're asking: how many SNPs are in common among these populations? Again, a ruby script:


require 'rubygems'
require 'mongo'

map = "function() { " +
 "var keys = [];" +
 "for ( item in this['kg'] ) { keys.push(item) }" +
 "emit(keys.sort().join(';'), {count: 1})" +
"}"
reduce = "function(key, values) { " +
 "var sum = 0; " +
 "values.forEach(function(doc) { " +
 " sum += doc.count; " +
 "}); " +
 "return {count: sum}; " +
"};"

db = Mongo::Connection.new.db("test")
coll = db.collection("snps")

result = coll.map_reduce(map, reduce)
result.find.to_a.each do |r|
  puts ['{', r['_id'], ':', r['value']['count'].to_i, '}'].join(" ")
end

This script takes 50 minutes to run using a mongo database on my MacBook laptop.

Counting things in a collection makes at least me think of mapreduce. And I had seen mapreduce mentioned on the mongodb website, so clearly wanted to give that a go. Unfortunately, you have to define the map and reduce functions in javascript, which is a bit unsightly within a ruby script, but so be it. For more info, see this blog post by Kyle Banker.

This shows us the following numbers of SNPs:

CEU alone: 2,954,418
YRI alone: 6,901,757
JPTCHB alone: 4,101,691
CEU/YRI: 915,476
CEU/JPTCHB: 926,406
YRI/JPTCHB: 1,105,797
all three: 4,836,814

As Pierre would say: that's it.

Trying out mapreduce - on the farm

2009-09-18T12:31:00.006+02:00

Received an email this week from Sanger helpdesk that they installed a test hadoop system on the farm with 2 nodes. Thanks guys! First thing to do, obviously, was to repeat the streaming mapreduce exercise I did on my own machine (see my previous post). Only difference with my local setup is that this time I had to handle HDFS.

As a recap from my previous post: I'll be running an equivalent of the following:

cat snps.txt | ruby snp_mapper.rb | sort | ruby snp_reducer.rb

Setting up
First thing the Sanger wiki told me was to format my HDFS space:

hadoop namenode -format

This apparently only affects my own space... After that I could start playing with hadoop.

Where am I?
It looks like hadoop installs its own complete filesystem: even if my path on the server would be /home/users/a/aerts, a

hadoop fs -lsr /

shows that in the HDFS system I'm at /user/aerts.

Preparing the run
First off: create a directory to work in. Let's call that "locustree" with two subdirectories, called "input" and "output".

hadoop fs -mkdir locustree
hadoop fs -mkdir locustree/input
hadoop fs -mkdir locustree/output

And copy your datafile to the input directory:

hadoop fs -put snps.txt locustree/input/

Running the run
Hadoop documentation mentions that any shebang line in the mapper and reducer scripts are likely to not work, so you have to call ruby explicitely. Provide both scripts as "-file" arguments to hadoop. Finally, from your local directory containing the scripts, run the following command to run the mapreduce job:

hadoop jar $HADOOP_HOME/contrib/streaming/hadoop-0.20.1-streaming.jar \
  -input locustree/input/snps.txt \
  -mapper "/usr/local/bin/ruby snp_mapper.rb" \
  -reducer "/usr/local/bin/ruby snp_reducer.rb" \
  -output locustree/output/snp_index \
  -file snp_mapper.rb \
  -file snp_reducer.rb

Et voila: a new directory snp_index now contains the file spit out by the snp_reducer.rb script.

Trying out mapreduce

2009-09-02T19:50:00.000+02:00

Photo by niv available from Flickr

I have long been interested in trying out mapreduce in my data pipelines. The Wellcome Trust Sanger Institute has several huge compute farms that I normally use, but they don't support mapreduce jobs. Quite understandable from the IT's and institute's point of view because it's a mammoth task to keep those things running. But it also means that I can't put a foot on the mapreduce path.

There are options to run mapreduce on your own, however. One is to install hadoop locally with the restriction that you can only use one node: your own machine. This walkthrough will get you started with that. Another approach is to use Amazon Elastic MapReduce (AWS EMR).

First of, I must admit that I'm far from familiar with mapreduce at this point in time. As far as I understand it's ideal for aggregating data when the algorithm can be run on a subset of the data as well as the complete dataset. The mapping step of the framework will create different subsets of that data, run the required program on it and return partial results. These results should be in the form of key/value-pairs and will serve as the input for a reduce step that combines all intermediate results into one.

What's held me back from using AWS EMR is the cost. It's practically nothing, but anything I run on their servers at this moment is paid for with my own money. So today I installed hadoop locally and started coding.

The pilot
One of the projects I'm working on is LocusTree: a way of indexing genomic features that allows for quick visualization at different resolutions. In this data structure the genome is first divided into bins of 4,000 bp (called level 0). The next step sees the creation of new "parent" bins that each hold 2 of the original ones (so each covering 8,000 bp; called level 1). You keep on grouping 2 bins into a parent bin until you end up with a single bin for a complete chromosome. For chromosome 1 that means 17 levels, I think.

Each of the features - for example a list of 15 million SNPs or log2ratios for 42 million probes - is put in the smallest node in which it fits entirely. A feature ranging from position 10 to 20 will fit in the first bin on level 0; but a feature from position 3,990 to 4,050 will have to go in the first bin on level 1. Subsequently a count and if applicable some aggregate value like sum, min or max is stored in that node as well as in every parent node.

So for each SNP/probe the indexing script has to find the smallest enclosing node and update the aggregate values for each parent up to level 3.

The setup
I'm not good at java so used the hadoop streaming approach for this indexing. The input file contains about 15 million SNPs and looks like this:

snp6216929   3   177718011   177718021
snp6216930   3   177718267   177718277
snp6216931   3   177718268   177718278
snp6216932   3   177718294   177718304
snp6216933   3   177718612   177718622
snp6216934   3   177718629   177718639
snp6216935   3   177718956   177718966
snp6216936   3   177719529   177719539
snp6216937   3   177719623   177719633

I want to get a list of nodes with the count of SNPs within that node and all child-nodes. Each node should also list the SNP names for which that node is the smallest enclosing node. So the output looks like:

3:0:44429   9   snp6216929,snp6216930,snp6216931,snp6216932,snp6216933,snp6216934,snp6216935,snp6216936,snp6216937
3:0:44430   1   snp6216938
3:1:22214   9
3:1:22215   1
3:2:11107   10
3:3:5553    10
3:4:2776    10
3:5:1388    10
3:6:694     10
3:7:347     10
3:8:173     10
3:9:86      10
3:10:43     10
3:11:21     10
3:12:10     10
3:13:5      10
3:14:2      10
3:15:1      10
3:16:0      10
3:17:0      10

First column is the node (chromosome - level number - start position), second column is the number of SNPs covered by that node, and third column lists the SNP names for the smallest enclosing nodes.

The scripts
Code for the mapper:

#!/usr/bin/ruby
CHROMOSOME_LENGTHS = {'1' => 247249719,
                    '2' => 242951149,
                    '3' => 199501827,
                    '4' => 191273063,
                    '5' => 180857866,
                    '6' => 170899992,
                    '7' => 158821424,
                    '8' => 146274826,
                    '9' => 140273252,
                    '10' => 135374737,
                    '11' => 134452384,
                    '12' => 132349534,
                    '13' => 114142980,
                    '14' => 106368585,
                    '15' => 100338915,
                    '16' => 88827254,
                    '17' => 78774742,
                    '18' => 76117153,
                    '19' => 63811651,
                    '20' => 62435964,
                    '21' => 46944323,
                    '22' => 49691432,
                    '23' => 154913754,
                    '24' => 57772954
                   }

def enclosing_node(chr,start,stop)
start = start.to_i
stop = stop.to_i
level_number = ((Math.log(CHROMOSOME_LENGTHS[chr]) - Math.log(4000)).to_f/Math.log(2)).floor + 1
resolution_at_level = 4000*(2**level_number)
previous_start_bin = 0
previous_level_number = level_number
start_bin = (start-1).div(resolution_at_level)
stop_bin = (stop-1).div(resolution_at_level)

ancestors = Array.new
while start_bin == stop_bin and level_number >= 0
  ancestors.push([chr, level_number + 1, previous_start_bin].join(":"))
  previous_start_bin = start_bin
  previous_level_number = level_number
  level_number -= 1
  resolution_at_level = 4000*(2**level_number)
  start_bin = (start-1).div(resolution_at_level)
  stop_bin = (stop-1).div(resolution_at_level)
end
return [[chr, level_number + 1, previous_start_bin].join(":"), ancestors.reverse].flatten
end

STDIN.each do |line|
name, chr, start, stop = line.chomp.split(/\t/)

encl_node, *ancestors = enclosing_node(chr,start,stop)
STDOUT.puts [encl_node, 1, name].join("\t")
ancestors.each do |node|
  STDOUT.puts [node, 1].join("\t")
end
end

Don't worry about the details. It basically takes the input from standard input and writes this to standard output:

3:0:44429   1   snp6216929
3:1:22214   1
3:2:11107   1
3:3:5553    1
3:4:2776    1
3:5:1388    1
3:6:694     1
3:7:347     1
3:8:173     1
3:9:86      1
3:10:43     1
3:11:21     1
3:12:10     1
3:13:5      1
3:14:2      1
3:15:1      1
3:16:0      1
3:17:0      1
3:0:44429   1   snp6216930
3:1:22214   1
3:2:11107   1
3:3:5553    1
...

The reducer script takes all this and summarizes it all:

#!/usr/bin/ruby
nodes = Hash.new
STDIN.each do |line|
node, count, ids = line.chomp.split(/\t/)
parsed_ids = ids.split(/,/) unless ids.nil?
if nodes.has_key?(node)
  nodes[node][:count] += count.to_i
  nodes[node][:accessions].push(parsed_ids) unless ids.nil?
else
  nodes[node] = Hash.new
  nodes[node][:count] = count.to_i
  nodes[node][:accessions] = Array.new
  nodes[node][:accessions].push(parsed_ids) unless ids.nil?
end
end
nodes.keys.each do |locus|
STDOUT.puts [locus, nodes[locus][:count], nodes[locus][:accessions].flatten.join(',')].join("\t")
end

Before running this under hadoop it's usefull to just manually pipe these together on the unix command line:

cat snps.txt | ruby snp_mapper.rb | ruby snp_reducer.rb | sort

Because this was the output I expected I could run the pipeline with hadoop:

$HADOOP_HOME/bin/hadoop
 jar $HADOOP_HOME/contrib/streaming/hadoop-0.18.3-streaming.jar
 -input snps.txt
 -mapper snp_mapper.rb
 -reducer snp_reducer.rb
 -output ./snp_index/
 -file snp_mapper.rb
 -file snp_reducer.rb

...and lo-and-behold: the snp_index directory contains the result I want.

Parting thoughts
This is the very first time I use mapreduce so the code is far from optimized. Note that the code above is not a reflection of my programming skills :-) I just wanted to focus on the hadoop pipeline rather than the code. Actually, running this thing on my own laptop lasted almost 4 hours and nearly brought it to its knees. It certainly beat it up badly: that hash in the reducer script is just way too big... This little mapreduce experiment is finished, but I might have another go later and optimize the mapper and reducer scripts (getting rid of that huge hash), or create some intermediate reducer step.

Unfortunately it's only using one compute node (my own laptop) so I couldn't check how fast it can be. Running a small subset on Amazon Elastic MapReduce shows that that should work as well.

I only kept track of the count as an aggregate in the nodes. In case I wanted to calculate sum, min and max, I'd probably use a json-formatted string in the key/value pairs:

{"count": 4, "sum": 17, "min": 2, "max": 12}

I will definitely use mapreduce in the future, given the premise that my work would either install it on a cluster locally or we get an institute account on AWS.

First test release of circular genome browser

2009-07-08T19:18:00.015+02:00

Worked a couple of days on pARP, the circular genome browser, and I think it's ready to be tested out by others. Consider this an alpha release: expect a lot of issues. It's easy to create regions with a negative length, for example. Also, I didn't focus yet on user-friendliness or general input files. Ways of interaction are not made clear to new users yet and the input files still need to have fixed names and be stored in a particular folder.

pARP is designed to be a genome browser for features that are linked to other features on a genome (e.g. readpair mappings). Using a circular display, lines can be drawn connecting these features.

pARP always shows the whole genome. You can zoom into selected regions but the rest is still shown albeit squeezed a bit more together. The reason for this is that I want to show the context at all times. Suppose you'd zoom into two regions A and B that are linked by a large number of readpairs. If the part of the genome that is not A or B is not shown any readpair that has only one of its reads in A or B will just not be shown. By showing the whole genome, even squeezed in a few pixels, you can at least see that some reads are linked outside of A and B.

I've put some information on the github wiki page, such as how to interact and what the datafiles should look like.

For a little taste: here's a very brief screencast:

A lot of things still need to happen:

Catch a lot of edge cases
Incorporate a library for fast loading of features (i.e. LocusTree, which doesn't exist yet)
Make interaction more straightforward: use mouse for panning/zooming for example
About 1,472 other things that I currently forget

Also: I'm looking for a new name for pARP. pARP stands for "processing abnormal readpairs" (which what is was meant for originally), but it's actually just a genome browser using a circular representation to show linked features. Suggestions I already got are encircle and SqWheel or Squeal (the last two based on sequence-wheel; Squeal was my own idea, so I like that most at the moment :-) ).

A very, very big thanks goes to Jeremy Ashkenas, the author of ruby-processing. With pARP I have been pushing the boundaries of what that library does, and he has adapted it for my needs as I went. See here for his ruby-processing library. Other thanks go to my colleagues Erin, Klaudia, Jon, Nelo and Chris for their ideas.

pARP can be downloaded or cloned from github. Mac, Windows and linux are available there as well.

LocusTree - searching genomic loci

2009-04-16T17:33:00.017+02:00

"Contigs should not know where they are." That's a phrase uttered by James Bonfield when presenting his work on gap5, the successor to gap4, a much-used assembly software suite. So you think: "Wait a second: you're talking about assembly, and the contigs should not store their position?"

This statement addresses a problem that we encounter often when working with genomic data: how to handle features. The approach often used is to give the feature a 'chromosome', 'start position' and 'stop position'. Seems reasonable, right? So if you want all features on chromosome 1 between positions 6,124,627 and 6,827,197 you just loop over all features and check if their range overlaps with this query range. Indeed: seems reasonable. Unless your collection of features goes into the millions. Suppose you have arrayCGH data with a resolution of 1 datapoint per 50 basepairs (the green/red bars in the picture). If you'd want to search for this region, you can focus on chromosome 1 (if you were so smart to create different collections per chromosome) and start comparing the start and stop position of each feature from the beginning. Chromosome 1 is almost 250Mb so would contain 5 million array probes. That's a lot of features to check.

What James alluded to, was that you should retain the position information at a higher level. He uses an adaption of an R-Tree, which is a datastructure that bins all raw data into bins of a certain size (say 25 elements). Those bins themselves are binned again into larger bins, and so forth until only one bin remains: the root. Each bin knows its start and stop position relative to the bin it's contained in. To search for the same region as above, you start at the root bin which (by definition) overlaps with that query. So you go to the next level and check each of the child bins of root. You can reject those that don't overlap with your range (the red bins) and only focus on the ones that do. Continue until you reach the bottom-most layer which contains your actual features.

This has two advantages for my work on pARP, the visualization tool for aberrant read pair mappings. First of all, it's faster to get those features in a given region. A major feature of the pARP tool is that you can zoom into any region. With more than 6 million values for readdepth, speed is a big issue. But that's not all. Think about this: if I have a screen resolution of 1200 pixels wide, why would I try and plot 6 million points next to each other? Twelve-hundred would be enough, isn't it? R-Tree to the rescue again. In building the tree, each bin stores some aggregate data from its members; for example the average readdepth. To draw the genome readdepth across the whole genome, I only have to start at the root of the LocusTree and go down to that level that contains 1200 bins. No need to load the complete dataset. At the moment LocusTree bins only take the average of whatever value its members have, but this can be expanded to contain any type of aggregation (e.g. number of objects).

LocusTree is not the fastest library around (being written in ruby, and by me), but I think it does what I need it to do. This means that it might not do exactly what you need it to do. But the code is at github, so feel free to fork and improve. I did add a list of to-dos at the bottom of the README file...

UPDATE: I have now realized that - due to DataMapper - this library does not work under jruby, which is necessary for pARP which uses ruby-processing. A better approach than the above would be to have C-based algorithms combined with an API, so that we don't need DataMapper.

The good and bad of genome viewers

2009-03-08T21:28:00.011+01:00

Back before the human genome was fully sequenced and NCBI, UCSC and Ensembl started working on visualization, it made a lot of sense to go for linear representations and use tracks for annotation. After all: chromosomes are linear. Using different tracks to show different types of annotation is the next logical step.

But there is not just one human genome on earth; according to Wikipedia there's about 6.76 billion copies as of March 2009. So instead of talking about "the human genome" in those browsers, we talk about "the reference genome". Each person on earth is different, and so is each human genome. (That putting the reasoning on its head, but never mind).

Differences between humans such as SNPs and microsatellites can still be shown in the track-based browsers.

Things get more difficult when you're looking at structural variation. Structural variation messes up the backbone of the linear genome browser: you can't show differences between individuals in one straight line. Suppose you want to investigate a copy-number variation (CNV) and consult UCSC. You'd find tracks such as this:

Although this does give you quite some information on the CNV in question, it's not an adequate representation of what the different alleles actually look like. It also highlights another issue: the concept of "the reference genome". As more and more genomes are getting sequenced, is the one that was picked first the best for visualization and indeed, the reference? To be able to handle the different MHC haplotypes in Ensembl, for example, the database contains a table called "assembly_exceptions" that contains the alternative assemblies for each haplotype.

I believe that further down the line (although it might be quite a while) we might need to forget the whole notion of a reference genome. Two options come to mind. First of all, we could create an artificial reference that contains all sequence and let each real sequence we want to look at well, reference, that artificial assembly. That would mean that the different MHC haplotypes for example would all be in the same sequence. Similarly, copy-number variants containing let's say 3 to 8 copies would include all 8 in the mock-assembly. Unfortunately this still cannot cover structural variation like inter-chromosomal translocations. We can't build a single artificial assembly that would incorporate those. So here's the alternative: deBruijn graphs. Instead of creating a single linear representation of a reference, just let's not. We could use building blocks to build up each individual. Take a look at this picture:

Suppose that each block is a part of a chromosome and the red and blue lines represent the path to follow to build up the chromosome for a particular individual. In this picture the red individual misses a part of that chromosome that is present in the blue individual, and another part is inverted. Notice that we don't make any (arbitrary) decision on what is the reference sequence. By dragging the blocks we can either place all red connections on one line or all blue ones, making them look like a reference.

If we'd then add annotations to this picture like genes, we'd be able to display fusion genes. Suppose that the densely-striped block is on chromosome 7 in the red individual but on chromosome 12 in the blue one. If there's a gene on the right breakpoints we end up with a fusion gene.

Time permitting I'm going to investigate how useful this will be in projects like CNVs in the 1000genomes project.

Who-o-o are you? Who who? Who who?

2009-01-13T15:20:00.001+01:00

Image by Danny McL via Flickr

There’s been quite a lot of discussions going on lately about author identification: Raf Aerts’ correspondence piece in Nature (doi:10.1038/453979b), discussions on FriendFeed, ... The issue is that it can be hard to identify who the actual author of a paper is if their name is very common. If your name is Gudmundur Thorisson (“hi, mummi”) you’re in luck. But if you are a Li Y, Zhang L or even an Aerts J it’s a bit harder. Searching PubMed for “Aerts J” returns 299 papers. I surely don’t remember writing that many. I wish… So if a future employer would search pubmed for my name they will not get a list of my papers, but a list of papers by authors that have my name. Also, some of my papers mention jan.aerts@bbsrc.ac.uk as the contact email. Well: you’re out of luck, I’m afraid. That email address doesn’t exist anymore because I changed jobs.

The idea exists to call into life a unique ID for each author similar to the doi (“digital object identifier”) for a paper. Thomson Reuters have created ResearcherID, but because doi’s are handled through a not-for-profit CrossRef, let’s call the unique author ID a dsi (“digital scientist identifier”). This dsi can then be used by that scientist to identify himself wherever he needs.

Here I’ll try to explain how I think this could work.

But first of all: what are the prerequisites for a dsi-based environment? Obviously, the journals would need to request the dsi of authors on submission rather than just their names and email addresses. They are able to get names and email addresses through the dsi. And secondly, we need a service that assigns dsi’s and where scientists can update their details and add information.

The service/website

Let there be a website (for argument’s sake http://www.dsi.org) that assigns new dsi’s to new authors (only one dsi per author). So I could for example be dsi.12345. This service should have additional functionality such as list of contributions, curriculum vitae, contact details, network. It should also provide a homepage or profile page for each scientist listing at least the name, affiliation and literature list (i.e. what you would get from a PubMed search). So if you’d go to http://www.dsi.org/dsi.12345 you’d see at least my name, the address of the institute I work and a list of papers that I co-authored.

Getting a dsi

It’s critical that one researcher only gets one dsi. This is less than straightforward because I believe many researchers will not be interested enough in the whole identity story to even remember if they already had a dsi or not. So if I were to go to the dsi website and request an ID, the website would ask for my name first. It’d also ask if I used different names in author lists (e.g. I’m a woman, got married and started using my married name instead of my maiden name). Using that information the service would then search pubmed for papers that are authored by someone with my name (who might be me). It could present that list to me and ask if I’m actually that same person or not. This way we’d build up a minimal list of papers. That minimal list would then be checked against the dsi database to see if there isn’t already someone with my name who has claimed these papers. Logically that person would be me and it would appear that I already have a dsi. If no dsi has this name and these papers associated the new dsi can be assigned.

Contributions

A central service like this would be ideal for collaborators and possible employers to find out about contributions of a specific researcher to science. Instead of asking for author names and emails (the latter change over time anyway), a journal would ask for the dsi of all authors. If the paper gets accepted that journal would notify the dsi service to add that paper to the researchers publication list. But it goes further than just the papers. It’s a shame that researchers virtually only get marks for their published papers (Publish or Perish) and not for other contributions to scientific research. What about people who submit data to genome annotation databases? What about contributions to discussion in comments to blog posts, FriendFeed, ...? Setting up public databases? Writing APIs for scientific data? Think of a browser-button with which you could sign certain contributions – anywhere. Signing a contribution would add a link in your list of contributions in the dsi system.

It should obviously be possible to log into the dsi system and edit or remove contributions that you made. That one little API you wrote 5 years ago seemed so important then but you’ve come to see it as insignificant now, for example.

Contact details

People change employer, email, address and even name. So there’s a problem inherent in only listing email address and institute on a paper. Using the unique dsi for the authors would always point to that researcher no matter how many times he or she moved jobs or contact information. When a researcher has his contact details changed he would log onto the dsi service (we’ll come to this later) and change those data. Other people would then see those details on the researchers dsi page (http://www.dsi.org/dsi.12345), or if the researcher wants to keep them hidden send a message through the dsi service itself. The researcher’s email address does not have to be visible to the outside world.

Network

Even though you might not want to make your email address visible for the whole world, you wouldn’t mind if the people you know would see it. Your network. I think that a dsi service should contain capabilities like those from LinkedIn. You should be able to build a trusted network (with people that you know well). This network is another important pilar in your contribution to science.

There would ideally be different personas you could set for your profile. The default would for example be that your profile page would only show your name and papers. But you might also have a full profile that is only to researchers who are logged into the service and are not further than two steps away in your network. That extended profile might show your contact details (including email), contributions outside of papers (e.g. comments on blog posts) and curriculum vitae.

OpenID

The above explains how I would like to see the issue of identification solved. But there is also the problem of authentication. How do I prove that I am dsi.12345? Ideally the dsi service would be an OpenID provider so that it let’s me prove that I own http://www.dsi.org/dsi.12345. Hopefully more and more websites (biomedcentral, nature, ...) would allow logging in using OpenID.

Apart from serving as an OpenID provider, the dsi service should obviously also be an OpenID consumer so I don’t have to remember another username and password but can use http://jandot.myopenid.com or http://saaientist.blogspot.com to log in.

I hope this gives a little bit of an idea of the environment I hope we’ll move to. Any comments welcome. Any progress even more…

To find structural variation, look at read pairs: introducing pARP

2009-01-06T18:18:00.027+01:00

Nextgen sequencing is making a huge impact on how research is done in the genomics field. One of the ways to discover structural variants in a genome for example is to create a clone library for an individual, sequence the ends of those clones and then map those ends to the reference genome. Suppose that the clones in the library are all 150kb large, then we would expect the ends of each clone to be mapped about 150kb from each other on that reference genome, in a forward/reverse direction. Any read pair that does not follow this pattern, might indicate a structural variation. There are of course numerous spurious mapping results, so we need to ignore those.

Suppose that the resulting data look like this:


1    1016287    1    1025027     FF     10
1    54809626   1    54814724    RR     20
1    65970649   1    67123551    DIST   32
1    143840263  1    143841351   RR     34
1    241524162  16   298176281   DIST   36

First two columns are the position of the first read from the pair; third and fourth columns refer to the second read from the pair. Fifth column is FF, RR or DIST: forward-forward, reverse-reverse or distance (i.e. >> 150kb). The last column is some arbitrary quality score assigned to the mapping of this read pair. Notice that the last of these lines shows a readpair where one end is mapped on chr1 and the other is mapped to chr16.

We can do two things: analyze and then create a picture, or create a picture and then interpret (see also one of my previous posts). In the first approach, you'd run a statistical analysis to see if certain regions have a higher prevalence of abnormally mapped read pairs. In the second, you plot the raw data and try to identify abnormalities by eye. Of course ideally you switch between both approaches.

To visualize raw read pair information I've written a tool called pARP (Processing Abnormal ReadPairs) and available from github. It's very similar to the display used by [edited] this paper by Hampton et al to display structural variation using Circos (see picture, taken from the circos website). But instead of just creating a static picture, pARP is meant to be an interactive tool to browse the data.

Below is a screenshot of pARP running on some test data. It doesn't look as nice as the above image, but remember that this is interactive and thus doesn't have minutes to calculate everything.

Some of the features:

pARP can display abnormal readpairs (forward/forward, reverse/reverse or wrong distance), read depth and other features (e.g. segmental duplications).
Circular display gives overview of between-chromosome mapped readpairs.
Chromosomes can be dragged from the circular display to the upper or lower linear display to show (a) more detail and (b) within-chromosome aberrant readpairs (note: none in the image above).
Visible readpairs can be filtered by quality score.
Readpairs that are close to the mouse position are highlighted.

Prefiltering of the data should be minimal, and only focussed on getting the amount of data down. For example, the readpair data file could contain all normal readpair mappings, but getting rid of those just makes the display much more visually clear and reduces the amount of data to be loaded by several orders of magnitude (obviously...).

The version just released (tagged v0.8) is workable, but not ready for prime time yet. At this moment the user has to run the tool using jruby instead of just loading it as an applet. Also the filenames to be loaded have to be changed in the parp.rb code itself. I hope to add functionality so that you can upload your own data into an applet, or use a URI to link to it. But can't promise because other work is waiting. So here's also a call for help: if you're interested in contributing, please do! There's a "features-yet-to-be-implemented" list further down.

Features not yet implemented:

pARP should be available as an applet/application.
User should be able to point to files or URIs representing files instead of changing filenames in the code itself.
Saving an image to disk (also from the applet).
Further performance improvements.
Fixing of not-yet-identified-but-definitely-present bugs.

And now for some technical stuff. To keep redrawing times low so that the interaction wouldn't suffer too much from the huge amount of data, I had to use a few tricks. First of all, pARP makes heavy use of buffers. Different parts of the image are stored on different buffers. When the user interacts with the display, only the relevant buffers are updated while the others are untouched. For more info, see the github wiki page on the subject. Secondly, I've found out how to use ruby threads to load some data asynchronously. In particular the readdepth data can be a huge hog on performance; there are >6 million datapoints for a genome window size of 500bp. So what happens is that (a) readdepth data for a chromosome is only loaded when that chromosome is displayed in the linear part of the image, and (b) the readdepth data is drawn onto a separate buffer that is only displayed when the thread is finished.

Many thanks to:

Ben Fry and Casey Reas for Processing
Jeremy Ashkenas for the ruby API to Processing

Update: reference changed for Circos picture

Visualize or summarize?

2008-11-12T18:25:00.012+01:00

Image by Kaeru via Flickr

I've recently started using raw visualizations to get an idea of what data looks like rather than writing scripts to summarize. And what I found is that presenting data visually in a raw format might be more useful than condensing everything down into just a few numbers. Trouble is that you need to know what you expect and make assumptions if you want to analyze the data. The best tool you have for identifying trends or non-randomness is yourself, not R or a scripting language.

Bought the Visualizing Data book by Ben Fry to help me with this. It explains how to use the Processing language to present data in a meaningful way. As far as I understand, Processing is a wrapper around the java language so that it becomes much more intuitive to use for simple people like me. The language is so easy that there was only a very small learning curve for me, even though I didn't know anything about java other than that it's an island and a coffee. In several of my projects I now start with writing a simple processing script and then throw all my data at it. No assumptions made. The fact that it's easy to interact with a display with mouse or keyboard makes it even more useful.

The processing code editor makes creating a java applet or application a matter of one click, so it's easy to make your displays available for other people.

It's only after having a look at the data that I write analysis scripts or help other people in deciding how to analyze/summarize. That analysis can then help to for example make a more opinionated display. So it's display, analyze, display, analyze, ... as in a hermeneutic circle.

One issue with this approach is that you have to be able to think of a meaningful display. And I must say that's often (but not always) the more difficult bit. I started following the RSS feeds of some visualization blogs like FlowingData as well as the Processing website itself to get exposed to different types of visualization, which does help.

Update: Follow the discussion on FriendFeed.

Data visualization

2008-09-13T13:00:00.000+02:00

Today is "Data management, mining, curation and visualization" day at the Genome Informatics conference in Hinxton. It might be one of the more interesting ones for me, because that's what I do: manage, mine, curate and attempt to visualize. And I must say the last bit the most difficult. It's not difficult to upload results into a genome browser, but is it the best way?

I say we have to break free from the track. Ninety percent of all visualizations today in genomics is track-based (add DAS tracks to Ensembl, upload BED files to UCSC or run your own gbrowse). It's ideal for showing features on a chromosome, but it's used even if it's not the optimal tool (a feature shared with Microsoft Excel, but let's not go there). Why's that? Because that's what there is, and they do provide very useful functionality. But at the same time, having them available tempers the search for new and innovative ways of visualizing data. Having a computer at hand doesn't help either, I think: it's just much easier in PowerPoint to draw a collection of squares than a rich multi-facetted picture. That's just more easily done by hand, but that's not what we do, is it?

One of the articles in Nature's Big Data issue calls for artists and visualization experts to be involved before all data are gathered. This idea got quite a few comments on FriendFeed as well. I do agree with the idea of visualization experts being involved in many projects, but that visualization expert should be you. Well... you don't need to be an expert, but still you should have an idea on how to show the gist of your results. I think that's one of the important things that's missing in MSc education (apart with a good introduction to data management): some course in visualization concepts. How do you visualize time-series? How do you visualize differences? And what about time-series of differences?

Small example: I've been asked to think about how to visualize copy-number variations between individuals. The most obvious to do is what's used on the UCSC and any track browser: show a box where the variation is. But it's a variation, right? So what does this box mean? That some individuals miss that bit? That it's duplicated? What individuals? Using a track-based genome browser, you must make one individual the reference.

This is not about tools (there's Processing), but about a mind shift.

As it happens, I hope to get hold of a small tablet in the next week or so to replace my mouse and relieve my RSI a bit, so that might be a good opportunity for me to at least explore a bit.

Keep drawing.

Using git to sync server with laptop

2008-08-26T20:55:00.004+02:00

After investigating git for the bioruby project, I started using it on basically every project I run. And what do I use it for? Two things: keeping track of changes (duh) and syncing between server and laptop.

I normally try to get IT so far to let me mount my server Documents folder on my laptop when I'm at work. So ~/Documents actually points to my network drive. That's nice, because I don't have to bother with keeping track of several places to store my documents. If I change anything on my network drive, it looks like it's been changed locally. And vice versa.

But: what if I'm at home (where I work just a bit more than the wife would like)? I can still SSH into the server and do some work, but I can't mount that network drive. So I started creating a ~/LocalDocuments folder on my laptop in which I copied any files I needed. But that obviously feels wrong as I now have more than one place to put my files: either on my network drive (~/Documents) or locally on my laptop (~/LocalDocuments).

...until I started using git...

When I start a new project, I create a new folder on the server: ~/Documents/Projects/some_new_project. Within that folder, I run "git init" and commit a README file. This creates the git repository. Next thing: clone it on the laptop.

On the server:


mkdir /path_to_directory/some_new_project
cd /path_to_directory/some_new_project
git init
touch README
git commit -a -m "First commit"

On the laptop:


git clone ssh://my_name@network_server/path_to_directory/some_new_project

Now I just work in ~/LocalDocuments, commit all changes in my local git repository and (very important:) push it back onto the server.


git push

NB: This setup has already saved me from not a small (well: medium) disaster. For some reason (no coffee yet?) the very first thing I did one morning was login on the server, go to the project folder I had been working on for about a month, and do a "rm -r -f this_project". Aaargh! After wiping away all that cold sweat I realized I only had to clone the repository on my laptop back onto the server.

Bioruby with git: how would that work?

2008-06-23T12:37:00.022+02:00

Disclaimer: This blog post is the result of several iterations of writing/discussion/rewriting from Anthony Underwood, Michael Barton, Matt Wood and myself, with additional help from Paul Thornthwaite.
Disclaimer nr 2: We are not yet git veterans ourselves, so if you see simpler ways of doing what we describe below (or spot any errors), please let us know so we can update this post and put it onto the bioruby wiki as well.
Disclaimer nr 3: This is a proposal. Bioruby has not moved to git yet. However, we are working on it and trying to get the support from the main developers. Update: bioruby has been converted to git (thanks, Anthony) and is not available on github. So you can clone or fork now. However, the official development is still on CVS.

Update: I have discovered a very good presentation on how to work and collaborate with git. If you're interested in using git, have a look at this talk. You can fast forward to 1hr10min27sec where he starts talking about the practical use. Very strongly recommended.

In this blog post, we try to give some guidelines on how people can contribute to the bioruby code if/when that code will become available on github. The rationale for what we describe here is very much based on the premise that the job for the maintainer(s) of bioruby should be as simple as possible. Their workload should be as light as possible; this means that there are some additional steps that any contributor has to go through.
What follows is only a proposal. This is not a standard operating procedure; it’s only a guideline. Feel free to digress from it or use a completely different workflow. But remember: keep it simple for the maintainers.

Git

Distributed source control. Git is a truly distributed source control system, and in contrast with CVS or SVN, there is no central repository. With CVS or SVN, every time someone checks out or exports the repository, his own copy is so-to-speak subordinate to the central one. Not so with git: every single clone is equivalent; none is more important than another. In technical terms, the copy of bioruby on your laptop is as important as the one that for example Toshiaki maintains. One of the big advantages is that continued support is more likely should a key developer move on to pastures new (or github goes up in smoke), since the community can simply elect a new "blessed" repository (see below).

A blessed repository. Noticed that I said “in technical terms”? In some cases, like for bioruby, we would obviously like to have some repository that we would consider the ‘true’ one. Enter the notion of a “blessed” repository. This is purely by convention: the community appoints one particular repository as the main one.
A good place to put this repository is Github. For bioruby, this blessed repository will start out to be http://github.com/bioruby/bioruby. Official bioruby builds will take place from there. However, development can take place in additional, personal repositories.

Forking. Any development of bioruby would happen in clones of this blessed repository. Using the “fork” button on Github not only creates a clone, but it automatically puts that clone on Github itself as well. (Forking has the added value of the github social aspect where the network of changes can be viewed.) So if I would want to contribute, I would fork from bioruby/bioruby (that is: username/projectname) which would automatically create http://github.com/jandot/bioruby.

Guidelines for contribution

There are several ways of contributing: you can either create a patch or use a fork/clone.

Here we’ll try to explain how contribution could work with forking for Bioruby, both from the individual contributor’s view as from the view of the person(s) managing the blessed repository. What follows is not a Standard Operating Procedure. You do not have to do it like this. However, it will make it easier on the blessed maintainers to merge your code.

A. Using patches

A.1 The contributor

The simplest way to contribute is to send in patches. RailsCasts has a great screencast explaining this.

Creating a fork
Click the “fork” button on the bioruby/bioruby page. This will create a new repository in your own namespace: jandot/bioruby. It’s on this clone that you will be working; you will not touch bioruby/bioruby itself.

Making changes
To actually start making changes (e.g. you want to add functionality for Ensembl cigar format), you create a local clone on your own computer (step 2 in the picture):

git clone git@github.com:jandot/bioruby.git

This will be your local master branch. The first thing to do after cloning your own fork, is to create an additional branch for the feature you want to work on: add_cigar_format (step 3)

The command to do this:

git checkout -b add_cigar_format

This will create the new branch and check it out so it becomes your active one. From the fluxbox wiki (http://fluxbox-wiki.org/index.php/Git_-_using): “Branching and merging is very powerful in git. You can create thousands of local branches, one for each bug you work on or feature you implement. It is good practice to do this because it safes your from accidentally pushing changes to another repository.”

So you’ll end up with 2 branches (do a “git branch”):

master: a reflection of the master branch of your remote repository
add_cigar_format: is where the actual work is done

The “git branch” should have a star in front of add_cigar_format because that’s your current branch. If master is starred, do a “git checkout add_cigar_format” to change to this branch.

Now you can edit and change to your heart’s content. The current branch you’re working on maintains an index of files that git is tracking. You can find the current status of the branch by typing

git status

which will list the current status of all the files. Changes can be committed to the local index by using the command.

git add file

The index is an intermediary between the working copy files you are editing, and the changes committed to the repositroy. Changes can be committed from the index to your local repository using the command

git commit

This command will also prompt you for a message describing the commit. Try not to do too much work before committing. A single commit should concern (part of) a single conceptual change with its tests. It’s good practice to commit often (and several commits per conceptual change), but do try not to mix different changes into one commit. This will make it harder afterwards if a commit has to be reverted.

Commits are applied to the only current checked out branch (i.e. add_cigar_format), and do not affect any other branches, or the original repository. Also, if you have to make site-specific changes (e.g. hard-coding a proxy server in one of the files), try to put those changes in one single commit. This will make it easier later to remove them.

Preparing the patch
When you think your change is ready for inclusion in the blessed repository (and you’ve included tests as well), you can create a patch file. To make sure that the blessed repository maintainers will have no problem merging your version, you will want to make the patch reflect the latest version of the blessed repository (step 5).

git remote add blessed git://github.com/bioruby/bioruby.git
git fetch blessed

So now you can check that the patch you will submit will only contain the changes that you want to be included in the blessed repository. One of the things to look out for is that there are not site-specific configurations in your branch (e.g. a hard-coded proxy or directory path, no “STDERR.puts”, ...). Hopefully, you put all those site-specific changes in a separate commit as described above. To get rid of them, you just revert that commit. “git log” will show you the SHA1 of that particular commit (the long crazy string), and you just run “git revert [that_SHA1]". After that, check your changes:

git log -p blessed..feature_c

When that’s done, you can create the actual patch (step 6):

git format-patch blessed..feature_C

This creates a file that you can send to the maintainer (step 7). And you’re done…

A.2 The maintainer

The maintainer gets an email from someone containing a patch. The first thing to do, is to create a new branch and merge the changes into that branch.

git checkout -b feature_c
git am <0001-feature_C_commit_message.patch

Of course he would want to check the changes by comparing the new version of the code with the one that is in the blessed repository (i.e. the master).

git log feature_c..master

If everything looks OK, he can then merge the changes into master itself and push it up onto github.

git branch master
git merge feature_c
git push

And he’s ready. Only thing left to do is remove the branch he created during the process.

git branch -d feature_c

B. Using a pull request

B.1 The contributor

This type of contribution starts out exactly the same as the one with patches: you fork/clone, create a feature branch and hack away.

Preparing the pull
When you think your change is ready for inclusion in the blessed repository, you will create a branch specific for this pull (e.g. called to_pull; step 5): “git branch -b to_pull”.

To make sure that the blessed repository maintainers will have no problem merging your version, you have to rebase your branch (steps 6 and 7).

git remote add blessed git://github.com/bioruby/bioruby.git
git fetch blessed
git rebase blessed/master
git checkout blessed/master fileA_for_user_environment_only
git checkout blessed/master fileB_for_user_environment_only

At this point, a “git log -p blessed/master..to_pull” can help you check that the differences between your _to_pul_l branch and the blessed branch only contain the changes that you intend to be pulled (e.g. getting rid of “STDERR.puts” statements).

When you’re satisfied, you can put the to_pull branch onto your remote repository so it becomes available for the maintainers of the blessed repository (step 8):

git push origin to_pull:refs/heads/to_pull

and push the “Send pull request” button on github.

After that, wait for any news if your change is accepted or not. When your remote to_pull branch becomes obsolete, you can remove it (step 10) with

git push origin :to_pull

B.2 The maintainer

The first thing the maintainer has to do, is get the latest version of his own (i.e. the blessed) repository.

git clone git@github.com:bioruby/bioruby.git

Then he can get a copy of your to_pull branch:

git remote add your_name git://github.com/your_name/bioruby.git
git checkout -b your_name/to_pull
git pull your_name to_pull

...and check what the change looks like.

git log -p master..your_name/to_pull

If he’s satisfied, he can merge your changes into the blessed master branch.

git merge your_name/to_pull

If there are no conflicts, he can then push the new version up onto github:

git push

Useful links

The thing about git

Using git within a team -> must-read

A large but very informative and simple document

Git repositories for maintainers

Using git for samba development

Git howto

GitGuide: intermediate and advanced git

Who we are

Anthony

Mike

Matt

jan.

Would you want to contribute to a small open-source project?

2008-06-03T11:23:00.010+02:00

Just a quick plug to see if I can find people interested in helping me out in some of my projects.

In the last 2 years, I started four open source projects (well: the last one was today...), each of which scratches my own itch and does what it needs to do for me. However, some features will have to be added and bugs be fixed to make them more useful to others (you, that is...).

If you are using one of these projects, please think about contributing. With that all-new fancy git version control system, it should be simpler than ever to get your own copy, tweak things a bit and send the improvements back.

1. Bio::Graphics
The Bio::Graphics library "allows for drawing overviews of genomic regions, similar to the pictures drawn by gbrowse" (from the homepage). I believe dgtized (Charles Comstock) and I have done not a bad job in creating something that is really useful, but of course some features are still missing.

Some things that come to mind that need help:

Find a good description on how to install it on a Mac (I've just changed job and am trying to do that, resulting in a major headache and no bio-graphics still).
Add new features such as a type of track to show continuous data (e.g. GC-content). We're thinking about ways to implement this, but additional ideas are welcome to get that done.
Is someone working on a gbrowse-like application in ruby?
Although I haven't found it a bottleneck yet, it'd might be a good idea to look at the performance of the thing.

A full tutorial and more information can be found here. The actual code repository has been moved from rubyforge to github and can be downloaded at http://github.com/jandot/bio-graphics. To get your copy:

git clone git://github.com/jandot/bio-graphics.git

2. Ruby API to the Ensembl database
In late spring 2007, I started the ruby API to the Ensembl database. This API relies on ActiveRecord and very much tries to copy the functionality of the perl API (including transfer of coordinates between coordinate systems). Just recently I was glad to hear people at the Sanger here are looking/have looked into it as well. At the moment, only the core database is covered but it would be nice if the other ones (funcgen, variation) would be added as well. In addition, the API was developed based on Ensembl release 45. With a new release coming out every few months, the API has to be tested against those as well.

What needs help:

Add an API for the other databases, including variation and funcgen.
Keep the API testing going for each new release of the database.

The full tutorial can be found on rubyforge, but the source control has also moved to github at http://github.com/jandot/ruby-ensembl-api. Get your copy using git clone git://github.com/jandot/ruby-ensembl-api.git

3. Ruby API to the UCSC database
I just started out the API for the UCSC database. There is data available in that database that cannot be found in Ensembl yet, for example copy number variations. So I started a new project (now solely on github) by copy-paste-modifying some code of the Ensembl API. Unfortunately, there are 3415 tables in the hg18 database (yes, that's three thousand four hundred and fifteen). Obviously, I only created interfaces for the tables that I will need at work.

What needs help:

Add additional tables to the API.
Think about additional functionality that might be added to some of the models.
A tutorial.

Again, get your copy from http://github.com/jandot/ruby-ucsc-api.

4. Simple Project Logger
As mentioned in my previous post, I use a simple rails application to keep track of the things I work on; a digital labbook, basically. That application is called Simple Project Logger (or its unix name sprolog). It's only task is to allow me to create tasks within projects and log what I've done for each task. Sprolog works good enough for me at the moment, but there is some bad bug I can't get out. In addition, it just looks ugly.

So if you're interested:

The authentication doesn't work yet. You can login using OpenID, but it is still possible to view anyone's projects and tasks by just typing the full URL. I know it will need a "before_filter :login_required", but that just breaks the thing.
Sprolog could definitely use some CSS-love.

You can download sprolog from http://github.com/jandot/sprolog.

I'll buy you a beer.

Keeping track of things: using a labbook for bioinformatics

2008-05-20T10:55:00.012+02:00

It's been a while since my last post. Left my last job, was unemployed for a month (while still chairing a session at a conference), and just started my new position here at the Sanger Institute.

On every job you're able to pick up some new things that can help you out later. One of the good ones from Roslin was how to keep a labjournal for bioinformatics. In the position before Roslin (at Wageningen University in the Netherlands), I remember having trouble remembering what I did to my data. So I was really happy to see that they actually had thought about those things in Roslin...

So what's the problem?
There are very significant parallels between bench-based labwork and computer-based data mangling. In both, you take some input (e.g. an eppendorf filled to the brim with DNA, or a data file downloaded from the internet), you perform some actions on it (e.g. PCR on such and such temperatures, or a grep followed with a sort and uniq) to get some output (e.g. an agarose-gel picture, or a number). In the wet-lab world, it's almost always mandatory to keep a lab journal in which you write down where you got the DNA from, which concentrations you used from which chemicals and what voltage you use for running the gel. However, for people doing a little bit of scripting to get some data out of a big set of files for example, there often is no such obligation. "I just played around with the data", you'll hear. But they will need a mighty good memory if they are to recall what they did after a couple of weeks. Bioinformaticians (i.c. those who manipulate data) have the same obligation as any other researcher: your work should be described in enough detail that other researchers can repeat the steps to get to the same result.

Enter the SOP for bioinformatics written by my former PI (little wave to Andy). It has some really good suggestions for people involved in data handling, mangling and mining. In this post, I will try to highlight some of them. Note that this is not about application or API development, but about data. (I hope to post a new blog entry about using svn and/or git later about that).

The central tool used for recording bioinformatics work at my previous job was RT Task Tracker, a web-based tool to record tickets and often used for keeping track of helpdesk tickets. I found it a bit too big and having too many features for my own purposes however and decided to write some little application myself that would do just what I need: the Simple Project Logger (sprolog, I can plug this in my own blog, right?). Mind that although I use it at the moment it's still in alpha and full of bugs.

The main requirements of the recording workflow are:

If you get data from somewhere/someone else than yourself, record where/whom you got the data from. Of course, there might be updates of the files you downloaded from that FTP server a couple of months ago even though those files have the same name. To be able to tell afterwards, md5sums should be made of any downloaded files and files that were sent to you by email.
Any mangling of the data should be recorded. Stuff like "my_script.rb <> output.txt" and "grep 'abc' input.txt > output.txt".

So how does that work in practice?
The sprolog application I wrote has the concepts of project, task and step. A project is a, well, a project. For example: "build my house" or "sprolog". A task is some distinct thing you have to do within the project, e.g. "place the roof" or "add authentication". Each task is then completed by a number of steps ("phoned contractor", "installed acts_as_authenticated").
When starting a new project, I give that project its own subdirectory under ~/Documents/Projects/. In turn, each task gets its own subdirectory within that project, named using the following convention: date + sprolog ticket number + short description (e.g. "20080513_T4-5_GenerateOligosForNewArray"). All work for that task is performed within that directory.
While performing the work, I copy/paste all necessary steps in sprolog. Typical steps look like this:


Step -> Recorded at Tue May 13 15:42:37 +0100 2008:

Saved attachment of John Doe as #{project_dir}/his_data.xls

Extracted tab-delimited version for each chromosome, changed newlines and added # before header.

MD5 (his_data_chr1.txt) = e5c38a91d8e5a666488863099fc5ef1c
MD5 (his_data_chr10.txt) = 9a702fb1f31bec42ec87089fc77efcc5
MD5 (his_data_chr11.txt) = e1a9a63e5c016cf93cb08ea6a5e425e5
MD5 (his_data_chr12.txt) = 8ab2bf7032f56df93b8b10c78bc2e1d4
MD5 (his_data_chr13.txt) = c1b2d609956edcf80657ed5f90b9469c
MD5 (his_data_chr14.txt) = fa6bfda1cd4e76f797ed8bd88d508448
MD5 (his_data_chr15.txt) = 46dbd8de0916dd69e81c519ac05671fe
MD5 (his_data_chr16.txt) = 302d920c6bef199a4bf40cfa2171348f
MD5 (his_data_chr17.txt) = 9aee0113c96f919c0603da3ccb9fca44

MD5 (his_data_chr8.txt) = e0d38c6804e39cae883dedfc648a2cda
MD5 (his_data_chr9.txt) = 94dfc7abc08d2b143f4eb13f29cadbdb
MD5 (his_data_chrUn.txt) = 5275595d7dfd4d4eb664e6bc9b08398c
MD5 (his_data_chrX.txt) = d90ba7f40b1019e1bbf981d894268dbc
MD5 (his_data_chrY.txt) = b991ff6f92869dcdc7b39da71d4d4b16

Step -> Recorded at Tue May 13 15:58:43 +0100 2008:

Venter: email boss gives the conditions on how to select deletions in reference genome.

Just to make sure I've understand correctly, if I want to identify
features for which there is >1kb of non-N sequence for which the
reference sequence has the allele then I identify all sequences in the Excel file and filter on those that have >1000 non-N bases.

Wrote script filter_records.rb to run this filter.

Step -> Recorded at Wed May 14 10:41:53 +0100 2008:

Renamed filter_records.csv to filter_records_on_non_n_bases.rb

ruby ./filter_records_on_non_n_bases.rb > filtered_records.csv

Number of lines in output file: 4411

Next step: repeatmasking

Problem: we don’t have the HuRef sequences, so those have to be downloaded first.

Downloaded assembled HuRef chromosomes from ftp://ftp.ncbi.nih.gov/genomes/H_sapiens/Assembled_chromosomes/ into ~/Documents/DataRepository/HuRef/

Step -> Recorded at Wed May 14 11:05:26 +0100 2008:

MD5 (hs_alt_HuRef_chr1.fa.gz) = 684e628536fa87b96343f1fea6219328
MD5 (hs_alt_HuRef_chr10.fa.gz) = 02ec433e2b00811db98c77a2fff3d161
MD5 (hs_alt_HuRef_chr11.fa.gz) = cca2a7098ed4d706dc8af7c58a2b9807
MD5 (hs_alt_HuRef_chr12.fa.gz) = 244bfd0f3f26cc109132c5518b2a1fb3
...

Lab journal
Once a considerable number of steps is performed or a task is completed, I print them out and glue them into a paper lab journal. That might look a waste of paper and completely unnecessary because you've got everything in electronic format anyway. I might change that behaviour later, but for the moment I just like to browse through physical pages when I need to know what I did rather then having to look at a screen. It's also easy to add annotations on those paper pages as well.

Note: if anyone is interested in helping develop sprolog, please let me know.

UPDATE: sprolog is now hosted on github at http://github.com/jandot/sprolog. Development on rubyforge will stop. Get your own copy by cloning it:
git clone git://github.com/jandot/sprolog.git

Note: picture taken from http://www.flickr.com/photos/cdnphoto/301083106/

Where did I get that data from?

2008-03-12T17:57:00.010+01:00

Did you ever have data lying around that you couldn't figure out where you got it from?

You downloaded and imported data from an FTP site into your database ages ago and you actually want to use it now. But if different records come from different sources, it can be really challenging to know what data to trust or how to retrieve additional information afterwards. Not keeping track of the source of the data breaks the chain of provenance.

I've seen it happen.

So what do you do? The obvious thing when you're working with a database is to have an additional column for that table where you can store the FTP URL or the people who sent you the data. Easy peasy. Things get a lot more complicated if different cells in a single row can have different sources. For example: suppose one of your tables (called markers) contains information on STSs.


id  name      fw_primer     rev_primer    PCR product length   gene
---------------------------------------------------------------------
1   marker_1  AACCGGACGA    GACCTCGGAGAC         241           CYP2D6
2   marker_2  TCAATGGAGG    GATTCGCTGACTC        183           BRCA2
3   marker_3  CGCTATGACTGC  AACTGCGTCATG         221           DAG1
4    ...         ...          ...                ...           ...

Let's say that you get the primers and STS length from a couple of inputs (e.g. dbSTS as well as designed by colleagues: marker_1 can come from dbSTS while the other two were created by your colleague Tim) and the gene information was added by two different colleagues (let's say Bert and Pat). In this case it becomes quite impossible to have that information stored in an additional column of the table ("primers: dbSTS; gene: Pat"??).

I've discussed this issue with a database guy during the last Perl Programming course at CSHL, and guess what: there's no easy solution. Having it play in the back of my head for the last couple of months, I finally took the effort to actually draw some possibilities on a piece of paper. And here's what might do the trick: just add an additional table. Let's call it marker_sources. This table has exactly the same columns as the markers table (apart from a foreign key to the data table).


marker_id  name      fw_primer     rev_primer    PCR product length   gene
---------------------------------------------------------------------------
   1       dbSTS      dbSTS          dbSTS            dbSTS           Pat
   2        Tim        Tim            Tim              Tim            Pat
   3        Tim        Tim            Tim              Tim            Bert
   4        ...         ...          ...                ...           ...

Result: the resolution of your provenance increases up to single cell level. This does mean additional tables, but for typical use they will not have to be queried if you're just interrogating the data. But at least you can get to that information if needs be...

Testing small scripts

2008-02-26T10:20:00.011+01:00

Seasoned programmers know this: testing should be an integral part of developing any script/program/software suite. Part and parcel is the unit test, where you test every little aspect of your program little by little.

For larger projects using a bunch of library files, the setup for testing basically always looks the same: there's your /lib/ directory with your class definitions and your /test/unit/ directory which holds yours tests. See here and here for introductions on full-blown unit tests.

That's all nice and fine, but as a bioinformatician you often just write small scripts for which it would be way to much hassle to create those different directories and separate files containing the classes from those containing the tests. So what do we do?

Often, you end up running your program and looking out for part of the output that you know should be correct. Let's take a very simple example. Suppose we have a file with just one column that has numbers in it. The same number can occur multiple times, and the ground-breaking script you'll write will just count the occurrences of each. Even though there are thousands of lines, you know from visual inspection that there are 7 1's and 15 2's. So a script could look like this:


occurrences = Hash.new(0)
File.open('data.txt').each do |number|
  number.chomp!
  occurrences[number.to_i] += 1
end
occurrences.keys.each do |k|
  puts k.to_s + "\t" + occurrences[k].to_s
end

So you'd run the script and check if you get the expected values for 1 and 2. If not: revise, rerun and check again.

But this looks like something that would be ideally suited for a unit test, if it weren't for the fact that it'd be too much hassle creating those different files and all. What if we could put the testing code in the script itself?

Actually, with a few adjustments, that's not a problem. Look at the following version of the code.


class Parser
  attr_accessor :occurrences
  def run
    @occurrences = Hash.new(0)
    File.open('data.txt').each do |number|
      number.chomp!
      @occurrences[number.to_i] += 1
    end
  end
end

if ! $test
  p = Parser.new
  p.run
  p.occurrences.keys.each do |k|
    puts k.to_s + "\t" + p.occurrences[k].to_s
  end
else
 require 'test/unit'
  class TestSimple < Test::Unit::TestCase
    def test_simple
      p = Parser.new
      p.run
      assert_equal(p.occurrences[1],7)
      assert_equal(p.occurrences[2],15)
    end
  end
end

So what happened here?

The original script, as so many scripts we write, actually does two things: (1) it parses a file to extract information, and (2) it prints some things out. In these cases, we can take the approach I outlined here. The biggest change you have to make, is to put your code in a class, otherwise you won't be able to run the unit test. Secondly, the if ! $test separates out the behaviour of the code based on if you want it tested or just run. I'll explain this line later. But if the if ! $test is true, the script just dumps the same output as the first version. However, when that statement is false, the script loads test/unit and runs two test: checking if the value for 1 is 7 and the value for 2 is 15.

How does that if ! $test work? If you call your script using


ruby -s my_script.rb -test

instead of


ruby my_script.rb

ruby will provide your script with an extra (global) variable: $test. See man ruby for more information.

So with this approach you can use test-driven development also in your teenie weenie scripts and not just in your mammoth software suites.

Making Bio::Graphics extendable

2008-02-06T12:39:00.000+01:00

One of the issues in a library like Bio::Graphics, is the plethora of glyph types that users will want. Here's a little showcase of what's provided by the library:

Features on a DNA sequences can be represented as filled boxes, open boxes, boxes with arrows, lines, triangles, ... In this post, I'll show you (and remind myself) how I came to a version of the Bio::Graphics code that makes adding glyphs straightforward both by myself and the user. WARNING: this post is going to be rather technical... Sorry about that.

First pass
Suppose we want to make it possible to create a picture like this one:

You basically have to tell your script that marker features should be drawn as triangles, and both scaffold and clone features as coloured boxes. The initial version of doing the actual drawing looked like this (only taking the relevant bits):


class Feature
  def initialize(glyph = :generic)
    @glyph = glyph
  end
  attr_accessor :glyph

  def draw
    case @glyph
    when :generic
      drawing.rectangle(left, top, width, height).fill
    when :line
      drawing.move_to(left,top)
      drawing.line_to(right,top)
      drawing.stroke
    when :triangle
      # code to draw triangle
    end
  end
end

This does work, but you see the issue, right? Whenever I or someone else comes up with another idea on how to represent a particular feature, the library code itself has to be changed. So far from extendable, that is...

Second pass: extracting the glyphs
To handle this issue for perl's Bio::Graphics, Lincoln Stein uses the Factory pattern. Which means that he creates a single GlyphFactory object that spits out different Glyph objects for each feature based on the configuration set at the Feature level. As I didn't know a thing about Design Patterns (i.e. before Russ Olsen's "Design Patterns in Ruby" arrived here at work) I had no idea how to set something up like that and just started coding away. As it turns out, I actually implemented it using a Strategy pattern.

What I basically wanted, is to delegate the actual drawing of a feature to a glyph. The Design Patterns in Ruby book gives a good example for formatting text. Here's the code:


class XMLFormatter
  def output_report(title, text)
    puts('< xml>')
    puts('  < title>#{title}< /title>')
    puts('  < text>#{text}< /text>')
    puts('< /xml>')
  end
end

class PlainTextFormatter
  def output_report(title, text)
    puts("***** #{title} *****")
    puts text
  end
end

This can then be used in e.g. a Report class like this (also from the same book):


class Report
  attr_reader :title, :text
  attr_accessor :formatter

  def initialize(formatter)
    @title = 'Monthly Report'
    @text = 'Things are going pretty well.'
    @formatter = formatter
  end

  def output_report
    @formatter.output_report(@title, @text)
  end
end

Looks a lot like what we need, isn't it? Translating this to our purposes, the library code could look like this:


class Glyph::Common
  def initialize(caller)
    @caller = caller
  end
  attr_accessor :caller
end

class Glyph::Generic < Glyph::Common
  def draw(left, right, width, height)
    @caller.drawing.rectangle(left, top, width, height).fill
  end
end

class Glyph::Line < Glyph::Common
  def draw(left, right, width, height)
    @caller.drawing.move_to(left,top)
    @caller.drawing.line_to(right,top)
    @caller.drawing.stroke
  end
end

And use it in the Feature class like this:


class Feature
  def initialize(glyph_object = Glyph::Generic)
    @glyph_object = glyph_object.new(self)
  end
  attr_accessor :glyph_object

  def draw
    @glyph_object.draw
  end
end

At least this approach splits out the actual drawing into different simple classes. But the extendability still isn't there: the user still has to open the library file containing all glyph definitions and hack away in there.

Third pass: loading glyphs automatically
It's be nice if we could add new glyph types on the fly just by creating a little file containing the code for that glyph's class. Convention over configuration to the rescue...

What I did, was create a folder (/lib/bio/graphics/glyphs/) that contains the description of all glyphs in separate files:
generic.rb


class Glyph::Generic < Glyph::Common
  def draw(left, right, width, height)
    @caller.drawing.rectangle(left, top, width, height).fill
  end
end

line.rb


class Glyph::Line < Glyph::Common
  def draw(left, right, width, height)
    @caller.drawing.move_to(left,top)
    @caller.drawing.line_to(right,top)
    @caller.drawing.stroke
  end
end

So ideally, the only thing to make a script work that asks for a feature to be drawn as a empty box (feature = Feature.new(:empty_box)), would be to add a file to that directory called 'empty_box.rb'. Several things have to be taken care of to make that happen:
* loading the new file
* translating the :empty_box to EmptyBox

To load all files in that directory is easy enough. Adding the following code to the main bio-graphics.rb file (which loads the whole library) does the trick:


glyph_dir = File.dirname(__FILE__) + '/bio/graphics/glyphs/'
require glyph_dir + '/common.rb'
full_pattern = File.join(glyph_dir, '*.rb')
Dir.glob(full_pattern).each do |file|
  require file
end

To translate the :empty_box symbol into the EmptyBox class takes a little more work: we need to convert the snake_case symbol into a CamelCase string, and then create an object of the class that has that name. To do that, I extended the String class a bit with these additional methods:


class String
  def snake_case
   return self.to_s.gsub(/::/, '/').gsub(/([A-Z]+)([A-Z][a-z])/,'\1_\2').gsub(/([a-z\d])([A-Z])/,'\1_\2').tr("-", "_").downcase
  end
    
  def camel_case
    return self.to_s.gsub(/\/(.?)/) { "::" + $1.upcase }.gsub(/(^|_)(.)/) { $2.upcase }.to_s.gsub(/\/(.?)/) { "::" + $1.upcase }.gsub(/(^|_)(.)/) { $2.upcase }
  end

  def to_class
    parts = self.split(/::/)
    klass = Kernel
    parts.each do |part|
      klass = klass.const_get(part)
    end
    return klass
  end
end

Now what happens here? The snake_case and camel_case methods should be not that difficult to understand and are not really where the magic happens. The String#to_class method however is a different story. As it happens, every class in ruby is also represented by a constant (the class name always start with a capital). To get to the class that has the name MyClass, all you have to do is retrieve the constant with that name: Kernel.const_get("MyClass"). Unfortunately, having namespaces (Bio::Graphics::Glyph::Generic) makes things a bit difficult. You can't just do Kernel.const_get("Bio::Graphics::Glyph::Generic"). To get to the Generic class, you have to call the const_get method on the Bio::Graphics::Glyph class, which doesn't exist yet. Therefore we have to look through all parts of the namespace and build up the class as we go.

With this code in place, I rewrote the Feature class to use this functionality:


class Feature
  def initialize(glyph = :generic)
    @glyph = glyph
  end
  attr_accessor :glyph

  def draw
    glyph_name = 'Bio::Graphics::Glyph::' + glyph.to_s.camel_case
    glyph_class = glyph_name.to_class
    glyph = glyph_class.new(self)
    glyph.draw
  end
end

Now all a user has to do to add a new glyph type to his application, is:
* create a file in the lib/bio/graphics/glyphs/ directory that defines the glyph
* make sure that the name he gives to that class is the CamelCase version of the symbol he wants to use (which should be snake_case)

There you go. As I warned at the start: technical. At the moment this setup works for what I need the Bio::Graphics library to do. There might be a chance that the approach is changed in the future as we need to handle subfeatures, subsubfeatures, subsubsubfeatures, ... more elegantly. But thats' something for another post.

What makes code beautiful?

2007-11-26T04:37:00.000+01:00

Saw this webcast a couple of weeks ago where Marcel Molina explains the notion of beautiful code. And I really recommend anyone writing code to have a look at it (totally irrespective of the fact he uses a ruby example...). If you look at the progress bar of the presentation, it looks like it's really long, but only the first half is the actual presentation, followed by a lengthy discussion. (The last few minutes have a nice wrap-up by Chad Fowler.)

The talk is basically split into two parts: (a) what is beauty, and (b) how does that apply to coding?

What is beauty?

After going into the history of what philosophers like Plate, Rousseau and others found beautiful, Molina focusses on the three rules set out by Thomas Aquinas: proportion, integrity and clarity. Things that are beautiful comply to these rules. The examples that Molina uses:

proportion: suppose the fingers on your hand would be twice or half as long. Wouldn't really be beautiful, would it? Think 'golden ratio' here.
integrity: a thing has to be fit for its purpose. Take for example a hammer made of crystal. It might look beautiful, but can't do what it's supposed to. The appearance of beauty is not necessarily the same as real beauty.
clarity: it should be clear what something is/does.

How does it apply to coding?

It's straightforward to apply these principles to programming code. First, let's take a look at proportion. The code for simple functionality should be short. Not long. To add two numbers should be a single call. Molina gives an example here of how that had to be done in assembly. Really way too much coding going on to just do that. Here at work, the phrase "Jan can do that in 2 lines of ruby" are often heard (fortunately a little bit less, lately (pfew)). I feel that's not because I'm a code hacker, but because Ruby makes it simple to do simple things. Perl has the same characteristic, but at a price (I'll come back to that later).
Secondly, there's integrity. Does your code do what it's supposed to do? That's the most widely-used notion of good code: do the tests run? But code that does what it has to do can be still be ugly if it doesn't comply to the other principles. Spaghetti, anyone?
Thirdly, there's clarity. How much of your code can your fellow programmers understand by looking at it for less than a day? The fact itself that you have to explain your code might sometimes indicate that you're violating this rule.

One quirk of applying these principles to programming code, is that there's a clear trade-off between them. Some hacked together code might be very proportionate to what it's supposed to be, but in doing so it might become obfuscated and therefore violate the clarity principle. I must say that I've often found it difficult to understand my own perl code after a couple of weeks of not looking at it. Most of this can be attributed to the fact that the language itself needs you to write less clear code ($_, the way objects are implemented, ...).

When you see smelly code, chances are it's one of these three that are the problem.

Why does this matter?

To paraphrase Rousseau: Good code is beautiful code in action. Creating beautiful code (by the above principles) is no guarantee for good code, but it's a pretty good start...

Note: There's a book on Beautiful Code as well, and apparently there's a chapter by Lincoln Stein on his Bio::Graphics library. Better get my hands on that book for ideas for the ruby version.

Named arguments in ruby

2007-11-05T17:02:00.000+01:00

One of the main disadvantages of using ruby that I bump into is the absence of named arguments (or keyword parameters). That's no problem for methods taking just two or three arguments, but it does get confusing when you have to be able to pass more than that.

For example, the Bio::Graphics::Panel::Track::Feature#new method takes six arguments: the track it belongs to, a label, a type, the location, an array of subfeatures and a url. (Note: there's only four in the released version. Label and subfeatures are worked on at the moment...) When you get to this number of arguments, deciding what the order should be is almost more difficult than writing the actual code... You really have to start thinking about how the class will be used most often. If the user only has a track and a url, it's a bit unsightly to make him/her create a new feature by typing


my_feature = Feature.new(my_track, nil, nil, nil, [], 'http://www.google.com')

Compare that to


my_feature = Feature.new(:track => my_track, :url => 'http://www.google.com')

The way to get this to work, is to pass a hash in the method definition. There's a few gotchas, however. For instance, you can't define the defaults in the argument list directly.

The following code snippet almost works, but not quite...


class MyClass
  def initialize(options = {:name => 'unknown', :size => 0})
    @name = options[:name]
    @size = options[:size]
  end
  attr_accessor :name, :size
end

But you can set the defaults like this:


class MyClass
  def initialize(options = {})
    @name = options[:name] || 'unknown'
    @size = options[:size] || 0
  end
  attr_accessor :name, :size
end

In addition, you have to check if all passed keys make sense. What if a user would use the MyClass#initialize method above with a :age key? It's easy enough to catch that in your #initialize definition, but you have to remember to do that...

Another major disadvantage is the fact itself that you need a workaround in the first place. This means that it makes it a no-no if you want to integrate your classes in a bigger framework. There's no way I should use this approach for Bio::Graphics if it would be integrated with bioruby later. That would result in inconsistency in calling different classes within bioruby, which is the last thing you want...

Let's just hope ruby 2.0 will allow for using keyword parameters. (Pretty please...)

The state of bioruby (or: how can bioruby grow?)

2007-10-10T17:49:00.000+02:00

A number of people asked me recently about the usability of ruby/bioruby and if it would be worthwhile for them to take the plunge and investigate bioruby more. So I thought writing up here would be a good idea...

First a disclaimer: this is my own personal view on bioruby, based on experiences in the last year-and-a-half. In addition, this is about the bioruby project, not the code or the people.

Let's first see what bioruby does. It's a library of ruby classes and modules that can be used in biological -omics research. Just like ruby itself, bioruby's origins lie in Japan. Version 0.5.0 was released in 2003 and we're at 1.1.0 now. A brief and incomplete overview of what's covered by the library:

sequences
locations
pathways
alignments
trees
databases: GenBank, RefSeq, Ensembl, KEGG, ...
applications: fasta, BLAST, HMMER, clustalw, sim4, spidey, ...
...

Now why this post? Because I believe adoption and use of bioruby could be much improved. Is bioruby dead? Far from it. I think it's more like it is growing out of its clothes as any toddler does when it's getting older.

So what's the problem? It's not the quality of the code. It's not that too much stuff is missing. It's a sub-optimal level of communication (between users, from users to core-developers and from core-developers to users) and the low visibility of the project.

How can bioruby be taken forward? Somewhere in March, I sent some suggestions to the bioruby mailing list in response to a post by a very frustrated Trevor Wennblom. What it basically boils down to, is to get organized and get bioruby much more to the fore-front. So what options exist?

Getting bioruby organized
First of all, it wouldn't be bad if there would be a (mixed American/European/Japanese) board-like little group of people (3 or 4) who would be able to take the executive decisions on releases and what new modules should be incorporated in the bioruby library (after discussions on the mailing list, obviously). This would take a lot of weight of the shoulders of Toshiaki Katayama who now has almost single responsibility (and stress) for this. Having this done by a small group of people would relieve him from some of that stress.

Secondly, we need something of a playground for experimental modules. Call it bioruby-edge or whatever. These could be any modules that are not ready to go into the core bioruby, but are already really useful. When they reach a good enough quality, they can be moved to bioruby itself. There is already a bioruby-annex project at rubyforge, but according to its description it's only meant to hold rails plugins. (However, I was told to put my Ensembl API there as well...)

In addition, it would be good if the rubyforge project website would be used for feature requests and bug reports. This would then be the one-stop shop for the development.

And of course there is the documentation. I think we did a good thing in that big push to document the API in 2006, but the community needs more. The bioruby website already hosts the BioRuby in Anger documentation written by Toshiaki and Pjotr Prins. That's great stuff for quick lookups and I often use that information (especially the sequence IO; I never can remember. Must be an APOE4 mutation.). It would be nice though if it would be worked out a bit more. Take a look at the BioPerl documentation. I've always found the howto's really helpfull: getting a bit deeper into how the code works as well.
The rubyforge system provides wiki functionality for its projects, which is apparently not activated in bioruby. There is bioruby-doc maintained by Trevor, but I think it would be good to keep the core things together: put the wiki on rubyforge.

Letting people know about bioruby
Tremendous things are happening in the bioruby code (e.g. the rails thing), but we just don't know about it. Let alone what might be in store for the future. In addition, we know that the library exists, but the community who uses it has up till now not been really talkative about what they used it for. What we need here is communication in all directions.
First of all: a paper in a medium/high-profile journal. And sooner rather than later. This could serve as the starting block for building a wider bioruby community.

Secondly, we need to let each other know what we're doing with bioruby and how we're using it. I'm talking blogs and social networks here. I hope the blog you're reading at the moment might be a small contribution. Both the end-users and the core-developers should get their thoughts and work out in the open. Reports by end-users will keep the developers a bit on their toes and can highlight things that can be improved in bioruby. Core-developers could on the other hand shed a light on what's in store for bioruby in the future. How do you yourself use the code? Are you contemplating something great? We'd like to know what you're planning... I got a reply from Toshiaki about writing a blog, and he mentioned that it's not straightforward to do that in English. I do understand that that's a hurdle, but what I'd say to everyone having that issue: no problem. So let it be English "with hair on" (ooh, hair-rising-on-my-back translated literally from Dutch, but you hopefully get the idea). It's about us getting the big picture. Not about reading poetry. If we get the meaning, that's the main point.

The core-developers have done and are doing a great job. Respect. The only thing is that this toddler is now grown enough to want to play outside and will need additional clothes for that.

Ruby has so much to offer for bioinformatics as it has tremendous functionality and is yet so simple to code in. It would be a shame if the bioinformatics community can not capitalize on that.

Of course, I'd be very interested in your comments. Let's start talking! Especially about how to start that social network.

Note: while writing this entry, there were actually two messages sent to the bioruby mailing list asking for more documentation and easier access to new users (October 10, 2007). One of these messages stated "...BioRuby docs should have a version of more readable/easy-to-use format for beginners apart from the API stuff". Quod erat demonstrandum.

Second note: there's been some comments on the mailing list about the fact that this post was too much of a criticism to the original contributors to bioruby. That's not what I intended to do. Instead, it was my intention to look at options on how to take bioruby forward and let it grow from its small niche today to a more widely accepted toolkit. I've changed some phrasing in the text to hopefully make sure that that intention is clear (including changing the title).

Using rake to manage your software project

2007-10-09T11:28:00.000+02:00

Do you have some of those projects where you have to be sure that you jump the same loops every time you edit some code? Take a look at the bio-graphics code. Every time I change anything in the code, I have to do the following things:

regenerate the RDoc documentation
regenerate the ruby gem
check SVN status
do an SVN update
perform the SVN commit
upload the new documentation to the website

That's a prime candidate for rake. Rake does the same as GNU make, which is dependency-based programming. The major advantage for us over GNU make is of course that it uses ruby syntax. With dependency-based programming, I mean that some tasks rely on other ones. GNU make is best know for managing the compilation of source files. But you can do other stuff with it as well: if I want to commit to SVN, I want to make sure that the latest RDoc has been generated as well as a new gem. Therefore, you can have the 'SVN commit' task depend on the 'generate RDoc' and 'generate gem' tasks. And the task 'generate RDoc' will depend on the freshness of the actual library files.

How's this work? You basically create a file containing tasks and tell rake to execute one or more of them, the Rakefile. There are several good tutorials on rake, like the one from Martin Fowler and from the Rails Envy guys. I'm not going into the nitty-gritty of how they're written. These tutorials are much better at that. What I will do here, is describe the Rakefile I use for Bio::Graphics. (Someone already asked in the comments on my post on using ActiveRecord outside of rails what the Rakefile was that I used. Actually, the one mentioned in that post was empty and just a place holder.)

Without further ado, here it is:


# 
# Rakefile.rb
# 
# Copyright (C)::   Jan Aerts 
# License::         The Ruby License
# 
require 'rake'
require 'rake/testtask'
require 'rake/rdoctask'

task :default => :svn_commit

file_list = Dir.glob("lib/**/*.rb")

desc "Create RDoc documentation"
file 'doc/index.html' => file_list do
  puts "######## Creating RDoc documentation"
  system "rdoc --title 'Bio::Graphics documentation' -m TUTORIAL TUTORIAL README.DEV lib/"
end

desc "An alias for creating the RDoc documentation"
task :rdoc do
  Rake::Task['doc/index.html'].invoke
end

desc "Create a new gem"
file 'bio-graphics-1.0.gem' => file_list do
  puts "######## Creating new gem"
  system "gem build bio-graphics.gemspec"
end

desc "An alias for creating the gem"
task :create_gem do
  Rake::Task['bio-graphics-1.0.gem'].invoke
end

desc "Check SVN status"
task :check_svn_status do
  puts "######## Checking SVN status"
  message = String.new
  message << "# SVN status requires manual intervention\n"
  message << "# For items with '?': either svn add or svn propedit svn:ignore\n"
  message << "# For items with '~': don't know yet\n"
  message << "# Please see http://svnbook.red-bean.com/en/1.4/svn-book.html#svn.ref.svn.c.status"

  output = `svn status`
  puts output

  allowed_status = ['A','D','M','R','X','I'] # See http://svnbook.red-bean.com/en/1.4/svn-book.html#svn.ref.svn.c.status

  output.each do |line|
    status = line.slice(0,1)
    if ! allowed_status.include?(status)
      raise message
    end
  end
end

desc "Check if SVN updates available"
task :check_svn_update do
  puts "######## Checking SVN update"
  output = `svn update`
  puts output
  if output !~ /^At revision [0-9]/
    raise "Please update your working copy first"
  end
end

desc "Commit to SVN repository"
task :svn_commit => [:check_svn_update, :check_svn_status, :create_gem, :rdoc] do
  puts "######## Doing SVN commit"
  system 'svn commit'
end

rake -T lists all available tasks:

rake bio-graphics-1.0.gem  # Create a new gem
rake check_svn_status      # Check SVN status
rake check_svn_update      # Check if SVN updates available
rake create_gem            # An alias for creating the gem
rake doc/index.html        # Create RDoc documentation
rake rdoc                  # An alias for creating the RDoc documentation
rake svn_commit            # Commit to SVN repository

The file_list at the top contains all files in the library itself, and will be used to check all timestamps. The file 'doc/index.html' task looks at the timestamp of the index.html file and if it's older than any of the files in file_list, it will regenerate the documentation. If it's newer, nothing happens. Same goes for bio-graphics-1.0.gem.

The check_svn_update and check_svn_status tasks just check if subversion needs some manual intervention before being able to commit. This should be able to catch conflicts in the working copy and the repository, or files that you forgot to add the SVN.

Note: why didn't I use the special Rake::RDocTask instead of the one I use here? Because the built-in RDoc task first removes your whole doc directory, also deleting the subversion metadata contained in it.

Bio::Graphics and rails

2007-09-27T10:57:00.000+02:00

As a follow up to my post on Bio::Graphics, I tried integrating this library in a rails application. After all, you'd get your data either from a file (like GFF) or a database. And let me tell you: it took me just 30 minutes or so to get a proof-of-concept running. This included installing rails itself, creating the rails app, creating the database, loading dummy data, and doing the coding itself. That 30 minutes was interrupted for a couple of hours, because I needed some advice from Kouhei Sutou, the author of rcairo, on how to write PNG images in memory instead of to a file.

So how do you do it? The proof-of-concept little database I created contained 3 tables:

chromosomes (columns: id, name, length)
tracks (columns: id, name, glyph, colour)
features (columns: id, chromosome_id, track_id, name, location, url)

Create some features for a couple of different tracks for a particular chromosome.

In views/chromosomes/show.rhtml, add the following line:


<%= @chromosome.to_png %>

My models/chromosome.rb looks like this:


require 'stringio'
require 'base64'
require_gem 'bio-graphics'

class Chromosome < through =""> :features

def to_png(width = 800, start = 1, stop = self.length)
 return %{}
end

def draw(width, start, stop)
 panel = Bio::Graphics::Panel.new(self.length, width, false, start, stop)
 track_container = Hash.new
 self.tracks.each do |track|
   if ! track_container.has_key?(track.name)
     track_container[track.name] = panel.add_track(track.name, track.colour.split(',').collect{|i| i.to_i}, track.glyph)
   end
 end

 self.features.each do |feature|
   track_container[feature.track.name].add_feature(feature.name, feature.location)
 end

 output = StringIO.new
 panel.draw(output)
 return output.string
end
end

UPDATE: Apparently, Blogger does not allow me to paste the correct code above. In the to_png method, replace the following ascii codes:

%7B with {

%28 with (

%29 with )

%7D with }

And that's it. I leave the integration of my ensembl-api, bio-graphics and rails as an exercise for the reader. We could make a ruby version of the Ensembl browser... and then: world domination. Mwahaha.

Graphics, genomics and ruby

2007-09-20T14:03:00.001+02:00

Having known and used the Generic Genome Browser (aka gbrowse, see here) for years now, it occured to me a while ago that it should be o so simple to create the same functionality with a much easier setup if we could use ruby instead of perl.

Gbrowse depends on bioperl's Bio::Graphics module. Although gbrowse has been instrumental for many people's research, it does take a bit of work to get it installed. Apart from bioperl, it depends on Apache for showing the results in a browser. Compare that to any Rails application, where you basically just need ruby and a "gem install rails". I've created rails applications in the past that contain exactly the kind of data that would typically be visualized by something like gbrowse. Takes no time at all to set up and you can even get away by virtually writing no code. And no Apache to be installed, or configuration files that you can't access because you're not root.

Such a rails application makes it possible to browse, edit and delete the data. The problem comes with the visualization bit. There's no bioruby graphics library (yet?) that automatically parses features on a reference and creates a nice picture of where your genes are on that chromosome. Of course, the genes should be clickable so you can link through to NCBI or Ensembl.

I've spend some time in the last year creating such a Bio::Graphics thing for ruby. I wanted it to behave the same as the one from bioperl: there one panel that has one or more tracks, and each track has features on it. Even though it was quite easy to create a proof-of-concept library, the most difficult part was actually finding the right backend.

What should I use to create the pictures themselves? As I'd worked with SVG before, that seemed the right way to go. Downloaded a library from http://raa.ruby-lang.org/project/ruby-svg/ and got a prototype running quite easily. Problem: I needed an SVG viewer or firefox to actually view the picture, and zooming in/out screwed up all text. So after weeks of digging around, I've found rcairo, a ruby-binding to Cairo. Migrating to this backend was easy peasy and the pictures look really nice (see at the top). Unfortunately, it's impossible to create clickable glyphs using Cairo itself, but that can be easily worked around by creating a html file with the map. That's exactly what gbrowse does as well, isn't it?

The picture at the top has been created using the following simple script:


g = BioExt::Graphics::Panel.new(800, 1200, true, 1, 610)

track1 = g.add_track('generic')
track2 = g.add_track('directed',[0,1,0],'directed_generic')
track3 = g.add_track('triangle',[0.5, 0.5, 0.5],'triangle')
track4 = g.add_track('spliced',[1,0,0],'spliced')
track5 = g.add_track('directed_spliced',[1,0,1],'directed_spliced')

track1.add_feature('bla1','250..375', 'http://www.newsforge.com')
track1.add_feature('bla2','54..124', 'http://www.thearkdb.org')
track1.add_feature('bla3','100..449', 'http://www.google.com')

track2.add_feature('bla4','50..60', 'http://www.google.com')
track2.add_feature('bla5','complement(80..120)', 'http://www.sourceforge.net')

track3.add_feature('piep','56')
track3.add_feature('bla','103', 'http://digg.com')

track4.add_feature('gene1','join(34..52,109..183)','http://news.bbc.co.uk')
track4.add_feature('gene2','complement(join(170..231,264..299,350..360,409..445))')
track4.add_feature('gene3','join(134..152,209..283)')

track5.add_feature('gene1','join(34..52,109..183)', 'http://www.vrtnieuws.net')
track5.add_feature('gene2','complement(join(170..231,264..299,350..360,409..445))','http://www.roslin.ac.uk')
track5.add_feature('gene3','join(134..152,209..283)')

g.draw('my_panel.png')

What happens here?
Line 1: Create a new panel for a sequence of 800 bp, with the picture being 1200 points wide. Make all glyphs clickable if a URL is defined (the true), and zoom into the region from 1 to 610 bp.
Lines 3-6: Create different tracks, each with a name, a colour (in RGB at the moment) and a type.
Lines 8-24: Add features to those tracks, each with a name, a locus and an optional URL to link out to external websites. Notice how it handles spliced features and features on the reverse strand?
Line 26: Create the PNG (and in this case: also HTML) file.

Here's a nicer way to produce the same type of output:


#Initialize graphic for a nucleotide sequence of 600 bp
my_panel = BioExt::Graphics::Panel.new(1000, 1200, false, 1, 600)

#Create and configure tracks
track_SNP = my_panel.add_track('SNP')
track_gene = my_panel.add_track('gene')
track_transcript = my_panel.add_track('transcript')

track_SNP.feature_colour = [1,0,0]
track_SNP.feature_glyph = 'triangle'
track_gene.feature_glyph = 'directed_spliced'
track_transcript.feature_glyph = 'spliced'
track_transcript.feature_colour = [0,0.5,0]

# Add data to tracks
DATA.each do |line|
 line.chomp!
 ref, type, name, location, link = line.split(/\s+/)
 if link == ''
   link = nil
 end
 if type == 'SNP'
   track_SNP.add_feature(name, location, link)
 elsif type == 'gene'
   track_gene.add_feature(name, location, link)
 elsif type == 'transcript'
   track_transcript.add_feature(name, location, link)
 end
end

# And draw
my_panel.draw('my_panel.png')

__END__
chr1  gene        CYP2D6      complement(80..120)
chr1  gene        ALDH        100..449
chr1  SNP         rs1234      107
chr1  gene        bla         complement(400..430)
chr1  SNP         rs9876      44
chr1  gene        some_gene   complement(join(170..231,264..299,350..360,409..445))
chr1  transcript  transcript1 join(250..300,390..425)
chr1  transcript  transcript2 253..330
chr1  transcript  transcript3 266..344
chr1  transcript  transcript4 complement(join(410..430,239..286,129..151))

If someone would actually be interested in getting the library behind this, just let me know. It should be really easy to incorporate this in a rails application where the data are actually stored in a database.

I wonder what if any role _why's Shoes thing would/could play...

UPDATE: This library has now been improved a bit and is hosted on rubyforge. You can find a tutorial and the whole API documentation at http://bio-graphics.rubyforge.org. You can find instructions on how to install and use it over there.

UPDATE TWO: Forget the previous update. I have moved the bio-graphics code to github. See http://github.com/jandot/bio-graphics. That should make it much easier to fork the code and get more input from other developers.

ActiveRecord - all vs all relationships

2007-09-04T15:08:00.000+02:00

Modeling genetics or genomics data presents its own challenges. One of the issues is that the actual definition of things change over time. A database system can only be based on the scientific knowledge at the time of conception. The prime example of course is the definition of a gene over the years. Before 1997, it was believed that the vast majority of these encoded proteins. As a result, 'genes' tables in databases typically had columns to store information on the start and stop codon. However, it became clear that many genes actually do not encode proteins, forcing the remodeling of biological databases. But that's not the topic of this post.

What is the topic here, is how relationships can be stored in a database. Suppose I want to store mapping data: markers mapped to linkage groups, clones mapped to physical maps, ... Markers are stored in a markers table, clones are stored in a clones table, linkage groups in a linkage_groups table; you get the point.

The database that I'm working with at the moment (and only have read-access to), stores the mappings in a mappings table which includes the following columns:

map_type
map_name
mapped_object_type
mapped_object_name

So records could look like:

 map_type       map_id  map_name      mapped_object_type  mapped_object_id  mapped_object_name
--------------+-------+-------------+-------------------+-----------------+------------------
chromosome     1       chromosome_1  marker              1                 marker_A
chromosome     1       chromosome_1  marker              2                 marker_B
physical_map   2       ctg1          clone               1                 clone_A
physical_map   3       ctg2          clone               2                 clone_B

To make things worse, markers can also be mapped to clones. This means that any clone can act as a marker, but also as a map at the same time.

 map_type       map_id  map_name      mapped_object_type  mapped_object_id  mapped_object_name
--------------+-------+-------------+-------------------+-----------------+------------------
clone          1       clone_A       marker              1                 marker_A

How can I model this in ActiveRecord? There's the concept of polymorphisms in ruby, which could solve this relationship nightmare if there would be only one thing in the mappings table that's polymorphic. But as it happens, there's two... Evan Weaver wrote this rails plugin has_many_polymorphs, which should do the trick (see here for a tutorial and background if it's unclear what I'm talking about). Unfortunately, as it is focussed on rails and not on ActiveRecord in general, it doesn't handle namespaces.

So here's what I've come up with:


module MyNameSpace
  class Mapping < ActiveRecord::Base
    # Relationships to feature-like things
    belongs_to :marker, :foreign_key => 'mapped_object_id', :conditions => ["mapped_object_type = 'marker'"]
    belongs_to :clone, :foreign_key => 'mapped_object_id', :conditions => ["mapped_object_type = 'clone'"]

    # Relationships to map-like things
    belongs_to :chromosome, :foreign_key => 'map_id', :conditions => ["map_type = 'chromosome'"]
    belongs_to :physical_map, :foreign_key => 'map_id', :conditions => ["map_type = 'physical_map'"]
    belongs_to :clone, :foreign_key => 'map_id', :conditions => ["map_type = 'clone'"]
  end

  class Marker < ActiveRecord::Base
    has_many :mappings_as_feature, :class_name => 'Mapping', :foreign_key => 'mapped_object_id', :conditions => "mapped_object_type = 'marker'"
    has_many :chromosomes, :through => :mappings_as_feature
    has_many :clones, :through => :mappings_as_feature
  end

  class Chromosome < ActiveRecord::Base
    has_many :mappings_as_map, :class_name => 'Mapping', :foreign_key => 'map_id', :conditions => "map_type = 'chromosome'"
    has_many :markers, :through => :mappings_as_map
  end

  class PhysicalMap < ActiveRecord::Base
    has_many :mappings_as_map, :class_name => 'Mapping', :foreign_key => 'map_id', :conditions => "map_type = 'physical_map'"
    has_many :clones, :through => :mappings_as_map
  end

  class Clone < ActiveRecord::Base
    # Relationships where the clone is the feature
    has_many :mappings_as_feature, :class_name => 'Mapping', :foreign_key => 'mapped_object_id', :conditions => "mapped_object_type = 'clone'"
    has_many :physical_maps, :through => :mappings_as_feature

    # Relationships where the clone is the map
    has_many :mappings_as_map, :class_name => 'Mapping', :foreign_key => 'map_id', :conditions => "map_type = 'clone'"
    has_many :markers, :through => :mappings_as_map
  end
end

The key here is to make the distinguish between mappings_as_feature and mappings_as_map. A marker object can only have mappings where it acts as a feature, while a clone can both have mappings where it acts as a feature and where it acts as a map.

Using this code, it's now possible to do:


clone = Clone.find_by_name('clone_A')
puts clone.mappings_as_map.to_yaml
puts clone.mappings_as_feature.to_yaml
puts clone.markers.to_yaml
puts clone.physical_maps.to_yaml

Voila (until further notice...).

UPDATE: Pratik blogged about has_many_polymorphs and lists the generated associations here.

A ruby API to the Ensembl database

2007-08-14T14:36:00.001+02:00

"Joy to the world, lalaa la laaaa." I can finally announce that I've released the ruby API to the Ensembl core database under the bioruby-annex umbrella. Go here for the release.

What is it?
The Ensembl database stores genetic and genomic data on a variety of species: sequences of chromosomes and positions of features such as genes and polymorphisms. This data is browseable using their genome browser, but is also directly accessible if you connect to their mysql database. A perl API to that database has been available from the start and is used by the ensembl people themselves to handle the data. A java implementation (called Ensj) is also available, but I don't know the status of that one. The ruby version should provide similar functionality to the perl API, albeit for querying only and not for writing to the database.

This API is aimed at the core database. Ensembl also provides the variation and compara databases, but these are not the focus of the current API implementation.

A minimal interface to the data of Ensembl was already available through Mitsuteru Nakao's ensembl.rb library in the bioruby project, and is based on the exportview functionality of Ensembl's web interface. Although very useful, it does not give the full functionality that can be achieved by accessing the database directly.

The ruby API basically provides two things: access to the data in the database, and transformations of those data.
Access to the data. (Virtually) all tables of the database are available through ActiveRecord, with all the automated query methods associated with that ('find_by_anything_you_like'). Say you want to get the object of a transcript with stable_id "ENST00000380593", you'd do

transcript = Ensembl::Core::Transcript.find_by_stable_id('ENST00000380593')

Transformations of the data. You might have the coordinates of a gene on the chromosome, but actually want them on a contig or supercontig. This is where the Sliceable#transform and Slice#project methods come in. In contrast to the perl API, there is no Sliceable#transfer method, because my interpretation of a 'slice' is slightly different from the perl implementation. Read the tutorial for more information.

Minimal script
Any script using the API would have to these steps:

require the library
include the Ensembl::Core namespace (not strictly necessary, but saves typing)
connect to the database
start doing stuff

So for example:

#!/usr/bin/ruby
require 'rubygems'
require_gem 'ensembl-api'

include Ensembl::Core

CoreDBConnection.connect('homo_sapiens')

transcript = Transcript.find_by_stable_id('ENST00000380593')
puts "5'UTR: " + transcript.five_prime_utr_seq

How to install
The API has been released as a gem file, which you can either download from the website and install using the command

gem install ensembl-api-0.9.gem

, or export from the SubVersion repository using the command

svn export svn://rubyforge.org/var/svn/bioruby-annex/ensembl-api

This gem depends on bioruby and ActiveRecord.

UPDATE: The code has been moved from rubyforge to github. Get it from http://github.com/jandot/ruby-ensembl-api

Documentation
Check the website at rubyforge, which will show the tutorial (based on the perl version) and the rdoc documentation. In addition, there are the tests in your gem directory, plus a sample script that shows all functionality of the perl-version of the API called examples_perl_tutorial.rb.

Credits
I owe a lot to the Ensembl core team for helping me out when I was at the Ensembl site as a "Geek for a Week"...

Call for help
If anyone would be interested in improving the API, don't hesitate to contact me. At the moment, for example, projections between coordinate systems only work if they're directly linked in the assembly table, and projections of the haplotype assembly_exceptions will now raise a NotImplementedError error. In addition, it would be very useful if we could add the variation and compara databases to the API.

ActiveRecord and mysql: show my databases

2007-07-30T14:05:00.000+02:00

Working on a ruby API for the Ensembl databases, I bumped into the issue of having to connect to a database without knowing its name.

The ensembl database server hosts databases for each species. Every two months or so, there's a new release which means a new database for every single species. To see what databases are there, you can login to the ensembl server with mysql:


mysql -h ensembldb.ensembl.org -u anonymous

(for more information, see here).

The command "show databases;" on the mysql prompt lists a total of 1035 databases at the moment, a short selection looks like this:


bos_taurus_core_41_2d
bos_taurus_core_42_2e
bos_taurus_core_43_3
bos_taurus_core_44_3a
bos_taurus_core_45_3b
bos_taurus_est_36_2
bos_taurus_est_37_2a
homo_sapiens_core_36_35i
homo_sapiens_core_37_35j
homo_sapiens_core_38_36
homo_sapiens_core_39_36a
homo_sapiens_core_40_36b
homo_sapiens_core_41_36c
homo_sapiens_core_42_36d
homo_sapiens_core_43_36e
homo_sapiens_core_44_36f
homo_sapiens_core_45_36g
homo_sapiens_core_expression_est_34_35g
homo_sapiens_core_expression_est_45_36g
homo_sapiens_core_expression_gnf_34_35g
homo_sapiens_core_expression_gnf_45_36g

To connect to the homo_sapiens_core_45_36g database, type "use homo_sapiens_core_45_36g;" at the mysql prompt. However, as all 'core' databases have the same database schema, the API applies to all of these species, and just has to connect to different databases. But how do you go about doing that? What you could do, is provide the full database name in the establish_connection statement. But having to memorize these full names, or having to open mysql connections prior to writing your scripts is, to say the least, far from optimal. But how do you query a database system without connecting to a particular database?

Basically, you make a connection to the host without specifying a database, and send the raw sql query "show databases;" over that connection. The code below does just that.


ENSEMBL_RELEASE = 45
DB_ADAPTER = 'mysql'
DB_HOST = 'ensembldb.ensembl.org'
DB_USERNAME = 'anonymous'
DB_PASSWORD = ''

class DummyDBConnection < ActiveRecord::Base
  self.abstract_class = true

  establish_connection(
                      :adapter => DB_ADAPTER,
                      :host => DB_HOST,
                      :database => '',
                      :username => DB_USERNAME,
                      :password => DB_PASSWORD
                    )
end

class CoreDBConnection < ActiveRecord::Base
  self.abstract_class = true

  def self.connect(species)
    db_name = DummyDBConnection.connection.select_values('show databases').select{|v| v =~ /#{species}_core_#{ENSEMBL_RELEASE.to_s}/}[0]


    if db_name.nil?
      warn "WARNING: No connection to database established. Check that the species is in snake_case (was: #{species})."
    else
      establish_connection(
                          :adapter => DB_ADAPTER,
                          :host => DB_HOST,
                          :database => db_name,
                          :username => DB_USERNAME,
                          :password => DB_PASSWORD
                        )
    end
  end
end

And then just have your classes (e.g. CoordSystem, SeqRegion, Gene) inherit from CoreDBConnection instead of ActiveRecord::Base.

To make the actual connection, start your script with:


CoreDBConnection.connect('bos_taurus')

I'm currently at Ensembl for a week ("Geek for a Week") to work on the full-blown ruby API, and am planning to give an introduction on how to use it in one of the later posts.

How do you process literature?

2007-07-25T21:23:00.000+02:00

A quick glance at the side of my desk reveals two stacks of manuscripts to read; each stack about 20cm high. Sounds familiar? There seems to be a major task in front of me to process all that.
First thing to do is to identify what caused those piles in the first place. The answer: no system that I'm satisfied with for reference management. Of course, there is software like Reference Manager and EndNote as well as websites like Connotea and CiteULike. But they all have one major flaw: they are not suited to store the knowledge gained from those papers. Entering a reference to those papers in the software is not the same as going through them and extracting useful information. Sure, they do have a notes field where you can jot down some short remarks, but often knowledge is much easier recorded and remembered in little graphs and pictures than in words. There's reference management, and there's knowledge management.

What do I want my system to look like? First of all, it should be searchable. The tagging system provided by CiteULike/Connotea seems good for that. Also (and this might seem illogical for a bioinformatician), the system should not be fully automatic or even electronic, but analog. Why? Just pressing a button to for example add the abstract from a paper to the system gives a sense of... what's the word in English: volatility? For some things you should use the help of a computer, and for some you shouldn't. There's a difference between using Excel to repeat the same calculation 50 times, and trying to use a pc for storing knowledge. It's me who needs to store that knowledge, not the computer. If that was the case, I could always go back to Google instead of making the effort of using a reference manager in the first place. I've played around with zotero and personal wikis in the past, and they just didn't do the trick: I still ended up just copy-pasting the information instead of absorbing it.

Another advantage of using an analog system, is that when you feel your productivity behind your computer is suboptimal, you can always take your cards, find yourself a quiet place, put your feet on a desk, and flick through the things you wrote down. Slippers and pipe are optional.

During my PhD a few years ago, I used a system that was exclusively based on index cards. The inspiration came from Umberto Eco's book "Come si fa una tesi di laurea" (Or "How to make a doctoral thesis") (1977) in which he explains how he handles the knowledge for his book research. For each manuscript, I'd make a new card. The front contained an identifier, the paper title, full reference and keywords. On the back I'd write down what I had to remember from that paper, including little graphs, schemas and stuff. I've got to admit that a drawback of using these cards was that they were not easily searchable, but linking them worked quite well with a bit of discipline.
During those years, I used the index card system both as reference manager and as knowledgebase. Although it did work to satisfaction, the role of reference manager should be fulfilled by a better tool.

Now how could I implement something like that into a workflow? Basically, any new paper to be read should be entered in CiteULike and tagged as 'to_read'. When I've got time to read it: see if it's necessary to print out, or preferably read from the screen (we want to be nice to the trees, don't we?). When I've read the manuscript and there is interesting information to remember, only then create an index card. In case it's a landmark paper and/or I've been adding a lot of comments and markings in the text: keep the printout as well, and mark the index card that I've got that printout as well.
Let's try this out for a few weeks and see where it goes...

BTW: for a knowledgebase system based on index cards taken to the extreme, see PoIC (Pile of Index Cards).

Documenting one-off parsers

2007-07-18T12:39:00.000+02:00

A lot of day-to-day work consists of parsing data files to transform the contents from one format into another or to create statistics. However, when you have to get back to those scripts at a later stage - you need something similar in another project or you notice that something along the way must have gone horribly wrong - it can often be quite hard to figure out what the script actually did. Having meaningless script filenames like ParseBlast.rb doesn't help either. (Parse BLAST into what?)

I must say that things improved a lot when I switched from perl to ruby. Trying to understand a perl script that I wrote a couple of weeks earlier was a real pain, while I normally have no problems understanding ruby code that I wrote months ago... Has a lot to do with the simple syntax and expressiveness (is that a word?) of that language.

But just understanding the code is not enough. To be able to assess if you could use a script from an earlier project in a new one, you often also have to get hold of information in that script that is not inherently encoded in its code: what was the project? What did the input look like? What did the output have to look like? That's where the script documentation comes in, because having that information easily available greatly reduces the time you need to assess if you should copy-paste-adapt the script, or just start from scratch.

What I try to do, is always use the same template when starting a new parsing script. Even when the parsing itself would only take 5 lines, I try to use this whole template. A quick walk-through:

Lines 2-32: The documentation, consisting of:

script name
short usage message
a description of what the script does
list of arguments that can be used
input format
output format
contact information

Lines 38-40: Definitions of classes that are used later on in the script itself
Lines 42-51: Parsing the arguments to the script. If the user (i.e. me when I run it) uses --help or an argument that does not exist, he automatically gets the documentation of the script. If I would use the -a tag here, the output to the screen would be


USAGE
-----
./this_script.rb [ -h | --help ]
                [ -i | --infile | < ] your_input.txt                  [ -o | --outfile | > your_output.txt ]

Lines 53-58: Logging of the options that were used in running the script.
Lines 60-70: Create the input and output streams.
Lines 72-75: The actual parsing code. This is the bit that does the work (using the classes described in lines 38-40).
Lines 77-79: Clean up.

Here's the complete template.


#!/usr/bin/ruby
# == NAME
# this_script.rb
#
# == USAGE
#  ./this_script.rb [ -h | --help ]
#                   [ -i | --infile | < ] your_input.txt #                   [ -o | --outfile | > your_output.txt ]
#
# == DESCRIPTION
# Description of what this script does...
#
# == OPTIONS
#  -h,--help::                 Show help
#  -i,--infile=INFILE::        Name of input file. STDIN if not defined.
#  -o,--outfile=OUTFILE::      Name of output file. STDOUT if not defined.
#
# == FORMAT INPUT
#   >gi|4531835|bla
#   ACTTACCGACCGACTGACTACTTATGCCA
#   >gi|4861534|blabla
#   CTACCCCATCTACCGGGGCTCGACT
#   ...
#
# == FORMAT OUTPUT
#   4531835    29
#   4861534    25
#   ...
#
# == AUTHOR
#   my full contact information

require 'rdoc/usage'
require 'optparse'
require 'ostruct'
require 'logger'

### Define classes here
class MyClass
end

### Get the script arguments and open relevant files
options = OpenStruct.new()
opts = OptionParser.new()
opts.on("-h","--help",
       "Display the usage information") {RDoc::usage}
opts.on("-i","--infile", "=INFILE",
       "Input file name") {|argument| options.infile = argument}
opts.on("-o","--outfile", "=OUTFILE",
       "Output file name") {|argument| options.outfile = argument}
opts.parse! rescue RDoc::usage('usage')

log = Logger.new(File.new('this_script.log', File::WRONLY | File::TRUNC | File::CREAT))
log_level = Logger::INFO # or: DEBUG, WARN, FATAL, UNKNOWN

log.info('Script this_script.rb started')
log.info('Options:')
log.info(options.to_yaml)

if options.infile
 input_stream = File.open(options.infile)
else
 input_stream = $stdin
end

if options.outfile
 output_stream = File.new(options.outfile,'w')
else
 output_stream = $stdout
end

### Actually do some stuff
input_stream.each_line do |line|
 output_stream.puts line
end

### Wrap everything up
output_stream.close
input_stream.close

(If I wanted this a bit DRYer, I'd have a constant with the script name at the top so I wouldn't have to repeat that filename over and over. To be done...)

If I wanted to over-organize, I'd create a little database with the descriptions of those scripts so I can search through them. Even though that would have it's use from time to time, that would be taking it too far. Having that documentation and a standard way of providing the command line arguments is enough for me for now.

I must admit I let this slip in the last couple of months, which doesn't mean it didn't work the months before that. That's just what happens when you go on holiday and completely forgot the habit of doing this afterwards.

UPDATE: InfiniteRed blogged about a similar approach later.

A choice of databases, or PostGres vs SQLite

2007-07-11T15:26:00.000+02:00

When managing data in bioinformatics, one of the main tools you use is databases to store the stuff. In many cases, flatfiles are sufficient, but sometimes you need the flexibility and reliance of a nice relational database. (Note: Excel is not a database). Actually, in quite a few cases it does make sense to use flat text files instead of databases, but that's the subject of another post.

As there is a plethora of RDBMS (relational database management systems) available, the first thing you have to ask yourself is what system to use. There's mysql, oracle and others. But for the purpose of this post, I'll focus on PostGreSQL (http://www.postgresql.org) and SQLite3 (http://www.sqlite.org), because that's what I normally use.

For the big projects with a lot of data, I normally use postgres, while smaller projects are typically served by sqlite databases. Both of these however have their good and bad things. Choice depends on a number of factors:

Installation. The main difference between postgres and sqlite3 is that one is based on running a server and the other is not. To be able to use postgres, you'll have to install the software, create a postgres user on your system and have that user start the postgres server. Next thing to do is create a user within the database system and grant it the right permissions (e.g. to create other users or databases). That's a bit more complex than installing sqlite3.

Database creation. To create a new database in PostGreSQL, you issue a createdb MyDatabase. This will add the database to your RDBMS. You can get into the database with psql MyDatabase. In contrast, sqlite3 databases are just text-files like any other. Create a new databasebase by calling it: sqlite3 my_database.s3db. This will create a file in your present directory called my_database.s3db. Note: you don't have to use the .s3db extension, but that makes it recognizable for the SQLiteAdmin tool (see below).

Use. Once you've got your databases created, their use is really similar. On the command line, use the psql or sqlite3 command. On Windows, you can e.g. pgAdmin and SQLiteAdmin.

I must say that psql is a bit easier on day-to-day typing, because it's got tab-completion, which I haven't been able to activate yet in sqlite3. I should spend some time one day to configure sqlite exactly like I want it to.

Locking. This is what bit me in the ankles with sqlite3... According to the website: "SQLite version 2.8 allowed multiple simultaneous readers or a single writer but not both. SQLite version 3.0 allows one process to begin writing the database while other processes continue to read. The writer must still obtain an exclusive lock on the database for a brief interval in order to commit its changes, but the exclusive lock is no longer required for the entire write operation." Indeed, I could load hundreds/thousands of records while simultaneously querying the database. But at what looked like a random moment, after hours of loading, it would choke, leaving me with a partially populated database. Any small child can tell you that this is not what you want to happen. I've never encountered something like that in postgres, and probably never will.

Backup. Because sqlite if file-based, making a backup of a database is nothing more than taking a copy and moving that copy to another hard disk or burning it on a CD-ROM. To make backups of postgres databases, we typically use pg_dump, which dumps the contents (and if you want the schema as well) of a database in a big-ass text file. So the backup is outside of the RDBMS itself.

So both RDBMS have their advantages and disadvantages. For huge projects, I tend to use PostGres, although I've started to use sqlite more and more lately. If only it wouldn't choke on the locking...

Bioinformatics and labbooks

2007-07-05T22:11:00.001+02:00

What I've really struggled with in the first years after I ended up in bioinformatics (still in denial back then), was how to record what I'd done. Researchers in a wet-lab environment typically use a labbook in which they write down what protocols they used for an experiment, together with the results and nice pictures of their PCR electrophoresis gels. This is a bit more difficult for bioinformaticians. Why? I think because the physical actions that are required to do the research are less pronounced. Filling out a PCR plate, putting it into a PCR machine and loading the result on a gel is a bit more effort than pressing Enter. When we code, what we actually do, is create the PCR machine hardware. The programming language acts as the nuts and bolts used to create that machine. The effort for the bioinformatician is in building the hardware, not running the tests themselves (that's the sweat and tears of the servers, not ours).

So I was pleasantly surprised when I arrived at my present job 2-and-a-bit years ago when they told me they even had an SOP on how to perform and record bioinformatics work. Now SOPs can be restrictive and more about paperwork than actually helping you out, but coming from an environment where everything (or rather: nothing) goes regarding logging bioinformatics work, this piece of paper was a big relief.

The main principle of a bioinformatics SOP is that you have to be able to record all steps you followed to transform one piece of information into another. Same as what a lab-oriented SOP does: "how do you get from a BAC library and a pair of PCR primers to knowing which BAC is positive for that marker". But of course there are other things to record than in a lab. The main principle is that we have to be able to prove that the data we use are what we say they are. Same goes for the scripts.

The input data. The startpoint for most of my dry-lab experiments is data that I downloaded from the internet. So what do I (have to) record? The URL of the file and the md5sum of the file once it's downloaded. That md5sum makes sure that, in case I have to reanalyze, I can check if that big sequencing centre has changed the contents of their files or not (of course without telling us).

The scripts. Same goes for the scripts. Apart from the fact that it's always a good thing to document your scripts, you also have to be able to prove afterwards that that actually was the script you used, and not just some other piece of code with the same name. Can be done with md5sums again.
It does happen that I logged the md5sum and still change the script afterwards (e.g. setting constants to another value). As long as that is logged as well, we can reconstruct the original script. Of course it is better to circumvent that problem by creating generic code, but sometimes the effort of doing just that outweighs the benefits gained. (Might be a subject of a later post.) Having the md5sum of the script acts as logging the manufacturer and type of PCR machine you used.

The output. The combination of identity-checked input and scripts should always produce the same output. For completeness sake, the SOP tells us to also log the md5sum again...

The system around it. So how do we actually do the recording itself? Where I work, the end product has to be a filled-in labbook. What we do, is use a task tracker (RT Request Tracker) to, well, track our tasks (duh). The moment we start a new project (let's define that here the GTD-way as anything that consists of more than one physical action), we create a new ticket (e.g. "Identify syntenic region between species A and B") and log everything in there: the project directory, background information, md5sums, workflows, interpretation of the results. When all is finished, we make a hard copy (well: print it out) and glue it into our labbooks.

In some settings it might be more sensible to log things in a wiki, as explained by Mike at bioinformaticszen (here), where he talks about using a hyperlinked document or a wiki to track what you've done. Of course it can make sense in many environments (i.e. if you don't have to care about audits and are tracking stuff merely for yourself), but the moment you have to be able to present stuff for audits or patent applications, it's critical that the documentation you generated is immutable. Don't get me wrong: I'm a big fan of the wiki-way (and use one at work as well), but not as a labbook. What do I find a wiki useful for, is operating guidelines, which are allowed to change over time.

Six months of Getting Things Done

2007-07-03T14:43:00.000+02:00

It's the start of July and about 6 months after I started implementing the Getting Things Done meme. And, man, did it change things... In this post, I'll share some of my experiences and the obstacles I encountered.

What is that Getting Things Done already?

I'm not going to explain the whole thing here. There are better resources for that (for a list of URLs, see the bottom of this post). Most importantly: get the book "Getting Things Done: The Art of Stress-Free Productivity" by David Allen. Basically, GTD is a way of turning the overwhelming amount of 'stuff' we have to remember to do into a manageable system and out of your head.

Implementing the GTD system has helped me a lot in the last half year. For starters, I got things done (more than I would normally have, I think) while at the same time it allowed me to have less on my mind. I no longer have that dreaded nagging feeling of "I should remember to do this" while knowing that I would forget something else by remembering it. The time that I couldn't get to sleep because I had to remember this and this and that and the other is now a distant memory. One of the effects that it also had was that it was easier to quickly switch between different jobs I had to do. I'm not yet at the stage of what the book calls a "mind like water" (i.e. neither over- or underreacting to anything), but it's started to get a bit fluid. Stuff like responding inappropriately to email, projects, thoughts about what I need to do (the over- or underreacting) leads to less effective results.

* What works / my system

What parts of the GTD system work for me? Up till now, it looks like I mainly focussed on the day-to-day tasks rather than the long-term goals and someday/maybe.

The hardware:

A good pen.

A filofax. After trying different low and high-tech ways to implement GTD, I finally ended up using a filofax. It's a black pocket-size (8x12cm) one called Identity. I renamed the different tabs in it to 'Next Actions', 'Projects', 'Waiting For', 'Someday', 'Calendar' and 'Lists'.

Project folders. I started creating project folders both at work and at home. I used to take the same notebook to all my meetings to take my minutes. That's now shifted to taking a single piece of paper to them, taking my notes and then archiving that paper in the project folder. It's always good to just be able to take a single project folder to a meeting containing all minutes and other supporting information.

How I use that hardware:

Next Actions and Projects. Always having the filofax with me, I instantly write down anything that comes into my head that I have to remember to do. At that moment, I don't make the decision if I have to put it in next actions or projects or someday/maybe. I just jot it down on my next actions list. The moment that I come to actually doing that thing, it might be sensible to convert it into a project or a someday/maybe item. At any one moment, I have 40 or so items on my next actions list.
David Allen talks about having different next action lists for different contexts: a list for stuff to do at home (@home), a list for stuff to do at the office (@office), while at a computer (@pc), while with a phone (@phone) or when going shopping (@shop). That definitely didn't work for me. I ended up having just one next actions list, using the type of bullet to distinguish between the different contexts. The different bullet point types I use are: a triangle (i.e. a stylized tear-drop for the tears of sweat I shed at work), a square with its bottom line missing (a stylized version of a house), a little phone and a dollar-symbol.

I noticed that it's really important to make a real distinction in how you phrase the next actions and projects. Projects are basically things you have to do that will require more than 1 step, while next actions are the indivisible counterpart. It really helps if you phrase the next actions as verbs ("write unit test"), while the projects are phrased as end-points rather than verbs ("API published" instead of "create API").

Calendar. I've been relatively successful in putting all my meetings in the filofax calendar. To remember to do stuff with a given deadline, I write a note in my calendar for example a forthnight before that says "add to next actions: write poster abstract". I do something similar for recurring events: I'll fill the next ten or so, and add a note with the last occurrence to fill in the next ten.

Templates. The filofax came with a bunch of empty templates to use. As it didn't take too long to get through those, I created a new template in OpenOffice with a grid and showing where to punch the holes. Just printing it out and making double-sided copies gives me all the empty sheets I want.

* What kinda works

The weekly review. This is one of the cornerstones of the system, where you take a little time once a week to go over your next actions, projects and other lists. Until recently, I've neglected this quite a few times. However, it's easy to pick up again and that's what I did. What I do, is purge the next actions list (taking the time to rewrite the actions on virgin sheets), check that I have at least one next action related to each project, and ask myself if I have to chase people on my waiting_for list.

Empty inbox. One of the things you _can_ do, is keep an empty inbox. I got to the stage of the empty inbox during a few weeks, but let it slip again. At the moment there are 1060 mails in it, but I know it will be relatively straightforward to do the big reorganize/purge exercise again.

* What doesn't work (yet)

Tickler file. The tickler file is a set of folders to organize paperwork that has to be done by a certain date (see the Wikipedia page for GTD). Although I won't implement this at work, we haven't really decided yet if this is really useful or not at home.

GTD software. Even though I'm a bioinformatician, I noticed that using software to keep track of my lists doesn't do it for me. I've tried ThinkingRock, BackPack (from 37signals) as well as a host of other little applications. In the end, not being able to use that system while I'm on the bus, in a shop or just sitting in the living room made it clear that I had to go for the analog version.

Literature. This is what I'm annoyed about the most at the moment. Even though it's not really part of GTD, reading literature should fit into the bigger system in some way or another. I haven't found out how yet. I like to print out the papers and make notes directly on them, but that's no way of organizing the information contained in them. I'm thinking about using CiteULike as a reference manager (don't like Reference Manager or Endnote) and making notes using Adobe Acrobat on the PDF. To be continued...

* The future

I noticed that the elements of GTD that I'm using are focused on day-to-day work. I'm not really happy with how I use the someday/maybe list in my filofax, do not have a list of big goals, ... So while I will continue to use this system in the next couple of months, I'll probably try working on the big picture as well. I might keep a couple of lists at home to do that. It's not necessary to carry the list with things you want to do one day or your major goals with you all the time.

Conclusion

Based on what I've experienced, it is really worthwhile to implement the Getting Things Done strategy to get to stress-free productivity. Start with reading the book. It's important to find your own implementation and decide what things you want to implement or not: next_action lists, empty inbox. You don't have to bring every suggestion into practice. There are many things in the book that I decided I wouldn't need.

It's really easy to start following the GTD philosophy (although becoming a black-belt requires a lot of work), and it's no problem if you fall off the wagon for some reason after a while without wanting to: it's straigthforward to jump on it again. Just start using those lists again.

Some really good tips on starting with GTD can be found on this blog entry.

And now back to some real work...

Links

* The book: see Amazon "Getting Things Done: The Art of Stress-Free Productivity"
* Introduction of concepts
* How to start
* Big list of GTD software

Naming conventions

2007-06-29T16:41:00.001+02:00

Naming conventions. You bump into them every single minute of the day. Naming new directories in your project folder, naming new tables in your database, ... Recently, the issue of naming convention came more to the foreground for me as I'm trying to write a ruby API to one of our databases (see later).

Two of the most-often-encountered naming schemes are CamelCase (ThisIsACamelCaseString) and snake_case (this_is_a_snake_case_string). And in the case of CamelCase: do you make the very first letter a capital or not? If I'm not mistaken, variables in java are often CamelCase, except the first letter (thisCouldBeAJavaVariable).

When thinking of names for directories and files (read also "Organizing yourself as a dry-lab scientist" on BioinformaticsZen), i.e. when there's no set naming convention that you have to follow (e.g. variable naming conventions), I tend to use different schemes for directories versus files for some reason. For naming directories, I use an underscore to separate different things in the same name. For example, I name my folders by concatenating the date they were created with the RT Task Tracker ticket and a description (the latter being in camelcase). For example ~/20070629_RT12345_ThisIsADirectory. For files, I normally use snake_case, except for scripts... Why? Maybe to distinguish those scripts from the data files. For example:


+- Documents
+- Projects
     +- 20070629_RT12345_ThisIsASampleProject
     |       +- input_file.txt
     |       +- output_file.txt
     |       +- log_file.txt
     |       +- ParseInput.rb
     +- 20070629_RT23446_AnotherProject
            +- input_file.txt
            +- output_file.txt
            +- log_file.txt
            +- ParseInput.rb

It would probably not be a bad idea to start to use a default data folder or something to keep all the input, output and other files, and a script folder for the scripts... I should take a look again at the BioinformaticsZen blog.

Of course things are completely different when you're coding or setting up databases. In these cases, your preferred programming language will have it's own conventions. In ruby, for example, classes are in CamelCase, while variables should be snake_case. All good and well, until they start to bite you in the you-know. I'm trying to create a ruby API to an existing database that would require polymorphic associations. This requires that the values in the something_something_type column should be class names, which are CamelCase and singular. But of course, the database has everything in snake_case. I found a workaround to get the thing up and running with snake_case, except that it requires the value to be plural, which it of course isn't. End result: I'll have to get the actual data values in the database changed to get this thing working.

Manual genome annotation tools

2007-06-26T12:04:00.001+02:00

An important part of genomics and genetics research is to know where your genes of interest lie on the genome and what the gene model looks like. In other words: to know where do the exons start and stop, what the UTR boundaries are, and where there are any polymorphisms in those genes. That's called genome annotation, that is.

With the sequencing of any new genome, the annotation of its genes is the logical next step. This often consists of 2 main phases: there's the automated annotation by BLASTing against known sequences or even de novo gene annotation. The second phase is the manual curation of those automated annotations: biologists or bioinformaticians with a biology background looking at those gene models and correcting things like "there should be an additional exon here", "this exon-intron boundary is 2 bp off" or "this actually is an alternative transcript".

Over the last two years, I've tried out several software packages to deal with the manual curation or annotation of sequences. And I must say: it hasn't been great. Let's walk through them:

VectorNTI:
Is a commercial package from Invitrogen. I didn't try this tool recently, partly because my experiences with it on my last job was less then impressive. From what I remember, it becomes unusable when you have to handle larger sequences or when you've got a high density of features. Things might have changed in the mean time, but remarks at the water cooler from colleagues do not support that hope.

Artemis:
The Artemis tool from the Sanger Institute in itself is a great tool with a lot of functionality. You can launch the thing using Java Webstart without having to install it on your own computer. However, I found it to be lacking in user intuitiveness and usability. In addition, it looked like it used non-standard ways to store my annotations. I found that any annotations that I made were stored in the original FASTA-file that I loaded as GenBank annotations at the top. Result: the FASTA-file itself became invalid, and it wasn't a GenBank file either. Still, this tool can be very useful for very small projects.

Apollo:
This is the biggy, developed by the FlyBase people in collaboration with the Sanger Institute in the UK. It's the tool that is referred to the most in the community. But not the tool that I'll use anymore. First of all, it is recommended that you have at least 2 Gigs of RAM when you try to run it. That's right: 2,000 Mb of the stuff. Not your average desktop PC, then... In addition, many people have reported that it crashed on them at random leaving their unsaved work, well, unsaved. The tool also has a lot of features. Too many, actually. Finding out how you can do something can take really long because you've got a haystack of things to go through.

What was the biggest drawback for me, was that I was not able to import my own externally generated results into the tool. You can only do that by creating a GameXML file, which is way to cumbersome to do.

Otterlace:
The HAVANA annotation group uses a tool based on AceDB for annotation. I can just be brief about this: I think it's a really good tool, but as it's not available to annotators outside the HAVANA group, has to be discarded as an option.

Argo:
And then I found Argo (from the Broad Institute). This finally looks like a tool that does what it should do: it has an intuitive interface and overview of your genomic region, it has import and export filters to GFF 1, 2 and 3 as well as GTF. It also allows you to quickly check for non-standard intron-exon boundaries and translations of your gene models and a whole lot of other features. Same as Apollo does, but without making it confusing and complicated. I would suggest this tool in combination with blixem (to check BLAST results in detail).

Databases and ruby (without rails)

2007-06-21T11:23:00.000+02:00

Just bumped into a really nice O'Reilly blog article that combines two of the things I like to work with: GTD and ruby.

As a bioinformatician, I often have to handle quite a lot of data, which I tend to put into databases. Ruby has a fabulous framework in Rails to access and manipulate data, but after having created rails applications for most of those, it became more and more clear that it would be preferable to just have the power of ActiveRecord without having to deploy the whole rails-thing.

Similar to what is discussed in the blog article by Gregory Brown, I created a directory template including the migration code, a connection configuration and a Rakefile. The directory structure looks like this:

+- config
|   +- project_config.yaml
|   +- load_config.yaml
+- db
|   +- migrate
|   |    +- 001_initial_schema.rb
|   +- import
+- lib
|   +- models.rb
+- Rakefile
+- README

Let's walk through this:

The project_config.yaml file contains the project name and connection settings to get to the database. For example

project:
name: MyFunkyProject
database:
adapter: sqlite3
name: db/my_funky_project.s3db

The load_config.rb file uses that information to connect to the database.


require 'rubygems'
require_gem 'activerecord'

class ProjectConfig
  attr_accessor :project_root
  attr_accessor :project_name
  attr_accessor :db_adapter
  attr_accessor :db_name
end

$config = ProjectConfig.new
$config.project_root = File.dirname(__FILE__) + '/..'

YAML.load_documents(File.open($config.project_root + '/config/project_config.yaml')) do |p|
  $config.project_name = p['project']['name']
  $config.db_adapter = p['database']['adapter']
  $config.db_name = p['database']['name']
end

$connection_settings = Hash.new
$connection_settings[:adapter] = $config.db_adapter
if $config.db_adapter == 'sqlite3'
  $connection_settings[:dbfile] = $config.project_root + '/' + $config.db_name
else
  $connection_settings[:database] = $config.db_name
end

ActiveRecord::Base.establish_connection($connection_settings)

The lib/models.rb file contains the... models, obviously. It should 'require' the load_config.rb file to get the connection.


require File.dirname(__FILE__) + '/../config/load_config.rb'

class Task < ActiveRecord::Base
  belongs_to :project
end

class Project < ActiveRecord::Base
  has_many :tasks
end

The db/migrate/001_initial_schema.rb is used to create the database.


require File.dirname(__FILE__) + '/../../config/load_config.rb'

class InitialSchema < ActiveRecord::Migration
  def self.up
    create_table :tasks do |t|
      t.column :description, :string
      t.column :project_id, :integer
    end
    create_table :projects do |t|
      t.column :description, :string
    end
  def self.down
    drop_table :tasks
    drop_table :projects
  end
end

The import directory will then hold a group of loading scripts, that 'require' the lib/models.rb file. They look something like this:


require File.dirname(__FILE__) + '/../../lib/models.rb'

File.open('my_file.txt').each do |line|
  line.chomp!
  # do_something useful
end

As you can see, this directory template - useful as it is for me at the moment - can be rationalized a bit, and I should add tests and stuff. Gregory Brown's article might just give me the right ideas to do that.

Anatomy of data integration

2007-06-19T17:38:00.001+02:00

The paper "Anatomy of data integration" by Brazhnik & Jones (J Biomed Inf 40:252-269 (2007)) gives a clear high-level overview of what is involved in the process of acquiring data from different sources and how to integrate them. Apart from talking about information pipelines and conceptual data models, it delves deeper into the concept of types of data elements (DEs). It really all speaks for itself, but it's nice to be able to name things in a meaningful way. Some good things that you already know, but are still helpful to write down.

Basically, in the context of a data source, you can distinguish between 2 types of DEs: ''focal'' and ''peripheral''. The focal DEs are usually mandatory and have a high quality, often because they are the primary reason for building the database in the first place (e.g. the assay sequence in a dbSNP record). In contrast, the peripheral DEs are often optional and are much more prone to error (e.g. the number of chromosomes sampled in a dbSNP record). As a remark: making peripheral DEs mandatory is asking for trouble. If I don't know the number of chromosomes sampled for a SNP but I'm forced to fill in some number, than that number will insert wrong data which is infinitely worse than having a nil value.

Meaningful integration can only occur between data sources that have a shared pool of focal DEs. You obviously don't try to integrate two different datasets based on the number of chromosomes sampled... This also means that in order to study correlations between two distant domains, you'd need to build what they call a multi-step integration staircase.

From the viewpoint of the integration itself, the focal DEs can be subdivided into two distinct groups: ''integration keys'' and ''informative elements''.
Integration keys are the backbone of the integration and a combination of DEs that identify exactly the same entity in two sources. They are chosen from the overlapping focal DEs. The informative elements represent to goal of the integration and contain the information that we actually want.

What if there are different data sources for the same data set? To make things worse, it might even be that particular data elements are similar but not the same. For example: you can get SNP data from NCBI and from Ensembl, but people have noticed recently that the same SNP can be annotated on a different strand depending on what database you look at. In this case, a practice of keeping all redundant data along with the information about the source becomes important. The quality of the data might vary between the different sources, and having all information available makes it possible for the researcher to make decisions based on all available evidence. "My SNP is on the forward strand according to NCBI, but on the reverse strand according to Ensembl. I trust Ensembl more, so..."