Showing posts with label scripting. Show all posts
Showing posts with label scripting. Show all posts

Wednesday, November 14, 2007

OS X 10.5 Leopard and MySQL

UPDATE: The Mysql.com supplied preferences pane now works. Just use that.

It seems that the once super-simple MySQL on OS X install is no more in Leopard, that is until the MySQL developer community catches up. Now there is a small amount of terminal-fu that must be done to enable a mysql server on OS X, which has been covered in some other blog posts. But I will repeat the procedure here, since I added a script or two to make things a tiny bit simpler for me.
First, if you did a clean install of Leopard, then you need to download and install the MySQL server for OS X from the mysql.com. Just install the main server, not the preferences pane or the start up item.
Next download this launchd file and copy it to "/Library/LaunchDaemons/com.mysql.mysqld.plist"
Next change the ownership to root. Might as well set proper permissions while you are at it:

sudo chown root:wheel /Library/LaunchDaemons/com.mysql.mysqld.plist
sudo chmod 644 /Library/LaunchDaemons/com.mysql.mysqld.plist


Next load the file into launchd using launchctl

sudo launchctl load /Library/LaunchDaemons/com.mysql.mysqld.plist

Finally, you can create some shell script to load and unload the server at will and place it where you normally place your binaries (mine is ~/bin/sqlserverctl ):
#!/bin/sh
case "$1" in
start)
sudo launchctl load -w /Library/LaunchDaemons/com.mysql.mysqld.plist
;;
stop)
sudo launchctl unload -w /Library/LaunchDaemons/com.mysql.mysqld.plist
;;
*)
echo "USAGE: $0 {start|stop}"
;;
esac

UPDATE: There was a small error in my start script, namely and extra space in the case statement, it has been modified and fixed above.

UPDATE: If you install MySQL via MacPorts, then use the "+server" varient and use this launchd file instead
/Library/LaunchDaemons/org.macports.mysql5.plist

Thursday, September 27, 2007

Mad plotz

A recent submission to a journal has caused us a few headaches over the past few weeks, as the editors sent the paper back to us stating that we did not meet the minimal reporting requirements for the type of experiment that was performed. Which is poetic justice in a way, since for the past few years I have been promoting the use of minimal reporting requirements and standard data formats.

I think this particular journal, however, has gone a bit too far in asking for annotated spectra for every identification in the result set. For most low-throughput experiments this is not such a big deal, but we had thousands of identifications and even wrote an algorithm to automatically assign a quality score to those identifications so that such manual validation of spectra should not be necessary.

But I digress. Instead of fighting it, we decided to give the editors what they want, annotated spectra for every hit. It turns out that this is not such a trivial thing to do. Even gathering all of the data was a tough job, since the experiment was performed many years ago on instrumentation that is nearing its end of life. A lot of file parsing and data reorganization had to be done, prior to any development effort to produce the images that the journal wanted.

A bit of background and some numbers will help us understand the enormous task we undertook. The experiment was a proteomics profile of two developmental stages in zebra fish. We used two methodologies, 2D gels and LCMS, to fractionate the samples and ran them through mass spectrometers. The 2D gels gave fewer identifications than the LCMS, but it was still a lot of data. For instance, just these results contaied of 30,000 peptide identifications! You can reduce that to about 2,000 proteins that the journal has asked for annotated spectra. Needless to say, the brute force method of taking screen shots of each spectra from the program would not work.

I wrote a few scripts and libraries to parse the raw data and the final result table to come up with the above figure. This is bringing mzXML, excel, and MGF files together with Ruby, C libraries, and the R statistical tool to produce the nice picture you see, but it took me two weeks to figure out the specifics. How on earth could a regular bencher do this?

I think the journal is in for a rude awakening once the backlash of angry rebuttals from paper submitters start flowing in. I would also like to see their reaction to the gigantic pile of spectra we are about to send them.

Tuesday, August 28, 2007

Gbarcode using GD script

BTW, here is the script I used to create the barcode in the previous post using the gbarcode and gd2 gems:


require 'rubygems'
require 'gd2'
require 'gbarcode'

include Gbarcode
include GD2

b = barcode_create("TEST1234567890")
barcode_encode(b,BARCODE_128)

w = 20
h = 100
x = 10

y1 = 10
y2 = h - 20

bars = b.partial.split(//).map {|e| e.to_i}
bars.map {|e| w += e}

i = Image::IndexedColor.new(w,h)
i.palette << c =" Canvas.new(i)" color =" Color::BLACK" font =" Font::Small" f =" File.open(">

Monday, August 20, 2007

The SNP per Gene count

IN my last post I related an example where a scientist came to me to parse a file for the number of SNPs per gene in an excel file. The simplest solution would be to use a hash keyed on the gene symbol and the value tracks the number of times you have seen a particular gene symbol. Here is the program:


require 'rubygems'
require 'fastercsv'

genecount= Hash.new()
FasterCSV.foreach(ARGV[0], :headers => true) do |row|
# headers => id,snp_id,genome_build,chromosome,coordinate,gene_symbol,priority,snp_per_gene
if (genecount[row["gene_symbol"]])
genecount[row["gene_symbol"]] += 1
else
genecount[row["gene_symbol"]] = 1
end
end


output = File.open("#{ARGV[0]}.rev.csv", "w")
output.puts("gene_symbol,snp_count")
genecount.each_pair do |g,c|
output.puts "\"#{g}\",#{c}"
end
output.close

Friday, August 17, 2007

Hacks Before Code

I often find that when you are trying to solve two problems at once, you do a poor job of both. Case in point, someone just came into my office asking how they would go about getting the number of SNPs per gene from some excel file they have. I start to explain set theory and databases and you could see visible signs of mental shutdown ensue (the slacking jaw, the glazed eyes). Trying a different tactic and showing them a script as I wrote it to create a hash keyed by gene and the value being the count of SNPs from the file gave equal results.

So I am trying out Something New. I am going to push that researchers learn to program in a context that is completely separate from science, and is hopefully fun enough that they stick with it for more than a few days. Enter Hackety.org, a project spearheaded by _why the lucky stiff that seeks to (insert Fake Steve Jobs "voice") re-instill the child like wonder back into learning how to program.

With HacketyHack, I hope that researchers are motivated to learn aspects of programming in an entertaining environment before they have to do any real work, which of course will suck some of the fun out of the activity.

I'll be putting together lessons to augment the existing 7 exercises of HacketyHack in the coming months with real but simple bioinformatics tutorials. So download that hack-box and get coding folks!