Archive for the ‘Uncategorized’ Category
Heatmaps are a common way to visualize matrices and R provides a variety of methods to generate these diagrams. One of the key features of a heatmap is the color scheme employed. By default the image method uses heat.colors which ranges from red (lowest values) to white (highest values). Other palettes include rainbow and topographical. However, I needed to replicate a color scheme used in a previous publication – dark red for low values and white for high values, passing through shades of pink. While RColorBrewer is an extremely useful package for generating sensible and informative color schemes it didn’t help me out here – the red palette was limited to 9 colors, whereas I needed a larger range. There’s also the colorschemes method from the ClassDiscovery package, but I haven’t checked that out.
So, I ended up making my own color scheme, realizing that shades of pink are a constant red, coupled with increasing amounts of blue and green (in the RGB scheme). Thus generating a matrix of R, G and B values, with a constant red value and increasing blue and green values and then converting the rows to hex colors (via rgb, gives me the required colors.
rgbs <- cbind(255,
seq(32,255, length=256)) / 255
cols <- rgb(rgbs[,1], rgbs[,2], rgbs[,3])
and an example of the color scheme
image(matrix(1:256, nrow=64, byrow=TRUE),
Some handy settings when running a query from the command line via sqlplus
set echo off set heading on set linesize 1024 set pagesize 0 set tab on set trims on set wrap off
-- might want to set column formats here
-- e.g.: column foo format A10
spool stats -- dump results to stats.lst
-- SQL query here
The move to Rockville is complete and I’ve started work at the NCGC. This week has been a bit hectic what with setting up house and getting up to speed at work. But I’m really excited with the stuff that’s going on here – lots of interesting projects in various areas of chemical biology and genomics. Also, being surrounded by some really smart people is invigorating (and leads to some nice lunch discussions). One aspect of the center that I really like is the commitment to Open Source – while there isn’t a whole lot at the moment, there’s some cool stuff coming down the pipeline. And a 30″ Apple Cinema display is sweet.
For now, back to reading up on HTS, RNAi and browsing for sofas.
In a previous post I described how one requires a custom RecordReader class to deal with multi-line records (such as SD files) in a Hadoop program. While it worked fine on a small input file (less than 5MB) I had not addressed the issue of “chunking” and that caused it to fail when dealing with larger files (the code in that post is updated now).
When a Hadoop program is run on an input file, the framework will send chunks of the input file to individual RecordReader instances. Note that it doesn’t actually read the entire file and send around portions of it – that would not scale to petabyte files! Rather, it determines the size of the file and ends start and end offsets into the original file, to the RecordReaders. They then seek to the appropriate position in the original file and then do their work.
The problem with this is that when a RecordReader receives a chunk (defined in terms of start and offsets), it can start in the middle of a record and end in the middle of another record. This shown schematically in the figure, where the input file with 5 multi-line, variable length records is divided into 5 chunks. As you can see, in the general case, chunks don’t start or end on record boundaries.
My initial code, when faced with chunks failed badly since rather than recognizing chunk boundaries it simply read each record in the whole file. Alternatively (and naively) if one simply reads up to a chunk boundary, the last and first records read from that chunk will generally be invalid.
The correct (and simple) strategy for an arbitrary chunk, is to make sure that the start position is not 0. If so, we read the bytes from the start position until we reach the first end of record marker. In general, the record we just read will be incomplete, so we discard it. We then carry on reading complete records as usual. But if, after reading a record, we note that the current file position is beyond the end position of the current chunk, we note that the chunk is done with and just return this last record. Thus, according to the figure, when processing he second chunk from the top, we read in bytes 101 to 120 and discard that data. We then start reading the initial portion of Record 3 until the end of the record, at position 250 – even though we’ve gone beyond the chunk boundary at position 200. However we now flag that we’re done with the chunk and carry on.
When another RecordReader class gets the next chunk starting at position 200, it will be dumped into the middle of Record 3. But, according to our strategy, we simply read till the end of record marker at position 250 and discard the data. This is OK, since the RecordReader instance that handled the previous chunk already read the whole of Record 3.
The two edge cases here are when the chunk starts at position 0 (beginning of the input file) and the chunk ends at the end of file. In the former case, we don’t discard anything, but simply process the entire chunk plus a bit beyond it to get the entire last record for this chunk. For the latter case, we simply check whether we’re at the end of the file and flag it to the nextKeyValue() method.
The implementation of this strategy is shown in the SDFRecordReader class listing.
In hindsight this is pretty obvious, but I was bashing myself for a while and hopefully this explanation saves others some effort.
I’ve been at IU for 3 years now, the last two as a visiting faculty member. However the time has come for a change. I have accepted a position at the NIH Chemical Genomics Center working in the informatics group and will start there in mid-May.
I’m pretty excited about this move, as the group is doing some pretty cool stuff with cheminformatics and high-throughput assays. The group releases their code into the public domain and you can see some of their stuff here. From a scientific point of view, this is an excellent place to do cheminformatics (as well as informatics and computation in general) because the biologists and chemists are right next door. This means access to fresh, raw data as well as the possibility of working with wet chemistry and biology to test out in silico analyses and predictions. (They also have a neat robotics set up, also next door).
As part of the informatics group, I will have fingers in many pies – modeling and computation, developing software, data management and project support. I’m particularly excited about the forthcoming RNAi screening infrastructure that is being set up at the Center.
What happens to the web services?
All services are currently running, but will not be available at the rguha.ath.cx URL’s after mid May. The services have been transferred to David Wild and I will try and redirect URL’s before I leave – if you find that a service you are using is not working, get in touch with David. However, the source code for the SOAP services are freely available on SourceForge, so you can always deploy them into your local Tomcat container. However, I do not plan on supporting these services any longer.
As I have described previously, I am shifting to REST based services. Thus all the CDK based cheminformatics services are converted to REST forms and are currently hosted of rguha.ath.cx, but after I move, this URL will no longer be available. However, running these services is very easy and is described here. If anybody is interested in hosting these REST services, please get in touch with me. In terms of future support, I plan to work on these REST services in my spare time. Note that these services replace a number of my mod_python based REST services (see the GitHub repository which I will updated shortly to remove deprecated services). In general, I plan to maintain the Python REST services. In particular, some of these services are still front-ends to SOAP services. Where possible, I will reomove the dependency on the SOAP service. The prediction service is something I plan to work on as it is a very light weight way to deploy predictive models built using R (and less hassle than going via the SOAP based R infrastructure). Unfortunately, until I work out some server solution, the predictive model service will likely not be publicly accessible (as a service, that is).
Over the last three years I’ve been maintaining a partial mirror of PubChem and the Pub3D database. After my move I will no longer be handling these, so questions should be sent to David. A number of services that I provide make use of these databases. Examples include, 3D structure search, InChI/InChI key resolution, synonyms and so on. These services will continue to run, though the URL’s will likely change. – as before, if you see something is down, David should be able to point you in the right direction. I will likely continue work on the 3D database (especially the similarity search aspects, which is currently flawed), but at this point I don’t know when or how new results will be posted.
I should point out that future work on the projects described here are spare-time projects and not part of my official duties at NCGC. So progress on these is not going to be particularly fast.
Finally, my website has moved to http://rguha.net and the old one should automatically redirect you to the new one.