b.atcg.us

Python multiprocessing - multiple producers, single consumer

2011-11-05T10:58:59-07:00

I'm working on some code that I initially wrote to write results to a mysql database. I chose mysql because it supports concurrent writes, and this program processess data in parallel having each process write it's results to the database when a given task is complete. This was a lazy way of getting the data stored relatively quickly and easily (let's forget about the overhead of any given process for a moment). That said, having mysql as a dependency of your code is a bummer, particularly for folks with one-off tasks who don't have the time or patience to install and configure mysqld. Additionally, portability of mysql data, while relatively easy using dumpfiles, is not as useful as several other options.

Because I'm working in Python, I can also use the exceptional sqlite3 module. It offers most of the database functions that i want/need, it's generally available on any platform, it is open source, and databases are portable between machines.

However, sqlite3 does not support concurrent writes - meaning that I need a general way to process my data in parallel while writing my results to a database using a single process. This basically translates to needing a multiprocessing model having multiple producers and a single consumer - something with few examples available on the interwebs.

So, after a bit of playing around, here's some test code that does just that:

(Relatively) Easily get coverage for velvet assemblies

2011-07-02T11:15:59-07:00

You can get kmer coverage from contigs assembled by velvet by parsing the kmer value from the output fasta header, but sometimes I want "actual" coverage for contigs or coverage across a specific subset of contigs.

Here is a way to do this relatively painlessly (requires that you first download and build the amos tools).

beast in the cloud

2011-04-09T12:58:34-07:00

We've been running beast and mrbayes on several data sets lately, generally using ec2 to help us run multiple analyses simultaneously. Along those lines, I was interested in getting beast (using the beagle-lib) running on ec2, to take advantage of their GPU HPC options (what a load of acronyms!).

Anyway, here are steps to get BEAST running on ec2. Also, expect something a little more formal in the (hopefully) near future describing this option, and even easier means of implementation.

454 assembly on ec2

2011-02-11T15:59:51-08:00

I recently came across a situation in which I needed to assemble some reads from a critter that we gathered from a couple of 454 sequencing runs. Unfortunately, the assembly of these runs using gsAssembler 2.5.3 required a bit more RAM than available on our local workstation running a Fedora image in VMware.

So, I decided to assemble the reads using Amazon EC2. The process was relatively painless once I found an OS image that would boot properly with the instance type I chose (m2.2xlarge. see here for a list).

Here's the steps to follow using a 64-bit Amazon Linux image:

Getting taxonomy information from NCBI

2011-01-29T11:01:51-08:00

Sometimes you need to get taxonomy information from NCBI, assuming that you know a particular species name. If you only are working with one species, then this is not very hard. When it comes to working with multiple species, however, attempting such a task using the web-frontend would be painful. Luckily, you can use the NCBI eutils to make the process somewhat easier, particularly when combined with BioPython:

Several things to remember here:

you must add your email address as the Entrez.email attribute
you should not request info. for more than 100 species during the work-day (see eutils restrictions)
you will need to make changes to get this to work with your particular input

sphinx + github with no submodules

2011-01-22T18:59:50-08:00

Previously, I detailed how I was using git submodules and gh-pages at github to host the html versions of my documentation generated by sphinx. Basically, the problem is that using git submodules for this is a real pain (for me, at least). You've always got to remember to sync up the submodule(s) just so, or things go horribly wrong and stop tracking each other.

Because I can never remember the correct order of what needs to be committed when, I searched for other ways. And searched. And searched. Finally, I ginned up something that I think I like. It basically uses a directory, within your project but ignored by .gitignore, to keep track of your gh-pages content. Because the .git directory tracks what you need to track, you can run make html for sphinx with impunity, and you can push your changes to gh-pages when you want - no need to remember what goes when and where.

Finally, you can track the raw rst of your doc files in the regular source tree, just like you should.

casting a numpy array of strings to int

2010-12-13T09:27:59-08:00

Sometimes you need to create an array from a string, and then you need to cast the array (which is of string type) into something more useful like int - for example when reading PHRED quality scores from a file. You can do this several ways, often using a list comprehension, perhaps like so:

import numpy

s = '40 40 40 40 40'
sl = s.rstrip().split(' ')
si = [int(elem) for elem in sl]
sa = numpy.array(si)

But, this is kludgy and there is a more efficient way:

import numpy

s = '40 40 40 40 40'
sa = numpy.array(s.rstrip().split(' ')).astype(int)

chunking a fasta file, part 2

2010-12-11T15:45:45-08:00

Well, it took me more time than I had planned to get around to wrapping this up... but, it is what it is.

I have completed some code that will use single- or multiple-processes to split a fasta or fastq file into a requested number of subunits. I have yet to test the speed of the code relative to something like split, but my guess is that it's rather fast, particularly for large files.

I provide a use case at the bottom of the file.

vmware fusion and the cli

2010-12-11T00:00:00-08:00

Today, I needed to shutdown an instance running on VMware Fusion so i could free up some resources. Normally, this is an easy task because I'm sitting at the machine in question. But, I ran out of my office the other day and forgot to shutdown the instance. So, i needed to quit it remotely. Digging into things a bit, I found vmrun, which is the CLI to VMware Fusion. Running vmrun without any options show how to make it work, and shutting down a VM is as simple as:

/Library/Application\ Support/VMware\ Fusion/vmrun stop \
"/Volumes/Data2/VM/Fedora 64-bit.vmwarevm/Fedora 64-bit.vmx" soft

a plethora of sequence tags

2010-11-20T00:00:00-08:00

Sequence tags can be attached to DNA reads of interest to let you track different pools of reads following a second generation sequencing run. The best way to generate tags for these reads is a matter of some debate:

you can "just make them"
you can generate them with a program
you can use Hamming distance tags
you can use Levenshtein distance tags

The problem with "just making" the tags is that sequencing errors can inadvertently turn one tag into another if there is, say, an erroneous substitution of a base into the tag portion of the sequence read. Error correcting tags using Hamming distance attempt to counter this effect, but are only robust to substitution errors, which can be problematic. Levenshtein distance tags are robust to insertion, deletion, and substitution error, but it is often hard to find available sets of Levenshtein distance sequence tags.

With all of that in mind, I offer several sets of Levenshtein distance sequence tags. These tags range from 4 to 10 nt and edit distance 3 to 9. The 10nt tags are somewhat slow to create (70 or 80 hours on a multicore machine), so you might as well just use these rather than generate a set, de novo. If you would like to check the tags, to ensure they are of the appropriate distance, you can.

For those interested in the nitty-gritty details, see the code, which is one program within edittag. Now, here are the tags:

extending the nextera indexing set

2010-11-09T00:00:00-08:00

UPDATE: For those of you unwilling to just jump at these, we have validated all of these extended indices and have not found problems with any of the full set of 24, across 6 lanes of multiplexed Illumina sequencing. You may also be interested in edittag, which provides a more generic solution to the barcode generation problem.

The Epicentre Nextera kit comes with aliquots of barcodes that are compatible with the Illumina multiplex system. However, they only send you adapters with 12 different indices:

#Index  Name
ATCACG  IDX_1
CGATGT  IDX_2
TTAGGC  IDX_3
TGACCA  IDX_4
ACAGTG  IDX_5
GCCAAT  IDX_6
CAGATC  IDX_7
ACTTGA  IDX_8
GATCAG  IDX_9
TAGCTT  IDX_10
GGCTAC  IDX_11
CTTGTA  IDX_12

For whatever reason (HiSeq, anyone?), we may want to run a number of indexed sequences in a single sequencing run larger than 12. Thus, we want to extend the set of 12 adapters to support a larger number of samples. Skip to the results, here are the indices/barcodes that extend the set of 12:

Extended Nextera Tags (6nt, edit distance = 2)

Now, as you'll see below, the default set of barcode IDX adapters are edit distance 2 apart from one another. I don't typically feel comfortable until barcode adapters are edit distance ≥ 3 from one another - so, if we want to extend the Nextera set (using most of the adapters provided), we can DROP IDX 1 and IDX 8 from the set of adapters that came with the kit, and then design additional adapters from the barcode sequences int the following (see below for slightly more detail):

Extended Nextera Tags (6nt, edit distance = 3)

I didn't have tons of time to spend on this, so what you can see over at github is slightly ugly. But, it is what it is. Also, I put the Methods posted over on github to catch any changes, but the <script> formatting that you get (notice below) is not quite WYSIWYG.

an alternative method to run colony2

2010-10-20T00:00:00-07:00

For parentage/sibship inference, I've started using Colony2 and MasterBayes in place of the venerable Cervus. However, one (of two*) things that is annoying about Colony2 are the scripts for running the program that are available on "alternative" operating systems (e.g. not Windows). These scripts are provided as part of the rcolony package.

Because the script to run colony through R (aka run.colony()) was causing me problems and because I'd rather do pretty much anything than hack on R code, I wrote some Python code to accomplish a function similar to run.colony() but that operates outside of R. Here goes:

* The second thing that can be a downer when using Colony on alternative operating systems is the formatting of the input file. I've got something in the works, there, too... it's just not at all on the front burner

compiling pslCDNAFilter on osx

2010-10-18T00:00:00-07:00

I need to filter some psl output from blat, and I'd like to do it using the pslCDNAFilter that's provided as part of the Kent Source tree. Although the Kent Source typically behaves well after you set your environment variables appropriately, this little bit was not behaving well. So, I undertook a horrific hack job to get it to compile (this basically requires altering the makefile and physically placing some text into pslCDNAFilter.c):

installing mira on osx

2010-10-16T00:00:00-07:00

I need to assemble some 454 RNA-seq data. I've evaluated a couple of assembly methods, and I'd like to give mira a shot. The output of ./configure indicates that this build can be problematic on osx, so here's what i did:

mysql date and time strings to datetime values

2010-10-14T00:00:00-07:00

Getting dates and times into MySQL can be a total pain. Sometimes, it's easier to just import the date values into a varchar column and do the same with time values - and then convert them to a datetime type (which is far more handy) than a date and a time column (since you can derive both from a datetime column).

So, in my case, I have this setup:

mysql> select usdi, capture_date, time from all_bands_temp 
    -> order by RAND() limit 10;
+-----------+--------------+------+
| usdi      | capture_date | time |
+-----------+--------------+------+
| 154187574 | 5/7/98       | 1312 | 
| 139107796 | 5/14/90      | 0820 | 
| 136141001 | 5/2/85       | 1209 | 
| 187140505 | 5/31/05      | 0750 | 
| 190134590 | 5/9/06       | 0000 | 
| 190134700 | 7/10/06      | 0000 | 
| 139105075 | 4/28/88      | 1135 | 
| 139107358 | 8/11/89      | 1017 | 
| 190134254 | 4/29/04      | 0000 | 
| 154189988 | 5/20/03      | 1322 | 
+-----------+--------------+------+
10 rows in set (0.12 sec)

To convert these into a datetime value, all i need to do is:

mysql> select usdi, capture_date, time,
    -> str_to_date(concat(capture_date, ' ', time), '%m/%d/%y %k%i') 
    -> as str_to_datetime
    -> from all_bands_temp order by RAND() limit 10;
+-----------+--------------+------+---------------------+
| usdi      | capture_date | time | str_to_datetime     |
+-----------+--------------+------+---------------------+
| 139106995 | 6/12/91      | 0707 | 1991-06-12 07:07:00 | 
| 149143944 | 6/2/95       | 1212 | 1995-06-02 12:12:00 | 
| 139107786 | 5/12/90      | 1730 | 1990-05-12 17:30:00 | 
| 147164369 | 5/24/91      | 1030 | 1991-05-24 10:30:00 | 
| 149143925 | 5/28/95      | 1140 | 1995-05-28 11:40:00 | 
| 154189201 | 5/5/02       | 1338 | 2002-05-05 13:38:00 | 
| 123156987 | NULL         | NULL | NULL                | 
| 190134380 | 6/11/04      | 0000 | 2004-06-11 00:00:00 | 
| 139109340 | 7/10/87      | 0830 | 1987-07-10 08:30:00 | 
| 154189917 | 5/3/03       | 1332 | 2003-05-03 13:32:00 | 
+-----------+--------------+------+---------------------+

Now that we know that, all we need to do is insert these values into a datetime column...

two sixty two eighty

2010-10-11T00:00:00-07:00

We have a nanodrop in our lab. It serves us extremely well, and i really appreciate it - i was initially "broken in" on an older spectrophotometer (Beckman DU-800) that basically required you to sacrifice half your sample to get a so-so measure of DNA concentration.

Because we have a nanodrop in our lab, people often ask me about 260/280 ratios. Well, more commonly, they ask me about 260/280 ratios when things go poorly.

First, i should state that the gold standard for whether or not a DNA extraction is "good" (i,e. assuming you are doing PCR) is whether or not you can get it to amplify... with or without trickeration.

Second, for those still interested in 260/280 ratios and the presence of contaminants, i offer this handy little publication (from nanodrop/thermofisher), which i'll just refer to as number nine.

starting a jekyll post from the CLI

2010-10-03T00:00:00-07:00

Starting a jekyll post is one of the minor annoyances of using jekyll (otherwise there are few). However, you can certainly script the redundant task. And, there are several scripts out there, but few worked well for me without throwing an error or two. So, I forked one script, added pieces of another, and the outcome is below:

chunking a fasta file, part 1

2010-10-03T00:00:00-07:00

I've been thinking for a little while about how best to chunk up a gigantic fasta file for distribution across several machines (more on that later). There are obviously several ways to do this - one of which would just be to sequentially read x number of fasta entries (say, 10,000 with some fasta parser) and split them off into a file. Then, rinse and repeat.

However, when we're talking about larger files (let's assume something in the GB range), we may actually want to generate a tuple of start and stop positions within a file (such that each block contains ~ 1 MB of data), and then use this tuple of positions to spread file splitting duties among more than one thread/processor - by having each process grab a set of coordinates, extract a particular part of a file, and write that piece to an outfile. Let's assume we're splitting something like the following:

1. this is my first line\n
2. this is my second line\n
3. this is my third line\n

For this sort of data, we can take the approach discussed in the widefinder benchmarks. Namely, we can write some code to seek through the file in 1 MB chunks, make sure the line we're on is actually a line, and record the start and end positions that split the file into 1 MB chunks/blocks.

However, fasta (and similar files like fastq) present an annoying challenge, by their very design... fasta files look like this:

>fasta-sequence-number-one
ATTTACCCTTTATTTGTCCAGTTGACGTTTTACTTTTATAGATAAATGTTTTTAAGAAGT
TATTGGGCGTTGCTACGAGTGATTGGTAAATACCTTATTGTTTTACTATGTCATGAAGTG
TGACGTACGTGTCATCCTATTTAAAACTTGTCAGTTGAATGTATCTGCATTCTTGGAGTT
>fasta-sequence-number-two
TGACGTACGTGTCATCCTATTTAAAACTTGTCAGTTGAATGTATCTGCATTCTTGGAGTT
ATTTACCCTTTATTTGTCCAGTTGACGTTTTACTTTTATAGATAAATGTTTTTAAGAAGT
TATTGGGCGTTGCTACGAGTGATTGGTAAATACCTTATTGTTTTACTATGTCATGAAGTG

What should be apparent, now, is the crux of the problem... the header line of a fasta file is followed by one, two, or many additional lines of sequence data - meaning that the header line (which begins with '>') and the trailing lines (which are not prepended with a symbol) are a unit that goes together. IN other words, this fasta unit could be 2 lines long (1 header line; 1 sequence line) or 20,000 lines long (1 header line and 19,999 sequence lines). Compounding the formatting nightmare is the additional fact that various fasta-formatting schemes will split the sequence portion of the read across an arbitrary length of characters, to make everything pretty.

Thus, the technique used in the widefinder benchmarks - assuming lines of a file are independent while splitting files into chunks of arbitrary size can cause grave problems when applied to DNA sequence data in fasta format: using this method, we could inadvertently separate (1) the header from it's sequence, (2) the header and line one from lines 2...n, (3) the header and lines 2-3 from lines 4...n, etc., etc.

So, what we need is a method to chunk fasta files that is fast (e.g. ideally faster than sequential read-write operations); efficient (e.g. chunking the file into smaller pieces); and accurate (e.g. not splitting a fasta mid-record).

With that in mind, I offer up fasta_chunker.py, the first of a series of steps that will read a fasta file, returning a iterator of start and offset positions of a chunk that is n (where default = 1) megabytes in size:

bulk sequence renaming with biopython

2010-08-31T00:00:00-07:00

Sometimes large, genome-size files are not named like you want them to be. Most of them are large and it's easier to rename them programmatically. Here is one way to do that, using BioPython.

from Bio import SeqIO

# create a dict to hold our new GI:names mapping, which looks like so, in 
# this case (from sheep)
#
# gi|289623201|gb|CM000885.1|   299839927   chr1
# gi|289623190|gb|CM000894.1|   94216033     chr10
# gi|289623189|gb|CM000895.1|   67137890     chr11
# gi|289623188|gb|CM000896.1|   86457535     chr12
data = {}

# read in this file, split it into a dict
for line in open('oviAri1 copy.info','rU').readlines():
    line_split = line.strip('\n').split('\t')
    data[line_split[0]] = line_split[2]

# open a file for the output
output_file = open('new_oviAri1_seq.fa','w')

# create an iterable to hold the new data
new_seq = []

# iterate over seq, updating the name.  This is going to give us something like:
#
# >chr1
# 
for record in SeqIO.parse(open('oviAri1.fa','rU'),'fasta'):
    new_record_name = data[record.id]
    record.id = new_record_name
    record.name = ''
    record.description = ''
    new_seq.append(record)

# write the whole thing out
SeqIO.write(new_seq, output_file, 'fasta')

restriction batches and biopython

2010-07-20T00:00:00-07:00

Sometimes, you need to get an idea of the sizes of fragments you are likely to get from a sequence cut with a batch of restriction enzymes. Here's how you would do that virtually.

This might be useful when you have no idea of the likely size of resulting fragments, but you'd like to have some idea - before you order a lot of enzymes.

# subset record #

test=seq[0:10000]

# take a look at the common restriction enzymes
Restriction.CommOnly

# do a restriction batch analysis on our test sequence with ALL common
# enzymes
Ana = Restriction.Analysis(Restriction.CommOnly, test.seq, linear=True)

# look that the enzymes that cut (BLUNT) the test sequence in more than #
# 10 spots
results = []
for r in Ana.blunt():
    if len(Ana.blunt()[r]) > 10:
        results.append([r, len(Ana.blunt()[r])])

# sort by the 2nd element of the list
for r in sorted(results, key=lambda record: record[1]):
    print r[0], r[1]
    
# or you can do this:
from operator import itemgetter
for r in sorted(results, key=itemgetter(1)):
    print r[0], r[1]

virtual restriction digests with biopython

2010-07-16T00:00:00-07:00

The other day, I needed to get an approximate idea of the size of fragments returned from a particular restriction digest. Generally, this sort of info is available on the intertubes, but I was not finding terribly much.

So, I decided to do a quick, virtual restriction digest using Python, BioPython, and numpy for some summary stats (in this case I "digested" chr1 of Arabidopsis thaliana):

from Bio import SeqIO
seq = SeqIO.read(open('chr1.fas','rU'),'fasta')

len(seq)
30427671

from Bio import Restriction
Restriction.HindIII.search(seq.seq)
sites = Restriction.HindIII.search(seq.seq)

# get the number of fragments produced by cutting
len(sites)
16963

# sort the sites by position (just to be sure)
sites.sort()

# get the distance between sites
dist = [sites[r] - sites[r-1] for r in xrange(1,len(sites))]

# generate some summary stats
import numpy
import math
dist_a = numpy.array(dist)

# mean fragment size
numpy.mean(dist_a)
1793.6787525056006

# 95 % Confidence invterval around mean
1.96 * numpy.std(dist_a, ddof=1)/math.sqrt(len(dist_a))
33.572399955020472

# get number of fragments btw. 100 and 500 bp - this could likely use some #
# cleanup, but I was in a hurry
less_than_five_bool    = dist_a < 500
less_than_five         = dist_a[less_than_five_bool]
greater_than_one_bool  = less_than_five > 100
one_to_five_hundred_bp = less_than_five[greater_than_one_bool]

removing a file/directory permanently from git

2010-06-29T00:00:00-07:00

I'm not entirely sure how many times I've needed this or how many times I've dug it up only to forget it later, but here is how to permanently remove a file or directory from a git repository. This comes by way of Dalibor Nasevic's site.

# blackhole the content
git filter-branch --tree-filter 'rm -rf my/folder' HEAD

# force the update
git push origin master --force

If you like to read up on these things, you can find more information in Chapter 6 of the ProGit book.

py2app hell, the second

2010-06-14T00:00:00-07:00

I have this funny relationship with py2app. I like what it produces but it is a real hassle to get working. Anyway, the basic situation is this:

i needed an application bundle for msatcommander
you cannot, technically, include the Apple-compiled python in an application bundle
this means you need to get one from somewhere else (e.g. roll-your-own)
pyqt4 requires a framework build of python (unless you hack, hack, hack)

Other interesting tidbits:

you need to keep the architectures of everything produced the same. i chose to stick with x86_64 and i386.
the prepackaged Qt is built for 3 architectures. This will cause you problems.

Therefore:

matplotlib hell, too

2010-06-14T00:00:00-07:00

For another thing I'm working on, I needed matplotlib to compile (this is on a slightly difference machine than that of the previous post). Bottom line is this: I needed matplotlib built for x86_64 and i386 and without too much other mumbo-jumbo. I wanted to build this in a virtualenv, because that is how i like things.

Needless to say, matplotlib building is pretty awful (see previous post). It sort of reminds me of py2app.

Of course, all of that being said, part of the problem is due to Apple switching architectures while maintaining backwards compatible support (at least up to OSX 10.5). Anyway, I was finally able to get a working matplot lib (with help from the stackoverflow post here, using the following OSX-specific makefile (note: i installed freetype and libpng from source into /usr/local/)

And the following edits to the ./Makefile that comes with matplotlib (not making these edits builds some shared objects for ppc and i386, which causes problems):

Finally, I built matplotlib with:

make -f osx.make mpl_build
make -f osx.make mpl_install

revisiting regex

2010-04-13T00:00:00-07:00

To follow up on a previous post, Google has today released re2, noting the release in a blog post. From the re2 website:

RE2 is a fast, safe, thread-friendly alternative to backtracking regular expression engines like those used in PCRE, Perl, and Python. It is a C++ library.

py2app hell, the first

2010-04-13T00:00:00-07:00

I ran into a problem with a program that I work on where wxPython wasn't cutting it any longer, for various reasons (notably their current lack of support for 64-bit archs - see here). I wasn't really in love with wxPython, anyway, so I decided to ditch it in favor of Qt4, in combination with PyQt4 (which depends on sip for its C++ bindings).

While the general process of rewriting my GUI interface was pretty straightforward (given the decent tools in the Qt SDK + pyuic from PyQt4), I figured that packaging up the release was going to be painful. Call me crazy, but the shift to 64-bit in Snow Leopard should cause some issues with py2app working well for building bundles. It generally turns out that I was correct (don't get me wrong, I like py2app, it's just touchy).

I searched around for specific instructions on how to fix lots of the things that were broken. These included what I think was an endianness error, a typo in py2app's subversion code, mach-o build mismatches between sip-PyQt4-PyQt, and failure of py2app to add a directory to the bundle you build (i think that's it). Of most utility in my problem solving, I came across this thread which details many of the additions made (in the subversion trunk) to get py2app working, so thank you Marc-Antoine Parent. You are awesome. This post by Piotr Maliński was also helpful.

Long story short, for the moment, on Snow Leopard (10.6.2) with the system python, you need to do the following to make things work:

sphinx documentation and github

2010-02-10T00:00:00-08:00

So, there's tools out there to get your documentation up on github (github-tools, sphinx-to-github), but neither totally jive with how I want to do things. So, I went searching long and hard for a way to:

keep the document-generating source in master
place the generated html in gh-pages, so it's served as static html
have everything just "work"

Luckily, I came across this post which gets you most of the way there (but assumes you're using jekyll, which isn't the case for me/us).

So…

climm and osx

2010-01-31T00:00:00-08:00

I've been looking for a CLI chat program that works pretty well. After digging around for some time, I came across climm. climm has a lot of dependencies, including:

tcl
iksemel
libgcrypt
openssl
pkgconfig
gloox
gnutls
libotr
autoconf
automake
libtool

You can get most of these from macports (if you install libgcrypt, you should pull down pretty much everything). You will still need to port install libotr and you will also need to download the source for gloox (1.0) as the version at macports fails to build.

Once all of that is done, get the source for climm, patch the source of configure.ac with the patch referenced below (climm-r2826.patch) and rm configure.

You also need to patch climm-0.7/src/io/io_dns.c replacing

#include <arpa/nameser.h>

with

#include <arpa/nameser_compat.h>

See climm-nameser_compat.patch. Now, you want to

make clean
./configure
make
sudo make install

Follow the instructions from the wizard. climm initially segfaulted after finishing the wizard, but appears to be "fine" on a restart.

Update: I should have mentioned that these changes were made to the 0.7.0 source distribution.

moving folder from an existing git repo

2010-01-16T00:00:00-08:00

Today, I needed to move some code from another project into it's own repo. I really wanted to keep the history of the changes to the code, during its development, so I went a googlin'. Here's what I found:

http://bit.ly/7fUfBy

Worked like a charm.

another mini next-gen sequencer

2010-01-15T00:00:00-08:00

Not to be outdone by Roche, Illumina has announced a new version of it's Genome Analyzer: the IIe. The IIe's specs are generally pretty reasonable (18-20 GB in about 10 days @ 2x100). Cost will be approximately $250,000.

Ilumina also annouced a big brother to the IIx, the HiSeq2000, which kicks things up a notch to 150-200 Gb (of 2x100) in about 8 days.

Wow.

pianobar - a CLI for pandora

2009-12-31T00:00:00-08:00

Just ran across this - Pianobar - which is a CLI for Pandora. This builds pretty easily on the mac. Following (with a modification) instructions from this thread over at reddit.com:

~~Install MacPorts or make sure it's squared away with this gist~~

Update: it's fantastic!

Update 2: just install it using homebrew:

brew install pianobar

temporarily fix wxPython on osx 10.6

2009-12-23T00:00:00-08:00

Just run this in the terminal to change to 32-bit Python for the session or set your environment variable in .bash_profile. Once you want to change back, just erase it. No fiddling around with defaults write com.blah.blah.blah:

python multiprocessing.JoinableQueue() example

2009-12-19T00:00:00-08:00

After a lot of pain and suffering trying to debug what was going wrong that appeared to be a stuck process, it simply boiled down to not closing my Queues:

Here's the full module's documentation: http://docs.python.org/library/multiprocessing.html

of parsers, regex and python...

2009-12-11T00:00:00-08:00

After re-reading the excellent Regular expression matching can be simple and fast by Russ Cox and noticing a reddit.com comment related to parsing expression grammars, I went on a google and came across the following handy summary of python parsers:

Python parsing tools

add google calendar items from launchbar

2009-12-11T00:00:00-08:00

I came across a post on Launchbar's forum, discussing a means of adding entries to your Google Calendar via Launchbar. I've made a small simplification to the script that runs the terminal command in the background (with less error checking, of course). Follow the instructions up to the script presented in the forum post (install GData and python-dateutil, download and place the gcalcli script somewhere) and copy the following into AppleScript Editor:

Then edit the first line to point to the gcalcli binary and save the script to ~/Library/Application Support/LaunchBar/Actions/.

I may eventually add some other modifications, although Applescript is not my bag.

Update: Google Calendar Quick Add Grammar

sql histograms

2009-12-09T00:00:00-08:00

via Michael Stassen @ http://lists.mysql.com/mysql/155540

rotate labels in ggplot2

2009-12-09T00:00:00-08:00

Rotating labels in ggplot2 can be mystery. Maybe this helps:

median and mode in mysql

2009-12-09T00:00:00-08:00

it's just so cute... Roche 454's gsjunior

2009-12-09T00:00:00-08:00

I may be behind the times a bit, but I’ve just run across the cutest little circa $100,000 (? my guess ?) DNA sequencer in the land (Press Release.)

You can see some opinion on the machine over at Pathogenomics. All points presented are certainly true:

library construction needs to be easier
reagents need to be less expensive
etc.

One thing that’s not mentioned, where I think this machine will be quite useful, is in the methods development and/or amplicon/sequence capture spaces - essentially those areas in which we are either trying to or succeeding at changing the initial focus of 454 sequencing from vertical (depth of coverage) to horizontal (breadth of coverage).

a smaller platform for testing new sequencing techniques would be super handy - particularly if/when library prep is simplified and costs are reduced. Overall, small-run efficiency should be greater as you’re using the entirety of the plate for each run (i,e. you won’t loose reads from areas of the plate covered by a gasket) and costs should be lower as a result of increased efficiencies and direct reductions in time and reagent prices. I guess that I am sort of thinking of it as a rapid application development tool, but within the genomics/genetics sphere. The assumption would be that once the method is working, you scale up to the “full” 454 run.
Within the amplicon/sequence capture sp here, the machine should be of sufficient size for daily to weekly runs and likely of sufficient size to track the workflow of several lab employees. Assuming that the gsjunior will return roughly 80,000 reads (from the 100k spec.), we get a theoretical coverage of 10X across 1000 loci from 8 individuals. Similarly, if we expand to 96 individuals, we’re going to get around 10X coverage of 80 loci (or @ 8X we can get roughly 100 loci per 96 individuals).