Bio and Geo Informatics

Adding bed/wig data to dalliance genome browser

2011-04-12T18:40:00.000-07:00

I have been playing a bit with the dalliance genome browser. It is quite useful and I have started using it to generate links to send to researchers to show regions of interest we find from bioinformatics analyses.
I added a document to my github repo describing how to display a bed file in the browser. That rst is here and displayed in inline below.
It uses the UCSC binaries for creating BigWig/BigBed files because dalliance can request a subset of the data without downloading the entire file given the correct apache configuration (also described below).
This will require a recent version of dalliance because there was a bug in the BigBed parsing until recently.

Dalliance Data Tutorial

dalliance is a web-based scrolling genome-browser. It can display data from
remote DAS servers or local or remote BigWig or BigBed files.
This will cover how to set up an html page that links to remote DAS services.
It will also show how to create and serve BigWig and BigBed files.

Note

This document will be using hg18 for this tutorial, but it is applicable to
any version available from your favorite database or DAS .

Creating A BigBed

Getting a bed file from UCSC

From the UCSC table browser choose
genome: Human

assembly: NCBI36/hg18

group: Genes and Gene Prediction Tracks

track: UCSC Genes

table: knownGene

output format "selected fileds from primary and related tables"

in text box, name it "knownGene.hg18.stuff.txt"

click "get output"

check kgXref under 'Linked Tables'

click 'Allow Selection From Checked Tables' at bottom of page.

check 'geneSymbol' from hg18.kgXref fields section

click 'get output' and a file named 'knownGene.hg18.stuff.txt' will be saved to your downloads directory. move it to your current directory.
To get this into bed format copy and paste this onto the command-line:
grep -v '#' knownGene.hg18.stuff.txt | awk 'BEGIN { OFS = "\t"; } ;
{   split($9, astarts, /,/);
    split($10, aends, /,/);
    starts=""
    ends=""
    for(i in astarts){
        if (! astarts[i]) continue
        ends=ends(aends[i] - astarts[i])","
        starts=starts(astarts[i] = astarts[i] - $4)","
    }
    print $2,$4,$5,$1","toupper($13),1,$3,$6,$5,".",$8,ends,starts
}' | sort -k1,1 -k2,2n > knownGene.hg18.bed
To create a BigBed from this, do (note if you're not on a 64 bit
machine, you'll have to find the 32bit binaries:
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/fetchChromSizes
wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bedToBigBed
chmod +x fetchChromSizes bedToBigBed
./fetchChromSizes hg18 > data/hg18.chrom.sizes
./bedToBigBed knownGene.hg18.bed data/hg18.chrom.sizes knownGene.hg18.bb

now knownGene.hg18.bb is a BigBed file containing both the UCSC and the common
name in the name column.

SQL

UCSC also has a public mysql server so the process of downloading to a bed can be simplified to:

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg18 -P 3306   -e "select chrom,txStart,txEnd,K.name,X.geneSymbol,strand,exonStarts,exonEnds from knownGene as K,kgXref as X where  X.kgId=K.name;" > tmp.notbed
grep -v txStart tmp.notbed | awk 'BEGIN { OFS = "\t"; } ;
    {   split($7, astarts, /,/);
        split($8, aends, /,/);
        starts=""
        sizes=""
        exonCount=0
        for(i in astarts){
            if (! astarts[i]) continue
            sizes=sizes""(aends[i] - astarts[i])","
            starts=starts""(astarts[i] = astarts[i] - $2)","
            exonCount=exonCount + 1
        }
        print $1,$2,$3,$4","$5,1,$6,$2,$3,".",exonCount,sizes,starts
    }' | sort -k1,1 -k2,2n > knownGene.hg18.bed

then proceed as the last steps above to create the big bed file.

Displaying A BigBed in Dalliance

From there, download dalliance:

$ git://github.com/dasmoth/dalliance.git
cd dalliance

and edit test.html, adding:

{name: 'UCSC Genes',
 bwgURI:               '/dalliance/knownGene.hg18.bb',
},

before the line that looks like:

{name: 'Repeats',

at around line 55.

Then edit your apache.conf (e.g. /etc/apache2/sites-enabled/000-default)
and put the following
(here i assume you cloned dalliance into /usr/usr/local/src/dalliance-git):

Alias /dalliance "/usr/local/src/dalliance-git"
<Directory "/usr/locals/src/dalliance-git">
 
    Header set Access-Control-Allow-Origin "*"
    Header set Access-Control-Allow-Headers "Range"
 
    Options Indexes MultiViews FollowSymLinks
    AllowOverride None
    Order allow,deny
    Allow from all
</Directory>

Then enable mod-headers apache module. On Ubuntu, that looks like:

sudo a2enmod headers

Then point your browser to:: http://yourhost/dalliance/test.html
And you should see the your 'UCSC Genes' track in full glory along
with the other niceties of the dalliance browser.

(bloom) filter-ing repeated reads

2010-10-22T09:12:00.000-07:00

In this post, I'll talk a bit about using a bloom filter as a pre-filter for large amounts of data, specifically some next-gen sequencing reads.

Bloom Filters

A Bloom Filter is a memory efficient way of determining if an element is in a set. It can have false positives, but not false negatives. A while ago, I wrote a Cython/Python wrapper for the C code that powers the perl module, Bloom::Filter. It's has a nice API and seems very fast. It allows specifying the false positive rate. As with any bloom-filter there's a tradeoff between the amount of memory used and the expected number of false positives.
The code for that wrapper is in my github, here.

Big Data

It's a common request to filter out repeated reads from next-gen sequencing data. Just see this question on biostar. My answer, and every answer in that thread, uses a tool that must read all the sequences into memory. This is an easy problem to solve in any language, just read the records into a dict/hash structure and use that to find duplicates and print out only the best or first or whatever. However, once you start getting "next-gen", this is less useful. For example, someone kindly reported a bug on my simple c++ filter because he had 84Gigs of reads and was running out of memory. And that is a bug. Anything that's supposed to deal with next-gen sequences has to deal with that stuff.

As a side note, my friend/co-worker came up with an elegant solution for filtering larger files: split into 4 files based on the first base in the read (A, C, T, or G) then filter each of those files to uniqueness and merge. I like that approach, and as the file sizes grow it could be extended to separate into 16 files by reading the first 2 bases. But, it does require a lot of extra file io.

Just Do It

So, I wrote a simple script that filters to unique FASTQ reads using a bloom-filter in front of a python set. Basically only stuff that is flagged as appearing in the bloom-filter is added to the set. This trades speed--it iterates over the file 3 times--for memory. The amount of memory is tuneable by the specified error-rate. It's not pretty, but it should be simple enough to demonstrate what's going on. It only reads from stdin and writes to stdout, with some information about total reads an number of false positives in the bloom-filter sent to stderr.
usage looks like:
python fastq-unique.py > in.fastq < out.unique.fastq
On my machine, a particular run with a decent sized file looks like this:

and here's the code:

filtering paired end reads (high throughput sequencing)

2010-09-20T16:34:00.000-07:00

NOTE: I don't recommend using this code. It is not supported and currently does not work for some sets of reads. If you use it, be prepared to fix it.

I wrote last time about a pipeline for high-throughput sequence data. In it, I mentioned that the fastx toolkit works well for filtering but does not handle paired end reads. The problem is that you can filter each end (file) of reads independently, but most aligners expect that the nth record in file 1 will be the pair of the nth record in file 2. That may not be the case if one end of the pair is completely removed while the other remains.
At the end of this post is the code for a simple python script that clips an adaptor sequences and trims low-quality bases from paired end reads. It simply calls fastx toolkit (which is assumed to be on your path). It uses fastx_clipper if an adaptor sequence is specified and then pipes the output to fastq_quality_trimmer for each file then loops through the filtered output and keeps only reads that appear in both. Usage is something like:

ADAPTORS=GAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGT,GAAGAGCGGTTCAGCAGGAATGCCGAGACCGATATCGTATGCCGT,GAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCG
pair_fastx_clip_trim.py --sanger -a $ADAPTORS -M 20 -t 28 -l 40 en.wt.1.fastq en.wt.2.fastq

Where the -a (adaptor) -M (length of adaptor match) -t (min quality threshold) and -l (min length after quality chops) options are copied directly from (and sent directly to) fastx toolkit. --sanger indicates that the reads have fastq qualities in the sanger encoding. If that option is not specified, qualities are assumed to be in illumina 1.3 format where the ascii offset is 64.
Output
This example will create 2 new files: en.wt.1.fastq.trim and en.wt.2.fastq.trim each with the same number of corresponding records that pass the filtering above.
Adaptors
As described in my previous post, sometimes there are multiple adaptor sequences in the reads. This script can filter out any number of adaptors--specified in a comma delimited option -a--in a single run.
Script
It's not too pretty, but it does the job:

As always, let me know of any feedback.

ngs / high-throughput sequencing pipeline

2010-09-12T18:21:00.000-07:00

This is the minimal set of preprocessing steps I run on high-throughput sequencing data (mostly from the Illumina sequencers) and then how I prep and view the alignments. If there's something I should add or consider, please let me know.
I'll put it in the form of a shell script that assumes you've got this software installed.
I'll also assume your data is in the FASTQ format. If it's in illumina's qseq format, you can convert to FastQ with this script by sending a list of qseq files as the command-line arguments.
If your data is in color-space, you can just tell bowtie that's the case, but the FASTX stuff below will not apply.
This post will assume we're aligning genomic reads to a reference genome. I may cover bisulfite treated reads and RNA-Seq later, but the initial filtering and visualization will be the same. I also assume you're on *Nix or Mac.

Setup

The programs and files will be declared as follows.

#programs
fastqc=/usr/local/src/fastqc/FastQC/fastqc
bowtie_dir=/usr/local/src/bowtie/bowtie-0.12.7/
samtools=/usr/local/src/samtools/samtools

#files. (you'll want to change these)
FASTQ=/path/to/a.fastq
REFERENCE=/path/to/reference.fasta

Most of the following will run as-is for any set of reads, you'll only need to change the FASTQ and REFERENCE variables above.

Seeing the Reads

If you have a file with 2GB of reads, you likely can't just read (pun intended) it and get an idea of what's going on -- though I've tried. There are a number of tools to give you a summary of the data including stats such as quality per base, nucleotide frequency, etc. While fastx toolkit will do this for you, I've found fastqc to be the best choice.
The command:

$fastqc $FASTQ

will write a folder "a_fastqc/" containing the html report in fastqc_report.html
Here's an example of the nicely formatted and informative FastQC report before quality filtering and trimming (scroll to see the interesting stuff):

from that, we can see that (as with all illumina datasets) quality scores are lower at the 3' end of the read. For this analysis, I'll choose to chop bases with a quality score under 28. In addition to the other information, we can see that there is an Illumina Primer still present on many of the reads, so we'll want to chop that out in the filtering step below.

Filtering the Reads

The fastx toolkit is not extremely fast, but it has a simple command-line interface for most common operations I use to filter reads (though it--like similar libraries--is lacking in support for filtering paired-end reads). For this case, we want to trim nucleotides with quality less than 28 from the ends of each read and then remove reads of length 40. In addition, we want to chip Illumina adaptor sequence. This can be done by piping 2 commands together:

# the sequence identified by FastQC above:
CLIP=GATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAAAAAAA
fastx_clipper -a $CLIP -i $FASTQ -M 30 | fastq_quality_trimmer -t 28 -l 40 > ${FASTQ}.trim

both fastx_clipper and fastq_quality_trimmer are provided in the fastx toolkit.
The first command clips the adaptor sequence from $FASTQ and pipes the output to the 2nd command, fastq_quality_trimmer, which chomps bases with quality less than 28 (-t) then discards and read with a remaining length less than 40 (-l) and sends the output to a .trim file.

Seeing the filtered Reads

We then want to re-run fastqc on the trimmed, clipped reads
The command:

$fastqc ${FASTQ}.trim

and the output looks like:

where we can see the quality looks much better and there are no primer sequences remaining.
This took the number of reads from 14,368,483 to 11,579,512 (and many shorter--trimmed--reads in the latter as well).

Analysis

Up to this point, the analysis will have been very similar for RNA-Seq, BS-Seq, and genomic reads, but you'll want to customize your filtering steps. The following will go through a simple alignment using bowtie that assumes you have genomic reads.

Aligning

The first step is to build the reference bowtie index:

${bowtie_dir}/bowtie-build --quiet $REFERENCE bowtie-index

That will create the bowtie index for your reference fasta in the bowtie-index/ directory. It can take a while so you may want to download the pre-built indexes from the bowtie website if your organism is available.
Then, you can run the actual alignment as:

${bowtie_dir}/bowtie --tryhard --phred64-quals -p 4 -m 1 -S bowtie-index -q ${FASTQ}.trim ${FASTQ}.trim.sam

which tells bowtie to try hard (--tryhard), assume quality scores are the latest schema from illumina (--phred64-quals), use 4 processors (-p 4), discard any reads that map to more than one location the the reference (-m 1), use SAM output format (-S) and then where to find the bowtie index and the reads. The output is sent to ${FASTQ}.trim.sam.
We've done a lot of experimenting with different values for -m, and can affect your results, but -m 1 seems a common choice in the literature. And it's clearly less sketchy than, for example, -m 100 which would report up to 100 alignments for any read that maps to less than 100 locations in the genome.

View the Alignment

From there, we want to view the alignment. Most tools can handle SAM formatted files, but will perform better with an indexed bam. To get that, we do:

# view with samtools filter out unmapped reads and converted to sorted, indexed bam.
${samtools} view -bu -S -F 0x0004 ${FASTQ}.trim.sam -o ${FASTQ}.trim.unsorted.bam
${samtools} sort ${FASTQ}.trim.unsorted.bam ${FASTQ}.trim
${samtools} index ${FASTQ}.trim.bam

where now ${FASTQ}.trim.bam is sorted and indexed and contains only mapped reads (the .sam file from bowtie contains unmapped reads).
You can have a quick look at the alignment stats with:

${samtools} flagstat ${FASTQ}.trim.bam

and you can see an awesome ascii curses view of the alignment with:

${samtools} tview ${FASTQ}.trim.bam

to get something that looks like:

But you'll probably want to use a "real" viewer such as IGV, MochiView, or Tablet, all of which I have had some sucess with.

Note that you may want to choose another aligner, because bowtie does poorly with indels. Something like GSNAP will be better suited for reads where you have more complex variants.
I welcome any feedback on these methods.

GSNAP

2010-07-14T12:50:00.000-07:00

Aligners

Since starting the methylcoder project, I've been using the bowtie short read aligner. It's very fast, uses very little memory, aligns Illimina, SOLID, and colorspace reads, and has enough options to keep you busy (including my favorite: --try-hard).
There's a new short-read aligner in my feed-reader each week. I wish, as a service, they'd tell me what they do that bowtie doesn't. There are 2 scenarios I know of that bowtie doesn't handle.

As with many aligners, bowtie creates an index of the reference genome. It uses a suffix array compressed with the Burrows-Wheeler Transform(BWT). It uses 32 bit integers to store the offsets into that index so it's limited to under 4 gigabases of reference sequence. This is more than enough to hold the human genome and, so is not a problem for most users, but, for BS-treated alignments, I need to store a converted (Cytosine to Thymine) forward, and reverse copy of the genome which doubles the size of the reference and puts the human genome past that limit. The "solution" is to split the reference and create multiple indexes. But, with that setup, when requesting reads that map uniquely to a single best location, it's only guaranteed that the mapping is unique to each index independently and post-processing is required. This is also a problem for any large genome that will not fit into a bowtie index.
Bowtie does not handle indels very well. Even the other popular aligners SOAP and BWA can only handle very short indels (up to about 3 bp) while MAQ and SOAP2 do not align reads with indels.

GSNAP

GSNAP: Genomic Short-read Nucleotide Alignment Program (from Thomas Wu at genentech) addresses both of those short-comings and seems reasonably fast. For the reference index, rather than use a suffix-array, GSNAP uses a hash of 12-mers sampled every 3 basepairs (by my reading, sampling complicates the search somewhat but reduces the size of the index). Since it can then use 12-mer seeds to span gaps, it's better at dealing with indels and can also be used to map RNA-seq data to a reference genome, with or without known, annotated splice sites! In addition, it can take known SNPs and add them to the 12-mer hash index so that reads mapping over a SNP will not be wrongly accused of having a mismatch.
Finally, it can map BS-sulfite treated reads to the genome without requiring the C=>T conversion--using the same re-indexed genome.
It's input format is a bit wonky for paired end reads. Normally, paired-end reads are specified in 2 separate files: file_1.fastq, file_2.fastq with the header (and order) indicating the pairing. GSNAP requires a single FASTA file (it will not accept FASTQ) with format:

>header
ACTCTCAGCGGGACGTTAACGCGACCGATTACGGTGACC
CCACGTGCCGACTTAGGCAGACCGACGTTACGCACCACA

where the first sequence is from the file_1.fastq and the 2nd is from file_2.fastq.
As of yesterday, MethylCoder supports GSNAP as an aligner, in addition to Bowtie. (This addition was quite simple since GSNAP does all the work anyway). MethylCoder includes a script to convert paired-end FASTQ files to GSNAP format. The script simply takes as arguments paired end FASTA or FASTQ files and puts them in a single file of the GSNAP paired end format. It's almost too simple to worry about, but I needed it so it's there...

The GSNAP paper is an interesting read. Ever heard of "galloping binary search"? I hadn't, but apparently it's used in python's timsort.

python tools for bioinformatics II: group-ing

2010-06-07T23:45:00.000-07:00

This is the second post describing python tools and idioms I've found useful in bioinformatics.

This post involves grouping items either using a disjoint set or python's
itertools.groupby. I was introduced to both of these
by the post-doc who sits behind me.

Grouper

Grouper is an implementation of a disjoint set, initially from this recipe by Michael Droettboom and included in matplotlib.cbook.
The slightly modified implementation we use is here.
Where the docstring gives a decent idea of useage:

"""
   This class provides a lightweight way to group arbitrary objects
   together into disjoint sets when a full-blown graph data structure
   would be overkill. 

   Objects can be joined using .join(), tested for connectedness
   using .joined(), and all disjoint sets can be retrieved using list(g)
   The objects being joined must be hashable.

   >>> g = Grouper()
   >>> g.join('a', 'b')
   >>> g.join('b', 'c')
   >>> g.join('d', 'e')
   >>> list(g)
   [['a', 'b', 'c'], ['d', 'e']]
   >>> g.joined('a', 'b')
   True
   >>> g.joined('a', 'c')
   True
   >>> 'f' in g
   False
   >>> g.joined('a', 'd')
   False"""

Basically the idea is if there's an explicit connection between `A` and `B` and an explicit connection between `A` and `C` the Grouper class will create a connection between `B` and `C`. In the implementation above, the explicit connections are created with `Grouper.join()` For genomics stuff, we can use this in many places, one simple use-case is finding local duplicates. These are also called tandem duplications and are identified as a group of genes with very similar sequence co-occurring in a local chromosomal region. We want to group these into a single "family" with all members, even if there is not an explicit connection between all genes--due to slight sequence divergence or sequence alignment (BLAST) oddities.
A snippet of code from here (more nice code from Haibao) shows the use of a grouper to find syntenic regions by checking if adjacent members of (an already sorted list) are within a certain distance (window) of each other:

after this, the elements that are grouped in the `g` Grouper object are in the same window (though not strictly syntenic).
This is a very nice abstraction for anything where you have groups of related objects. It reduces the accounting invovled because once you have `join`ed all the elements, querying the Grouper object with any element will return all elements in the group to which it is joined.

groupby

Python's itertools has a number of useful features. `groupby` is a means to partition an iterable similar groups of adjacent items, where `similar` is determined by a function you supply to groupby. The example from the docs is:

but it's also possible to specify an arbitrary grouping function, so grouping capital letters together:

So the `k`ey tells what the grouping function returned and the values or groups are grouped with other adjacent items from the list that also pass the test. So in some cases, we often want to sort the input before using groupby and then we can get groups and avoid an extra nested for-loop that we would otherwise need for filtering.
An example use is grouping BLAST hits by query and subject. Here, BlastLine is a simple object that takes a line from a tab-delimited blast and converts the e-value and bitscore to float and the starts and stops to ints, etc.

and since the blasts are sorted beforehand, this gives a useful set of blasts in blast_iter--containing only blasts with a matching query and subject. To do this without groupby requires quite a bit more code and reduces readability.

I recently saw a post about using `groupby` to parse fasta files. It's a good idea as a fasta file consists of pairs of a single header line, followed by N lines of sequence. A header line always starts with ">" so the test function for groupby is clear, it's simply getting them out in pairwise fashion, here's one way building from the post above:

groupby-grouper

Here's a final example, both of these used together to find tandem duplications given a genomic Bed file (contains chromsome, start, stop) and a list of blasts. Again this is code from Haibao from here.

So, this creates a `simple_blast` that's easier to sort and group, it iterates for the groups based on the subject
and then groups hits if they're on the same chromosome an within a given distance (where distance is measured in genes).
I like this example because I'd previously written code to do the same thing without Grouper or groupby and it was longer, slower and less readable.

python tools for bioinformatics I: convolve for moving average

2010-05-30T09:19:00.000-07:00

This will be the first post of a few describing python tools or idioms that I have found work to well for doing bioinformatics -- or more specifically genomics.

Convolve

numpy makes a lot of things easier--for example, moving average. We often want to know the GC-content for a given sequence, which is the proportion of DNA sequence in a given window that is made of of (G)uanine or (C)ytosine. You can often see moving average code done with a nested for loop--loop over each nucleotide and then over each surrounding nucleotide in that window. Below is the gist of the simpler, faster method using convolve. It does require you do get your sequence data into a numpy array, but that's easily done from any DNA (or protein) sequence string doing "np.array(sequence_string, dtype='c')".
So the workflow is to

get the sequence into a numpy array of dtype 'c'
convert the boolean array to int so True is 1, False is 0
decide on the window size and make a kernel of that shape that sums to 1
profit

with a quick addition for plotting:
the first 100 kilobases of gc content for 50 and 250bp windows look like:

where, as expected, the 250 basepair window is more "smoothed".
From there, it's possible to do some analysis, for example grab the regions with the highest gc content.
That's some fairly dense code to pull out the values centered in the highest GC-content windows and then show the mean of those windows. (Arabidopsis thaliana, which I used for that example, has a fairly low genome-wide GC-content, so the 0.61 average is quite high.)

Another Example

There was a recent question on biostar where someone wanted to find regions with an abundance of a couple amino acids above some cutoff.
The main part of that looks like:

    abool = (seq == 'N') | (seq == 'Q') # convert to boolean
    mw = np.convolve(abool, kern, mode='same')
    if mw.max() < cutoff: continue

Where that operates on 'N' and 'Q' amino acids in a protein sequence instead of 'G' and 'C' as in the example above.

Kernel

Finally, this approach is much more flexible than I've shown, the kernel does not have to be uniform, it can be guassian, or even only taking values to the left of each nucleotide, or weighting the nearest 10 nucleotides at 1 and the rest at 0.2. Even given those changes, once the kernel is chosen, the rest of the "workflow" remains unchanged.

Notes

In cases where the input array is numeric--not sequence--there is no need to do the conversion to a boolean, simply run the convolution on the original array with the chosen kernel.

In the examples above, I've shown only np.convolve with mode='same', for most cases dealing with sequence, I think this is a good choice, but it's best to consult the documentation for each specific case.

Finally, in cases where the kernel and the input array are very large, it may be better to use fft convolution from fftconvolve in scipy. I haven't used this much, I think it may require doing your own padding and chopping at the edges...

numpy + GDAL = agoodle

2010-05-19T17:44:00.000-07:00

It's been a while since I've posted anything geo-related. Just want to let folks know about agoodle a project that makes it easy to access raster geo files as numpy arrays (thanks to GDAL), and to query them with polygons. The simplest way to see it in action is to go to landsummary, a project Josh Livni and I worked on a couple years ago, and query the map with a polygon (click "Draw Polygon" on the right side of the map first). Then you'll see a nice summary and a pie-chart (due to some google-chart-api usage by Josh) of the land-use types in the polygon you queried. The backend of that is agoodle. The map is, of course, openlayers.

You can see the usage on the github page but basically, it'll be something like: which gives you a dictionary of raster grid code keys to the area values-- meaning, for each landclass, it tells you the sum area of cells of that class that fall inside the polygon. Note that summarize_wkt can also take keyword arguments for raster_epsg and poly_epsg, so you can query with a polygon that's in a different projection than the raster.

The library could use some work. There's sketchiness about "inside", so the sum of the areas in the returned dictionary will not match the area of the polygon. Internally, we use matplotlib's nxutils which does a very-fast point in poly test for every cell in the raster that's within the bounding box of the polygon. So any cell that passes that test will be included in the results--all or nothing, there's no corrections for the case where a cell is half in the polygon. This is not a problem unless the polygon is small relative to the size of the raster cells (or the perimeter of the query polygon is large). Patches Welcome.

Check out the code, we made gdal, numpy and nxutils do all the hard work, but it comes together pretty well. If there's other opensource projects out there that do this well, let me know and i'll link to them. I'm posting this because as far as I know, no others exist.

fileindex

2010-04-10T17:27:00.000-07:00

Disclaimer: This is about a generic indexing method that fits in < 50 lines of code, while it's cool, I'm not suggesting it be used in place of a real project. The code is here and explained below.

Everyone's indexing big sequence files. Peter Cock has a bio-python branch to index via sqlite. Brad Chapman writes up a nice approach using tools from bx-python. And there's my own pyfasta for fasta files. This morning I set out to use another one: screed from the pygr guys to index fastq and fasta files via sqlite db. Titus wrote about screed and posted to the biology-in-python mailng list, which is how i originally heard about it.

Screed and the biopython branch use sqlite to get quickly from name to thing--random access. This is a nice approach because sqlite comes with python and it's easy to use and quite fast. Thinking simply about an index, all it really does is get you from some id (e.g. a fasta header or fastq name) to the thing (the fasta sequence or the fastq record).
In the case of flat files, it's a mapping from and id or name to the fseek file-position at the start of a record. Given that, it's possible to make a generic file indexer that creates an index given a file and a function that advances the file pointer (fh) to the next record and returns an id.

So for the case of the FASTQ format, which contains a new record every 4 lines, the first of which is the 'id', the parser could be:


def fastq_next(fh):
    id = fh.readline() # save the id.
    # advance the file handle 3 lines.
    fh.readline(); fh.readline(); fh.readline()
    return id

So regardless of the file format, this is the interface. The function takes a file handle and 1) advances the file position to the start of the next record and 2) returns the id.
Given that, the indexer call looks like:

FileIndex.create('some.fastq', fastq_next)

all that call has to do is repeatedly send a filehandle to fastq_next, accept the the id returned by fastq_next, and save a mapping of id to the (previous) file position. An implementation detail is how that mapping is saved. I use a tokyo-cabinet BTree database.
Once the index is created, usage is:


fi = FileIndex(somefile, handler)
record = fi[somekey]

where handler is a callable (including a class) that takes a pointer and returns a thing. For fastq, it could be:


class FastQEntry(object):
    def __init__(self, fh):
        self.name = fh.readline().rstrip('\r\n')
        self.seq = fh.readline().rstrip('\r\n')
        self.l3 = fh.readline().rstrip('\r\n')
        self.qual = fh.readline().rstrip('\r\n')

so usage for fastq looks like:


>>> fi = FileIndex('some.fastq', FastQEntry)
>>> some_id = '@HWI-EAS105_0001:1:1:7:1680#0/1'
>>> record = fi[some_id]
>>> assert isinstance(record, FastQEntry)
>>> record.seq
TATTTATTGTTATTAGTTATTTTANTANAAATANTNGANGGGGAGGAAGGNNNNNNTNNNNNNNNGANNNNANGAG

FastQ

So, from any name, there's direct access to each record. Given the implementation of FileIndex setup, it's actually possible to create and use an index for a fastq file with 9 total lines of code, 6 of which are the class itself:


class FastQEntry(object):
    def __init__(self, fh):
        self.name = fh.readline().rstrip('\r\n')
        self.seq = fh.readline().rstrip('\r\n')
        self.l3 = fh.readline().rstrip('\r\n')
        self.qual = fh.readline().rstrip('\r\n')

# note the 2nd argument takes a file handle and returns a name.
FileIndex.create('some.fastq', lambda fh: FastQEntry(fh).name) 
fi = FileIndex('some.fastq', FastQEntry)
print fi['@HWI-EAS105_0001:1:1:7:1680#0/1'].seq

SAM

That's decent, but what makes it cooler is that this same interface can be used to implement an index on a SAM (Sequence Alignment/Map) format file:


class SamLine(object):
    def __init__(self, fh):
        line = fh.readline().split("\t") or [None]
        self.name = line[0]
        self.ref_seqid = line[2]
        self.ref_loc = int(line[3])
        # ... other SAM format stuff omitted.

f = 'some.sam'
FileIndex.create(f, lambda fh: SamLine(fh).name, allow_multiple=True)
fi = FileIndex(f, SamLine, allow_multiple=True)
print [(s.name, s.ref_seqid, s.ref_loc) for s in fi['23351265']]

# output: [('23351265', '2', 8524), ('23351265', '3', 14202495)]

That's it. Note one extra thing: In an alignment represented in the SAM file, I used the read id (the first column) as the id for the indexing. That does not have to be unique, so I specify allow_multiple=True and it returns an array of records for every key. (So although using tokyocabinet is an implementation detail, putcat() functions make it much easier to support multiple entries per id.).

Performance

For the SAM file I tested, the original file is 2.5G and the index created by FileIndex is 493M. For a large fastq file of 62.5 million lines (about 15.5 million records) the file is 3.3G and the index is 1.2G (compared to 3.8G for screed). It takes about 5 minutes to index the FASTQ file (compared to 11 for screed).
On my machine, with that FASTQ file, I get about 17,000 queries per second (including object creation). For the same fastq, screed gives about 10,000 queries per second (which matches the image in Titus' post).

Implementation

The full implementation right now is 44 lines, including (sparse) comments and is shown in the gist at the end of this post. It lacks a lot of nice things (like, um tests) and I'm not recommending anyone use it in-place of a well thought out project like the biopython stuff or screed, still, i think the concept is useful--and I am using it for SAM files.
I've put the full code in bio-playground on github, including the examples from this post. But, this is generic enough that it could be used to index anything: blast-files, fasta, ... whatever. As always, I welcome any feedback.

About Tokyo Cabinet

In the course of this, I found major problems in 3 different python wrappers for tokyo-cabinet. I reported bugs to 2 of them. So, kudos to py-tcdb. While I don't like that I have to explicitly tell it _not_ to pickle what I send it, it's well tested, and the code is very nice and ... best of all, I haven't been able to make it segfault. It is also easy_install'able.
Another thing I learned today is that you can't use tokyo-cabinet hash for more than a couple million records. A web-search shows that lots of people have asked about this and the answers always have to do with adjusting the bnum bucket parameter. That does not work. If you can create a tokyo cabinet hash database with over 5 million records that does not slow on inserting, please show me how. Until I see it, I think one has to use the BTree database in TC.

writing and building a lua c extension module

2010-04-04T12:43:00.000-07:00

[Update: 2009-04-09]
This package is now on luarocks under it's new name: "stringy". It now includes startswith and endswith methods (more coming). Under advice of lua gurus, I no longer add to the string base methods. And finally, it's available in my lua github repo.
[/Update]
I've been messing with Lua programming language lately, mostly because I wanted to try to use love2d as a visualization tool. I got side-tracked into building a c extension for lua. The C-API is much different from python.
In lua, all the c functions have a signature like:

int c_function(lua_State *L)

and you use the lua c API to pull values off the lua_state stack thingy -- L.
And then to return values, you just push them back onto the stack. I don't grok this fully yet, but it seems to handle all the memory allocation for you.
Anyway, it's hard to find a full example of creating a C extension for lua 5.1. It actually seems more common just to provide patches for the lua distribution itself. There are some docs but they were difficult for me to find and it's not clear which docs are for lua-5.1--the current version. So I'm including a full example with sourcecode, simple tests, Makefile, and luarock spec.
The full gist is here.
The C code (shamelessly stolen from the wiki and book -- not my own code) actually implements 2 very useful functions string.split and string.strip. These are otherwise easily added using lua's regular expression searches, but these are faster as they're not using the regexp machinery:
Note the functions are included in a struct and "registered" with the luaopen_stringext function.
The Makefile then builds the shared library: The shared library stringext.so is then immediately usable from the lua interpreter. A test session looks like:. Where string.split() returns a table (array) of the tokens. Another cool thing about lua is visible in that script. The added functions actually become methods on Lua strings! So after importing stringext, all strings now have strip() and split() methods! This is because of the line in stringext.c:

luaL_openlib(L, "string", stringext, 0);

which tells it to add the methods to the "string" module.
Finally, luarocks... Luarocks are the equivalent of python eggs or ruby gems (you gotta love all these clever names). They take a rockspec file (equiv of python's setup.{py,cfg}). Mine looks like this.
I've requested that be added to the main luarocks repository, there doesnt seem to be a way to upload your own rock directly. Still once you write the rockspec, you can build the C extension without a Makefile by typing

luarocks make

and it handles all the appropriate flags and such.
Anyone interested can download a tar-ball of the entire distribution here.
I add that the lua community seems very active and helpful. I asked a question about building the extension and received quick, helpful replies.

organizing genomic data

2010-02-16T07:06:00.000-08:00

In some cases, using sqlite, a DBM, or just a python pickle is as much of a pain as dealing with something like gff (generic feature format) which is comprehensive but requires a lot of manipulation. Relational db's are slow and have extra overhead when I don't need the relational or ACID parts and DBM's, for the most part, only allow querying by keys. These are 2 cases where I've found a nice alternative.

Big ol' Matrix

because it was already nicely pre-processed, I downloaded some Arabidopsis expression data from here. They keep the data in a hierarchy of files. I want to be able to quickly find the coexpression for any gene pair. 22573 * 22573 == 509 million entries is a bit large for a hash (aka python dictionary)--and even disregarding speed concerns it's more than I want to push into an sqlite db--but fits in a 1.9 gig file of 32 bit floats. So each entry is a measure of coexpression between 2 genes. In that sense, this is a key-key-value store. [for more details on how the coexpression is calculated, see the link above, I'm using this dataset specifically so I don't have to think about the raw expression data which is difficult to deal with].

I'm using numpy to create and then read that binary file (though it's just 32bit floats so the format will work for any language). The index of a particular gene in the array can be set as it's alphabetic order. So, to find the level of coexpression between the 13656th gene and the 12355th gene is arr[13657,12356]--or using pointer arithmetic if you're not using numpy. To find the average level of coexpression of the 13656th gene with all others is arr[13657].mean(). Looking up the index of each gene is a simple matter of reading in the file of gene names (already ordered alphabetically) and creating a dictionary lookup of gene => index. The one line of code to do that looks like this:


idx_lookup = dict((accn.rstrip(), i) for i, accn in enumerate(open(path + "/names.order")))

which can take an accn name like "At2g23450" and return its integer index.
Since the binary file is 1.9G of data, it's not ideal to read the entire thing into memory. But, it's simple to memory map it with numpy into a big ol' matrix/array:


L = len(idx_lookup)
coexp = np.memmap(path + '/coexp.bin', dtype=np.float32, shape=(L, L), mode='r')

so then, given a pair like: AT1G47770,AT2G26550
the code to find the coexpression is:


ai = idx_lookup["AT1G47770"]
bi = idx_lookup["AT2G26550"]
ce = coexp[ai, bi]

this is stupid fast. To memory map the array, read in the idx_lookup and do 6294 queries takes 0.17 seconds on my machine--and much of that is the one-time cost of reading in idx_lookup. The code for all of this is readable here.
and svn'able via:

svn checkout http://bpbio.googlecode.com/svn/trunk/scripts/coex/

That includes the code to get the data (get.sh) convert it to the .bin file (to_array.py) and serve it as a web script (coexp.wsgi).

Flat

In most eukaryote organisms there are on the order of 30K genes (varying from maybe a couple thousand up to 60K or so). Each gene can have multiple sub-features: CDS coding regions, mRNAs, UTR's, introns, etc. So it's a pretty small amount of data, but large enough and with enough structure that it's nice to have it in an efficient data structure. After five years, I think am finally converging on a good abstraction. It holds 1 feature per row in a subclass of a numpy array. This hides some complexity arising from the fact that each gene feature has multiple exons and each gene can have a different number of exons (gff uses one row per cds * gene * mRNA * exon). Numpy has typed data "fields" which are something like columns. Common types are floats and ints of various precisions and python Objects. So, the exons start, stops are stored as a list of python integers in an object field of a python array and that python array is stored as a 'locs' field in each row of the numpy array. The .flat file format looks like this:

#id  chr accn    start   end    strand  ftype   locs
1   1   AT1G01010   3631    5899    +   CDS 3760,3913,3996,4276,4486,4605,4706,5095,5174,5326,5439,5630

where the columns should be mostly understandable by the header. `start` is the left-most position of the feature on the genome and `end` is the rightmost. `ftype` tells the type of the feature(s) in the following column, `locs` which is a list of start,stops. most often ftype is `CDS` (or exon) for coding sequence--or the part of the gene that gets translated. From there, it's possible to take advantage of the speed and syntax of numpy:


 >>> from flatfeature import Flat
 # note also attaching the genomic fasta sequence.
 >>> flat = Flat('data/thaliana_v8.flat', 'data/thaliana_v8.fasta')
 # so can get coding sequence.
 >>> cds_seq = flat.row_cds_sequence('AT1G01370')
 # find all the gene names in a given chromosome, region:
 >>> [gene['accn'] for gene in flat.get_features_in_region('1', 5000, 7000)]
 ['AT1G01010', 'AT1G01020']

 # how many features on chromosome 4 are on the - strand:
 >>> flat[(flat['seqid'] == '4') & (flat['strand'] == '-')].shape
 (2502,)

The row_cds_sequence, while an ugly name, shows the advantage of tyeing the fasta to the Flat object--it allows extracting the genomic sequence from the fasta based on the coordinates in each row. It also highlights a problem with this data structure--what about alternative splicings which are very common? For our lab, we always use the union of all coding sequences so I report the resulting intervals as the 'locs' column. I intend to improve the flexibility of how that is handled in the code, if not the file.

As a side note, the final example in that session shows how you can think of the numpy slicing syntax as a sort of SQL-like selection. so: flat[(flat['seqid'] == '4') & (flat['strand'] == '-')] reads like: "select * from flat where seqid = '4' and strand = '-'", but there's no database, it's all in memory, and the work is all done in C, very quickly. Actually, I think flatfeature and the huge matrix could both classify as NoSQL.

The code for flatfeature is on github it includes an example .flat file for arabidopsis thaliana in the data directory. I'll be tinkering with that for some time. As with all my code, any forks, patches, suggestions or ridicule will be welcome (in that order of preference). Finally, to add that for many cases genometools works great, they do a nice job of wrapping GFF for a number of languages including python. But for some things GFF is too much trouble.

Loopless programming (for calculating methylation types)

2010-01-18T10:02:00.001-08:00

DNA Methylation is used in plants as a means of epi-genetic regulation. Bi-sulfite sequencing is a method used to determine the methylation pattern of a given set of DNA. Methylation occurs at 'C'ytosines. In plants, there are 3 types of methylation, determined by the nucleotides that follow the 'C':

CG: a 'G' follows the 'C'

CHG: anything but a 'G' follows the 'C' and a 'G' follows that

CHH: no 'G' in the 2 subsequent nucleotides.

So, programmatically, it's a trivial matter to calculate the type of methylation that can occur at each cytosine in a given sequence. Though, once the edges and edge-cases and reverse-complement are handled, it becomes a few lines of code full of loops and if's. For python (and most scripting languages) that's slow and loopy and iffy. The rest of this code explains how I utilized numpy to make a fast methylation-type calculator without loops or ifs.

Loopless

Given a numpy array `seq`, it's possible to find all the 'C's with:

np.where(seq == 'C')

. From there one can calculate the methylation type without looping or if'ing with:

which potentially takes more memory but is extremely fast, and (IMO) nicely shows how to take advantage of the things numpy is good at using the slicing and indexing.
From there, I put that (slightly modified to differentiate between + and - methylation) into a function named _calc_methylation, and call it from this function:

which does some work to handle the ends of the sequence and reverse-complementing. The result, as shown in the docstring, is a numpy array where the values indicate the type of methylation as:


0: not a C or G
1: CG at a 'C' (+ strand )
2: CHG at a 'C' (+ strand )
3: CHH at a 'C' (+ strand )
4: CG at a 'G' (- strand )
5: CHG at a 'G' (- strand )
6: CHH at a 'G' (- strand )

The full implemenation is here. (even the reverse-complement is done without loops or ifs)

needleman-wunsch global sequence alignment -- updates and optimizations to nwalign

2010-01-10T09:06:00.000-08:00

I've written previously about nwalign, a python package I wrote to do fast (thanks to cython) global sequence alignment. Over break, I spent some time fixing some bugs and improving performance.

Ack

It's actually nice to get bug reports for a relatively obscure bit of software like this as it shows it's getting used. Thanks especially to R. Christen for his patience in (repeatedly) showing me places where nwalign was not doing the right thing.

Bugs

Placement of Gaps

Some of the "Bugs" were actually just "unspecified behavior". For example, given input text of "AGEBAMAM" and "AGEBAM", the alignments:

AGEBAMAM
AGEBAM--

and

AGEBAMAM
AGEB--AM

Have the same score. However, the first version is generally more, um, polite. Previously, the behavior was deterministic, but depended on the length and order of the input sequence. As of version 0.3, nwalign will always put matches together when given a choice between 2 (or more) paths through the scoring matrix.

Gap Extension

nwalign allows use of a scoring matrix to lookup the score/penalty for conversion from one letter to another or a simpler version where the user specifies gap_open, gap_extend, mismatch penalties, and match reward. For the simpler non-matrix path, earlier versions of nwalign did not heed the gap_extension--using the gap_open penalty regardless of the previous entries in the dynamic programming matrix. It's common to use only a gap penalty without a separate gap_extend penalty, nwalign no longer assumes that's what is preferred -- if it is, one can simply specify a `gap_extend` penalty that is equal to `gap_open`.

Scoring

Previous versions of nwalign only returned the alignment and didn't provide access to the score of the alignment. Recent versions allow the user to get the score via:


>>> nw.score_alignment('CEELECANTH', '-PELICAN--', gap_open=-5,
...                     gap_extend=-2, matrix='PAM250')
11

where `11` is the score of the alignment given, the gap_open and extend_penalty, and subsitution scores in the file 'PAM250' (which is in the NCBI substitution matrix format).

Performance

The earliest version of nwalign was literally about 100 times faster than the perl version from the BLAST book. But, there were a few places where performance has improved even more since that. The most dramatic was in the lookups for the scoring matrix values. Given, again the 2 short strings: "AGEBAMAM" and "AGEBAM", the alignment algorithm has to do the inner loop 48 ( 6 * 8 ) times, in each of those loops, it has to look up the substitution score in the matrix. For example in the first iteration, it has to find the score for "changing" from an "A" to an "A"; in the PAM250 matrix, that score is +2. In the substitution matrix, the row and column "keys" are the amino acids and the values are the scores for changing from the row key to the column key. Previously, I stored the keys in their own array, then did a linear search to find the index. So, for both the row and column, an amino acid letter is translated to an index with a function like:


aa_keys = ['A', 'R', 'N', 'D', ... ] # or 'ARND...'
def findpos(amino_acid, aa_keys):
    i = 0
    while i < len(aa_keys):
        if aa_keys[i] == amino_acid:
            return i
        i += 1
    return - 1

except in C(ython) so it was much faster, the alignment then uses the return values for row and column amino acid to look up the score from the substitution matrix. When the string are long enough-- as they will be with real proteins, this is a huge speed bottleneck. So, it trades memory for speed. Instead of using a 25 * 25 matrix to store the substitution matrix, it now uses an X * X matrix where 'X' is the ord value of the largest amino acid. Usually this will be less than 'Z', so the matrix will be 90 * 90, stored efficiently as a numpy array. From there, the lookup is directly with matrix[ord(x_amino), ord(y_amino)] which is extremely fast in C because the ord is unnecessary as a char can index an array. This gave at least a 3X speed improvement. I could reduce the memory used by the matrix by subtracting ord('A'), but that would be trading speed for memory since it would require 2 extra subtractions per inner loop.

Also, the latest version uses less memory in other areas; in order to do the dynamic programming, the algo has to save N * M arrays of direction "pointers" (in the left/right/diagonal sense (not the c sense), scores and gaps. However, the gaps are not actually needed for the entire trace-back, they are only needed 1 step back to determine if a current gap is a gap_open, or gap_extend. So, now, instead of an N*M matrix for the gap matrix, it's N*1 and M*1 arrays. For large N * M, this is enough to offset the increased memory incurred by how the substitution matrix is stored.

Running a benchmark test script with nwalign-0.1 takes 12.13 seconds and the same script on nwalign-0.3.1 takes 3.41 seconds (and actually runs in 2.99 seconds with unladen-swallow... but that's another story). The script does an alignment on 1200 * 1600 basepair sequences 100 times. As with most things, I figured most of this by banging my head against the wall long enough that the stars aligned, so to speak, anyone who cares to look at the code and offer suggestions would be much appreciated.

updates to pyfasta

2010-01-04T18:24:00.001-08:00

At the end of last year, I put in quite a few updates to pyfasta. One of the nicest is the new flatten stuff. In order to provide fast access to the sequence data, pyfasta creates a separate flattened version of the sequence file containing no newlines or headers for any file that it interacts with. That flattened file is used as the basis for the index which allows fast random-access. This is an additional file, nearly the same size as the original, and can be more space overhead than one would like to incur when dealing with large files. The new "flatten_inplace" keyword arg to the pyfasta.Fasta() constructor will remove all newlines but keep headers. This will leave the fasta file in a valid FASTA format that BLAST or any sequence tools will understand, but will also allow fast access via pyfasta, since pyfasta only needs to know the file position where each sequence starts and ends.
With this option, a file like:


>a
actg
actg
actg
>b 
aaaacccc
aaaacccc

will overwritten (flattened in-place) to:


>a
actgactgactg
>b
aaaaccccaaaacccc

simple. When dealing with large files (where this is actually useful), the flattened file does not behave well when opened in an editor because the editor will attempt to read a number of lines into the buffer and a single line may be 200 mega-bases. So, this is a pain if you're planning to sit down with a cup of joe and read through the genome, but otherwise, the fasta file should be un-affected.
This method is not currently the default (though it may be so in future versions). But, it's possible to use the commandline:

pyfasta flatten some.fasta

which will create the flattened fasta (and the index file) and a placeholder some.fasta.flat, containing the text "@flattened@" as a marker to pyfasta that it's ok to use the original (now-flattened) fasta. Once the file is flattened, there is no performance loss compared to having a separate flat file containing no headers.

pyfasta was a fun project for me in 2009. It's a ridiculously simple little module, but when I started it, I didn't think there was a good alternative. (Though discriminating Fasta-ers should look at the sequence module in pygr, and the Bio.Seq module in BioPython which I think has improved quite a lot recently). It has over 100 tests and very close to 100% test coverage for the modules in pyfasta, and much of the code is run once for each of the 4 backends.

the source is on bitbucket.

genome scrubber : mask repetitive sequence

2009-12-23T06:46:00.000-08:00

This is to describe a simple tool I've made available (svn repo) for masking repetitive sequence.

rice (Oryza Sativa) version 5 sequence looks something like below when run through pyfasta info.


rice.fasta
==========
>1 length:43596771 
>3 length:36345490 
>2 length:35925388 
>4 length:35244269 
>6 length:31246789 
>5 length:29874162 
>7 length:29688601 
>11 length:28462103 
>8 length:28309179 
>12 length:27497214 
>9 length:23011239 
>10 length:22876596 

372.078M basepairs in 12 sequences

So, it's not huge (still only 1/10th the size of human) but, it can be difficult to deal with the entire genome because of the large amount of repetitive sequence and transposable elements. This is sometimes mistakenly referred to as "junk DNA", while that's not true, it does make whole-genome analyses a pain as a the output is dominated by repetitive sequences matching their own families. Doing a blast of the rice genome with this command:

/usr/bin/blastall -K 100 -a 8 -d rice.fasta -e 0.001 -i rice.fasta -m 8 -o rice_rice.blast -p blastn

nearly causes my machine to run out of memory (16G), takes a couple of days to run, and results in a blast output of 5.1G and 84 million rows--that's 84 million blast hits with an e-value below 0.001! By definition, that output is dominated by the repetitive elements. Repetitive elements are interesting, but in the case were we want to look at synteny, we have to wade through that 5.1G of stuff to find the very small chunk of data we need. This adds time to run the sequence comparison, time to parse, time to plot, time to analyze, and data to store, etc...
The solution we use in our lab is to create a "masked" sequence where we 'X' out any DNA sequence occurring more than a given number of times in the original blast. So, from the output of the blast above, any basepair that is covered more than 50 times is hard-masked with 'X'. The single script to run this is here. Available from svn via:

svn checkout http://bpbio.googlecode.com/svn/trunk/scripts/mask_genome

when run as a script with no options, it will print some help text. To mask rice at 50X, it is run as:

python mask_genome.py -b rice_rice.blast -f rice.fasta -o rice -c 50 -m X

This hard-masks with "X". To soft-mask (changing basepairs to upper-case except those occurring more than 50 times which are lower-cased.) one can use the commandline option -m SOFT.
The resulting file: "rice.masked.50.fasta" has the same number of chromosomes, each with the same number of basepairs, but with repetitive regions masked. To show the difference, I re-ran the full-genome blast, except this time on the masked genome:

/usr/bin/blastall -K 100 -a 8 -d rice.masked.50.fasta -e 0.001 -i rice.masked.50.fasta -m 8 -o ricemasked.blast -p blastn

. The resulting blast file is only 128MB and 2 million rows (contrast to 5.1G and 84 million rows) or a blast file that's 2.4% of the original size. This makes things simpler, and makes doing genomic analyses more enjoyable. This is a different tool than the dust filter that comes with BLAST or many others, because it masks based on the global coverage at each location (there are similar tools).

Implementation

Counting

In order to do the full-genome mask, mask_genome.py stores an array for each chromosome, with a length equal to that of the chromosome in a pytables/hdf5 structure which acts like a numpy array. So to increment the counts for a blast hit, the code looks like:


        cache[qchr][qstart - 1: qstop] += 1
        cache[schr][sstart - 1: sstop] += 1

where cache is the hd5 structure (except for some cases it's read into memory to improve speed), 'qchr' and 'schr' are the query and subject chromosome keys with the values being the count arrays, and 's/qstart', 's/qstop' are the bounds of the subject and query for a particular blast line. This sends much of the work to numpy, instead of iterating over the range of qstart, qstop in python. That is repeated for every row in the blast file.

masking

Once all of the counts are created, the script uses that data to mask out the original sequence for any element in count that is greater than the cutoff. That is achieved in a single line iterated per chromosome:

masked_seq = np.where(numexpr.evaluate("hit_counts > %i" % cutoff)
                                      , mask_value, seq).tostring()

where `np` is from "import numpy as np", `seq` is the original fasta sequence in upper case, `mask_value is either a constant, like 'X', or in the case of soft-masking, a lower-case copy of `seq`, and `hit_counts` is the array of counts. The only tricky bit there is the use of numexpr which makes things a bit faster. The `masked_seq` is written to file after it's corresponding header, and the masking is performed for the next chromosome. The entire masking takes under 10 minutes on my machine.

Plotting

Generally, we used the masked fasta file as the final output, but there's useful information in the hd5f file--a sort of measure of repetitiveness at each basepair in the genome. also in the svn repo is a short matplotlib web script (for mod_wsgi/mod_python but can be modified to run as cgi or command-line) to generate an image given a request like:

?org=rice&xmin=12&xmax=90000&width=800&seqid=2

where 'org' is the name of the organism sent to the mask_genome.py script, `xmin`, `xmax` specify the extents in basepairs, and `seqid` the chromosome.
Here's an example from our genome viewer (built on openlayers) where I've added copy count as a layer (with admittedly poor cartography). The nice gene renderings and basepair location ruler are courtesy of Eric Lyons.
The height of the blue line indicates the copy count. The flat red line is at a copy count of 50. So, the genes in this region are repetitive, but it's possible to see that they are surrounded by long terminal repeats, and that there is some repetitive sequence in the introns (the gray between the big green exons). Scrolling along, openlayers-style, it's possible to see patterns, and spot transposons.

web

I've made a few genomes available and will be adding more. There's also a (currently undocumented) service to plot the genome. So if the organism is brachy and the chromosome of interest is Bd1 then:

http://copymask.syntelog.com/plot/brachy/Bd1/?xmin=12&xmax=38345&width=600

will result in:

where the image is described above.

If you're doing full-genome analyses, give it a try and let me know any suggestions. There's a full example starting from only a genomic fasta, including the BLAST command in the README.rst included in svn. If you get a genome that's too big to run through a full-genome BLAST (maize, ahem). You can split with "pyfasta -split -n 4", run 4 blasts, concat the results, and send through mask_genome.py.

rst2s5 with syntax highlighting

2009-10-23T16:28:00.000-07:00

Restructured Text to S5 Presentation

(lots of caffeine today so 2 posts in one day)
The stub example presentation I'll be talking about is viewable as a presentation here (click on that page to advance the demo slides).

There's a nice browser-based tool for presentations called S5. In recent python docutils, there's a tool called rst2s5.py which converts restructured-text to an s5 presentation. However, it's not obvious how to get syntax highlighting for code blocks to work.
So pygments, a python library that will highlight syntax for many programming languages comes with this file which they recommend you use as a starting point. That's what I did, and I've created a stub example project accessible via subversion:


$ svn export http://bpgeo.googlecode.com/svn/trunk/rst2s5_template/

with a build script and a couple of example slides (and a nice theme). It's possible to change the theme by editing rst-directive.py (included in the source) and changing the STYLE global to a theme that pygments knows about (being comfortable with my machismo, I chose "fruity"). One way to find other themes is to check the drop down box on the paste bin at the pygments site.

To use the template, just edit the index.rst, then run build.sh and view the resulting index.html in a browser. The S5 "theme" it's using is specified in the build.sh script and contained in the ui/ sub directory, you can find more themes on the s5 site and others that come with the s5 install.

Check out the pretty example project here (with python code 3 slides in). As with any S5 slideshow, you can click to advance slides or use the controls that appear in the bottom right when the mouse is in that area.

As a bonus, if the index.rst file contains python shell sessions (doctests) like the example, you can check them with nose using:


$ nosetests --with-doctest --doctest-extension=.rst index.rst

some python ctypes stuff in Rtree

2009-10-22T17:55:00.000-07:00

I've been working with and on the Rtree python module. It's some cool work done by Howard Butler (hobu) (originally with Sean Gillies) to make python bindings for the Spatial Index C++ library by Marios Hadjieleftheriou which provides various tree structures for efficient spatial searches in n-dimensions. Hobu has written a C API for that along with a new ctypes wrapper to that API which appears in Rtree 0.5 and greater. There is some cool ctypes stuff in there which I'm starting to understand.
From the website:

ctypes is a foreign function library for Python. It provides C compatible data types, and allows calling functions in DLLs or shared libraries. It can be used to wrap these libraries in pure Python.

as a simple example of how ctypes works we can pretend there's no atan() in python's math module and access the one from libm in c like this:

import ctypes
libm = ctypes.CDLL('libm.so.6')

# the following 2 lines correspond to the c signature: double atan(double)
libm.atan.argtypes = [ctypes.c_double]
libm.atan.restype = ctypes.c_double

print libm.atan(0.22)
print libm.atan(ctypes.c_double(0.22))

where line 2 tells ctypes how to find the library (have a look at Rtree or shapely source code to see the cross-platform way to do that). lines 4, 5 tell it the input types (argtypes) and return type (restype) respectively, and lines 7, 8 call the c function by way of the ctypes wrapper. Here, it's calling the version of atan with a double precision number. With simple types, you can let ctypes wrap a python value in the type or you can do so explicitly as in the last line.

Things get more interesting with more complicated return types. For c function with a char * return type, e.g. this contrived example:

// does not need to be freed.
char* fn_char(){
    char *s = "asdf";
    return s;
}

the ctypes invocation -- with 'ccode' being this contrived library as loaded by ctypes.CDLL -- looks like:

ccode.fn_char.restype = ctypes.c_char_p
print ccode.fn_char()

which returns "asdf" as expected and does not leak. (Note that ctypes.c_char_p is "c char pointer" or "char *".) If you get a copy of a char * and are responsible for freeing it's memory, e.g. from the c function:

// needs to be freed.
const char * fn_const_char(){
    return (const char * )strdup("asdf");
}

the ctypes for that looks like:


def get_and_free(achar_p): 
    s = ctypes.string_at(achar_p)
    libc.free(achar_p)
    return s


ccode.fn_const_char.restype = get_and_free
print ccode.fn_const_char()

where libc is the standard c library defining the function to free() memory. In this case, it takes advantage of the feature that .restype can be a callable which takes the pointer return from the c code. In get_and_free(), ctypes.string_at() turns that pointer address into a python string. Then the char * pointer is free'd, and the python string is returned, and "asdf" is printed as expected.
It's also possible to do more rigorous error checking with errcheck in which case the ctypes looks like:


def err_check(char_p, fn, args, fn, args):
    s = ctypes.string_at(char_p)
    libc.free(char_p)
    return s

ccode.fn_const_char.restype = ctypes.POINTER(ctypes.c_char)
ccode.fn_const_char.errcheck = err_check

where err_check gets 3 arguments, the result of the function, a reference to the function, and the args sent to it (in this case args is empty). Note that in this case, we have to specify the restype ctypes.POINTER(ctypes.c_char) so that we still have the pointer address--which we then free. When the restype is specified as ctypes.c_char_p (char *), then ctypes automatically gives us the python string and we can't (as far as I know) free the memory and a leak occurs. Also, in the case above, I haven't actually added any extra error checking, in Rtree hobu has a few functions to check the return values, see that code here.

This post has been pretty basic, there's also good ctypes code in shapely, geodjango, and libLAS. Next post I'll talk about callback functions -- calling a C function that expects a pointer-to-a-function with a python function as an argument.

biostuff

2009-10-14T09:14:00.000-07:00

I've been trying in 2009 to write less throw-away code. I'm not sure how successful I've been at that, but at least I'm writing more code that I keep around.
Previously, I stuck anything of at least marginal quality and re-usability into my google code project bpbio. As of yesterday, I've moved a lot of stuff from there to bitbucket. "Biostuff" is where I'll put modules that are well documented and tested in hopes that using a distributed VCS and a project that doesn't contain my initials will foster any contribution. Currently, all the modules on bitbucket are also on pypi.

pyfasta provides pythonic access fasta sequence files. Previously, it had been a part of genedex (which I've stopped supporting since @hobu has done so much good work on Rtree that genedex is now pretty much obsolete) but it's been pulled out and simplified and improved. Check out the docs on pypi.

nwalign is a command-line or python interface to the Needleman-Wunsch global sequence alignment which I've blogged about previously. Whenever I need to do stuff with cython and numpy I use nwalign.pyx for reference (though there's probably better material out there).

simpletable, as the name suggest is a wrapper around pytables to remove some of the boiler-plate in creating a table and dataset and Description. That's it.

skidmarks is a small module to check for runs in data (get it skidmarks, runs?). It implements Wald-Wolfowitz, autocorrelation, serial, and gap tests. Each function (implementing one of those tests) returns a p-value which indicates the level of support to reject the null hypotheses that the sequence is random and the chi-square or z-score value as appropriate. I've been using this and monte-carlo simulations to see if runs in genomic data could be explained by random events.

Any contributions, suggestions, or bug reports are welcomed--the interface at bitbucket should make this easier to do, just fork and fix and pull-request.
Meanwhile, my less documented/ more crappy code will continue to live on google code--at least until it matures. I've got a couple modules in the pipeline that will be added once they're cleaned up and documented.

starting haxe. (stuff i want to remember)

2009-07-26T14:33:00.001-07:00

I've been tinkering with a flash project recently. Actually haxe, so it's only linux tools -- VI, and the command line -- not the GUI interface people normally associate with flash. This post is a summary of how to get started with haxe using only the command line, and a project containing flash stuff I want to remember.
To start, here's a gist of shell commands that will set up haxe on an ubuntu machine. The installers from the haxe website work fine for windows and mac (and I think 32 bit linux).

Haxe has a slightly different syntax from actionscript 3, but for most things it is identical. These docs are very good, and better than the adobe site, I have that page open always when working with haxe. I also grabbed an "actionscript.vim" from the internet somewhere and put it in ~/.vim/syntax/ for syntax highlighting and added this line to my .vimrc:

autocmd BufRead *.hx set filetype=actionscript

Then compilation and code is simply a matter of following this.

If you only look at haxe/flash occassionally, some stuff can be less than obvious. I've started a github project as a testing ground where I can record the stuff I figure out.
Currently it has examples of:

+ call javascript from flash (ExternalInterface.call())
+ call flash from javascript (ExternalInterface.addCallBack())
+ add an image from a (local) url
+ keyboard events
+ interact with bitmap data of image.
+ style a text field.

all in a single .hx file. Most of those are fairly simple, they just require the correct incantation. There's a version of this running on my work machine here to demo this stuff. The possible interactions are mostly listed as instructions in either the HTML or the flash movie in that page.
To build the flash(9) movie, just type 'haxe build.hxml' in the learnflash/ directory. Commenting out the '--no-traces' line in build.hxml will cause any trace() call to be sent to firebug--this is extremely useful for debugging.

I'll probably add more as stumble upon it.

displaying and serving big(ish) data with numpy and memmap

2009-05-23T11:37:00.001-07:00

In this case, "big" is only 8 million rows, but I have it working for 40+ million extremely fast and it will work for anything that can fit in memory (and larger with a bit of extra work). So the idea is to create simple web "service" (where "service" is just the URL scheme I pull out of my .. uh genes ... ) implemented in a wsgi script.

A web request containing a width=XXX in the url will show the user wants an image. so a url like:
http://128.32.8.28/salktome/?seqid=1&xmin=491520&xmax=499712&&width=512
will give an image:

(geo-hackers will recognize that with a URL scheme like that, it'll be simple to put this into
an openlayers slippy map.) where the each horrible color represents a different tissue, and values extending up from the middle represent the top (+) strand (W) while those on the bottom represent the - or Crick strand. The heights are the levels of expression.
Without the width=XXX param, the service returns JSON of values for the + and minus strand for each basepair and tissue in the original file--which is less useful than having a summary of mean value for each tissue over the range:
http://128.32.8.28/salktome/?seqid=1&xmin=491520&xmax=499712&summary=t
gives:


{"+": [0.008325556293129921, 0.0072462377138435841, 0.0074741491116583347, 0.0017609233036637306, 0.010895536281168461], "-": [0.0088664013892412186, 0.017081102356314659, 0.009236418642103672, 0.0043859495781362057, 0.021286465227603912]}

which are the mean values from basepair 491520 to 499712 for the + and - strands and each of the 5 tissues. A single tissue can be requested by number as:
http://128.32.8.28/salktome/?seqid=1&xmin=491520&xmax=499712&summary=t&tissue=2
which gives:


{"+": 0.0074741492195734898, "-": 0.0092364182547917447}

Likewise, one can request and image with only a single tissue:
http://128.32.8.28/salktome/?seqid=1&xmin=491520&xmax=499712&&width=512&tissue=2

Transcriptome data often has enough rows of data that it's possible, but not really convenient to stick it into an RDBMS. the data I'm using is available here and looks like:


26 C 125.05 106.18 119.80 159.78 101.40
26 W 141.90 139.07 137.50 151.91 171.18
51 C 131.85 108.57 71.00 156.41 133.79
51 W 123.83 129.97 122.50 100.23 139.09
76 C 136.84 104.98 120.30 88.88 151.64
76 W 119.35 112.79 160.00 119.51 130.29
126 C 138.93 100.92 111.00 84.82 103.48
126 W 87.82 118.30 110.30 157.17 126.42
151 C 184.50 71.24 157.80 396.25 86.42
151 W 136.76 119.95 107.00 98.10 108.60
176 C 135.75 86.98 115.80 188.41 93.53
176 W 147.57 158.19 131.00 117.45 164.77
[... millions more ...]

where the first column is basepair position, the second is either 'W' or 'C' for the Watson or Crick strand of the (double-stranded) DNA sequence. Columns 2 and on are measurements of transcription (see the wikipedia article) for various tissues in the plant Arabidopsis Thaliana which we use because it has a small genome, it's well annotated and doesn't have too much repetitive DNA or too many transposons.

The data is somewhat irregular in that not every 'C' row has a 'W' counterpart and sometimes vice-versa, for those cases, I create a row to fill the missing row with the same values as the existing row. The code below runs through and creates numpy .npy files of the data:

It saves the positions and the actual transcriptome data into separate files because for the web service, the transcriptome files will be memory-mapped while the smaller position files will be read into memory. This code is run at startup (and so not done for every request):


from numpy.lib import format
for ichr in range(1, 6): # for 5 at chrs
    schr = str(ichr)
    # memmap these
    seqids[schr] = format.open_memmap(os.path.join(path, \
                               'at.tome.data.%s.npy' % (schr,)), mode='r')
    # read these into memory.
    posns[schr] = np.load(os.path.join(path, 'at.tome.posn.%s.npy' % schr))

So that the seqids dict has values that are the memmaped contents of the expression data in the .npy files. For each web request, it then uses numpy's searchsorted to do a binary search so that when a user requests &xmin=1234&xmax=4567 searchsorted() finds the position in the array where 1234 would fall (exactly like python's bisect module). and like wise for 4567. The pair of array indexes from the searchsorted() calls with xmin and xmax give the indexes into the tome array for which to grab the data:


    minidx = posns[seqid].searchsorted(xmin)
    maxidx = posns[seqid].searchsorted(xmax)

    data = seqids[seqid][minidx: maxidx]
    data_idx = posns[seqid][minidx: maxidx]

where data contains the values of the transcriptome data and data_idx contains the associated base_pair positions. I could have saved the basepair positions in the seqids/data numpy array, but since searchsorted() will likely traverse many chunks of the array, I think it's best to have that part in memory rather than memory mapped.
The rest of the code is a bunch of if statements (which I should really spread across multiple functions...) deciding whether to show an image with matplotlib or grab a summary of the data and return it via simplejson. The gist of the entire wsgi script is at the end of this post.
An example of what this looks like in an openlayers map with gene annotations is in the image below.

where one can see some correlation between where the (green) genes are and the higher transcriptome levels. In summary, it's nice to be able to take data in this format and write a 120 line script which provides full, fast access both to the raw data and to images which we can use to find patterns to explore.

Needleman-Wunsch global sequence alignment

2009-04-21T18:09:00.000-07:00

[EDIT 01-01-2010] this is now available at bitbucket.
I've written a simple, fast, python version of Needleman-Wunsch as I couldn't find one to use. It uses Cython and specifically cython-numpy goodness. It's easy-installable as:

 sudo easy_install -UZ http://bpbio.googlecode.com/svn/trunk/nwalign/

or via svn from:

svn co http://bpbio.googlecode.com/svn/trunk/nwalign/

it will put an executable 'nwalign' into /usr/bin/ which when run will give this message:


Usage: 
    nwalign [options] seq1 seq2 
    

Options:
  -h, --help           show this help message and exit
  --gap=GAP            gap extend penalty (must be integer <= 0)
  --gap_init=GAP_INIT  gap start penalty (must be integer <= 0)
  --match=MATCH        match score (must be integer > 0)
  --mismatch=MISMATCH  gap penalty (must be integer < 0)
  --matrix=MATRIX      scoring matrix in ncbi/data/ format,
                       if not specificied, match/mismatch are used

where the matrix is optional but can be the full path specifying the transition score from row to column as here.
If the matrix is specified, match and mismatch are not used. If the matrix is not specified, match and mismatch are used. a typical run looks like:


$ nwalign alphabet alpet
alphabet
alp---et

$ nwalign --matrix /usr/share/ncbi/data/BLOSUM62 EEAEE EEEEG
EEAEE-
EE-EEG

And usage from a python script can be seen in test.py.

I wrote this for a colleague who was using a perl script. This is over 100 times faster than the perl script for long sequences--and the perl script had the BLOSUM62 matrix hard-coded as a hash. There are a couple places that could still be sped up, but I think the improvement will be relatively small. This kind of script is perfect for the numpy-cython mix as I can just access the contents of the 2d array at c-speed without having to do pointer arithmetic. If there's any obvious optimizations I missed, I'd be glad to know them.

python object initialization speed

2009-04-16T21:32:00.000-07:00

On the Cython mailing list, I saw this mentioned for avoiding init overhead, so i wrote up some code to try it. Basically, instead of using an __init__, it uses the PY_NEW macro (which I don't pretend to understand fully).
I ran a benchmark with 5 cases:

PY_NEW macro (still has python overhead for each call to the creator function)

regular python init

python init using __slots__

cython init (cdef'ed class)

batch PY_NEW: calling PY_NEW from inside cython to avoid python call overhead

batch init on cython class

the timings look like this:


PY_NEW on Cython class: 1.160
__init__ on Python class: 30.414
__init__ on Python class with slots: 10.242
__init__ on Cython class 1.185
batch PY_NEW total: 0.855 , interval only: 0.383
batch __init__ on Cython class total 0.998 , interval_only: 0.540

So, the PY_NEW is .383 compared to .540 for using a __init__ on a Cython class, but both are much faster than python. I was surprised that using slots gives a 3x speed improvement over a regular python class. That Cython is faster is no surprise.
Stefan Behnel explains better than I could.

All the code is smashed uncomfortably into this gist.

for kicks, i tried with unladen-swallow. It comes out almost 2x faster on the python times both with and without slots. I didn't use the optimization stuff. Cython even works with unladen-swallow--just have to rebuild the .so--and the timings are the same as Cython-with-CPython.

apache mpm-worker with php on low memory servers

2009-04-12T18:34:00.000-07:00

I'm partly writing this because I think such valuable info is too hard to find, I'd read a lot that if you want php, you have to use mpm-prefork but it's not true!

I've been using slicehost for dev server for almost a year now. Since my 1 year deal is up, I decided to switch to their affiliate mosso. I really like slicehost, but Mosso "cloud servers" seem a good fit for a server that goes through spurts of development and use followed by weeks of non-use. So now, I can keep it as a 256MB instance at about $10/month and update to a larger instance when doing real dev.
I built it today as a 1024MB instance -- installed all my usual stuff, and updated my build script for ubuntu. That's here.
The machine I'm on is extremely fast, normally I set GDAL building and leave, but it finished before I had a chance. After all was built, I resized it to a 256MB server -- that took 12minutes, but my instance was accessible for at least 10 of those.
After that, I log right back in and all works fine. On this dev server, I have to run PHP (at least I dont have to code in it). I've been getting tired of having to use apache's mpm-prefork because ubuntu won't let me have php and mpm-worker. (For more from someone who actually knows what they are talking about, Graham Dumpleton, the author of mod_wsgi, has a good write up here.) So, every apache process takes up a good chunk of the available memory and things are sloooow. Even just running trac, it starts swapping.
I found a good post for using fcgi for php and followed it blindly and all works!
Using worker_mpm, Trac (via mod_wsgi) is so fast I doubled checked to make sure my instance was still at 256MB. After setting ThreadStackSize to 1.5MB trac takes up only 16% of the available memory. The relevant bits of my apache config look like:


<IfModule mpm_worker_module>
    StartServers          2
    MaxClients          40
    MinSpareThreads      5
    MaxSpareThreads      20
    ThreadsPerChild      20
    MaxRequestsPerChild   5000
</IfModule>

MaxKeepAliveRequests 1000
Timeout 30

# i keep this fairly high for serving image tiles
# from mapserver
KeepAliveTimeout 40

# default is 8 per child process, only use 1.5MB
ThreadStackSize 1500000

AddHandler fcgid-script .php
FCGIWrapper /usr/lib/cgi-bin/php5 .php

and then in 000-default, i added +ExecCGI to the DocumentRoot Directory Options.
If any of that config is not sane, let me know.

I should say that I have no idea how PHP will do under fast-cgi, if it uses more memory, or not, or if it will have any problems, on quick look, everything seems ok.

Apparently, you can also use mpm-worker with php if you build php from source. Since using fcgi worked so easily, I didn't try that.

biohash

2009-02-20T09:51:00.000-08:00

This is a quick project I did a while back but, I've seen people interested in similar ideas, so I'll post my implementation here.

Geohash encodes latitude/longtitude locations into a string such that "nearby places will often present similar prefixes" and the longer the string, the greater the precision. Using this python implementation by Schuyler as a reference, I ported the concept to a "biohash" which can encode intervals. It works in a similar fashion, starting with the extremes and halving the space until it finds the smallest space that contains the interval.

The use to allow efficient search of intervals using a BTree index, as in any relational db. It's implemented with only a dumps() and loads() function after the pickle interface. The dumps function takes start and end args and returns a 1/0 encoded string. The loads takes a 1/0 encoded string and returns the tightest interval it can given that string. Both functions take a rng kwarg, which can be as small as the maximum end value. If all the intervals are small, and the rng is very large, the biohash will not help much. The rng used to load must be the same as the one used to dump or the values won't be correct.

I had plans to finish up a set of SQLAlchemy models for this that would save the hash and use it to do range queries behind the scenes, but haven't finished that up yet. The code is in my google SVN. It will even pull in a fast cython version of the encoder.
If anyone wants to use it, improve it or get it going with SQLAlchemy, it available from SVN with:
svn checkout https://bpbio.googlecode.com/svn/trunk/ihash/
it has tests in __init__.py

Stupid things I did as a Bioinformatics Programmer in 2008

2009-01-09T08:00:00.000-08:00

In 2008, I was good enough at programming to get my ass kicked by hard problems. I think that's the most positive way to say it.
My main bioinformatics project was a long annotation pipeline. It takes days to run, often using 8 CPUs. It's driven by a big ol' Makefile. I made the mistake of passing data between steps in un-structured text files or python pickles.

I'd create one at the beginning of the pipeline and not notice it was messed up in a way that affected other parts until the entire pipeline was done, days later. Toward the end of the pipeline, I'd need something simple, like the strand of a BLAST hit, but I'd have to parse through an entire GFF file, or load some huge pickle into memory just to get to that. Then I'd need some annotation, and I'd have to add a slow step of doing a lookup in a script that'd otherwise run very quickly.

I was passing around data in arrays and tuples, so then when I changed the order or added another element in the script that was creating the tuple, a downstream script that was using the tuple would be using the wrong index to access data. If I was lucky, my code would fail, if not, it'd be using the strand when what it should have had was the start location.

I hit problems where I'd run out of memory. At one point, I ran out of disk space (it's a big series of datasets), hit bugs of software I was using.

When I should have run just 1 chromosome to test the pipleline in 1/100th (comparing 10 chromsomes to 10 chromosomes) the time, I ran it over an entire
genome.

When I should have taken the time to really fix small mistakes as I found them, I instead worked around them, making the code unnecessarily complex as a result. If I had fixed instead of fudged in those cases, I would have been more productive.

I did write tests, but not enough, and I didn't set up the project in a way that it was really testable. I'm still learning how to do that. All the other stuff may just be discipline, but the testing is very difficult for me in bioinformatics. I've extracted what I could into tested libraries and added checks for the intermediate data, and every time I found a dumb error, I'd add an assert. (there's a discussion of pretty much exactly the problems I'm describing here: http://ivory.idyll.org/blog/sep-08/the-future-of-bioinformatics-part-1b.html )

I'd assume that it was ok to use a tuple, because I was only going to store start and stop, but then later I'd need to add chromosome, then strand, then score, and pretty soon i'd have code elsewhere like:

if cns[2][4] < cns[2][5]: 
     ...

where i have no idea what i'm comparing there. That one's simple to fix, just use objects, or dicts at the very least, but a 2-tuple is so tempting, and
a 3-tuple not so bad, and ...

On top of all this, the project was changing as I was working on it, so I was changing what went in to the start of the pipeline, what we were getting out, and the steps I was doing. And because of all the extra little hacks in there, I would be stuck in some function far away from the data that I needed (*).

It sucked. I can make those the sort of mistakes on a project that runs in 5 seconds and has a very simple, and relatively known output. But, not for complex pipelines.

It's been painful watching myself do such stupid stuff, and then reading the code afterward, but I think I've learned a lot. Much of that code still sucks, but I've moved to more rigid data-structures. Every time I make a change now, I do it for real, it's not just hacked in. I have to deal with someone pacing around, waiting for the results and asking me questions like "why is it so hard to just add in the RNA?" or whatever. Also, statistically speaking, I'm
probably running out of mistakes I can make...

And actually, the most difficult things have been inter-personal relationships, but that's for another post--as are the more positive, awesome things I did with projects I set up correctly, and had good test coverage...
And finally, I did make a lot of mistakes, but this has been genuinely a difficult problem.

* see this post, especially the last 4 paragraphs. It's good to hear that someone with 43 years of programming language has the same problems as I do.

landsummary.com

2008-12-17T19:31:00.000-08:00

One of the things I like the least about my real job, and much of the contract work that I do is that i'm usually the only programmer working on each task. So, it's been very fun to work on a project with Josh Livni (His writeup). We got together one afternoon, and by the time we left, we had a reasonable start of what we call landsummary, we've since put in a fair bit of work sprinkled here and there. Josh set up an AWS server--it's nice for me to have fewer sys-admin duties too.
What's actually on display is fairly modest. What it does is takes a user-drawn square, circle, or arbitrary polygon, and uses that to summarize the NLCD dataset along with some census data. The things that make it more than lame are that it's very fast, it can easily be extended to summarize any raster dataset, and we have a sorta cool API (not documented) which allows us or anyone to query the data with a WKT Polygon and request a particular service -- currently nlcd, population, weather.precip, and a couple of services with environmental engineering application. All of them use the same libraries--postGIS for doing the census related stuff, and GDAL, gdal_array for doing the raster (currently just NLCD) queries. Josh handled all the census data, I know very little about that, except that whenever I've tried previously, it's been a pain to work with, and now, Josh has a nice set up for it. I took the lead on the raster summaries. For that, I wrote a little library that wraps gdal_array, so you can take a GDAL datasource:


>>> g = AGoodle("something.tif")
>>> a = g.read_array_bbox([xmin, ymin, xmax, ymax])

And then 'a' is a numpy array with all the niceties that entails. So, if we want to get just the food cells, which have values of 81 and 82 in the NLCD dataset, it's just:


>>> a[(a == 81) | (a == 82)]

For arbitrary polygons, we use a surprisingly fast function from matplotlib to mask anything that's not inside a list of verticies (polygon). Then do any summary stats on the masked array.

Thanks to Josh, we have a fairly nice django project structure, with separate apps for each little analysis we've added. In my previous django projects, I've dumped everything into a single app and hacked away, the structure we have now makes it easier to keep what's needed in my brain. Also, when hacking with someone, I'm less likely to put in total crap code. Josh has already had a good laugh at some code where I found 25 closest weather stations using:


SELECT * from stations ORDER BY ABS(lat - ?) + ABS(lon - ?) LIMIT 25

then sorted those 25 using geopy.distance to make sure it was the real distance. In my defence, I really wanted to use vanilla sqlite and so didn't have postGIS at my disposal--also, it was quite fast for only 6,000 stations. We've since dumped it all into postGIS. There's probably a couple of other gems in there.

So, back to the modest functionality part... Actually, it turns out, this is a fairly difficult thing to do in McClick software--time consuming in user and processor time. So, having a way to click a point and see land-use stats and population data appear in about a second is pretty cool--and it's on the web. We've already found a couple folks with interesting applications, and we're interested in finding more--the original motivation was 'foodmiles'-- from this post. And there's a couple things we'll probably add in from that, people I happened to hear talking in a cafe today were talking about foodmiles and seemed interested in incorporating the carbon foot-print of exporting / importing food vs. growing locally. My friend, Megan also has lots other ideas for things that firms commonly do McClick style with the NLCD data.
There's more info on the about page, but suffice to say we make full use of all the usual open-source GIS, science tools.

python interval tree

2008-11-03T18:46:00.000-08:00

EDIT: added a couple points inline.

I'm obsessed with trees lately -- of the CS variety, not the plant variety. Although we are studying poplar, so I'll be using trees to study trees.
I'd tried a couple times to implement an interval tree from scratch following the Wikipedia entry.
Today I more or less did that in python. It's the simplest possible form. There's no insertion (though that's trivial to add), it just takes a list of 'things' with start and stop attributes and creates a tree with a .find() method.
The wikipedia entry that baffled me was about storing 2 copies of each node's intervals--one sorted by start and the other by stop. I didn't do that as I think in most cases it won't improve lookup time. I figure if you have 1 million elements and a tree of depth 16, then you have on average 15 intervals per node (actually fewer since there are the non-leaf nodes). So I just brute force each of those nodes and move to the next. I think that increases the worst-case, but makes no difference in actual search time--with the benefit of halving storage.

EDIT: now the version in my repo keeps the intervals sorted by start, so it can avoid doing the brute for search at each node during a search when search.stop < node.intervals[0].start. This did improve performance.

The tree class takes a list of intervals and calculates a center point. From there it partitions them into left, overlapping, and right in terms of their relation to the center point. Overlapping are assigned to the current node, and left and right are recursively partitioned in that fashion until there are only `minbucket` intervals per node, or the specified `depth` has been reached AND there are fewer intervals than `maxbucket`. So a tree can have a greater `depth` than requested if it would otherwise have more than `maxbucket` intervals in a single node. The Wikipedia version doesn't have maxbucket or minbucket...

EDIT: the maxbucket actually only works on leaf-nodes, and has no effect otherwise.

I'm sure that's painfully obvious for anyone who's ever taken a CS course, but it was foggy at best for me until I implemented. Below is the entire implementation:


class IntervalTree(object):
    __slots__ = ('intervals', 'left', 'right', 'center')

    def __init__(self, intervals, depth=16, minbucket=96, _extent=None, maxbucket=4096):
   
        depth -= 1
        if (depth == 0 or len(intervals) < minbucket) and len(intervals) > maxbucket:
            self.intervals = intervals
            self.left = self.right = None
            return 

        left, right = _extent or \
               (min(i.start for i in intervals), max(i.stop for i in intervals))
        center = (left + right) / 2.0

        
        self.intervals = []
        lefts, rights  = [], []
        

        for interval in intervals:
            if interval.stop < center:
                lefts.append(interval)
            elif interval.start > center:
                rights.append(interval)
            else: # overlapping.
                self.intervals.append(interval)
                
        self.left   = lefts  and IntervalTree(lefts,  depth, minbucket, (left,  center)) or None
        self.right  = rights and IntervalTree(rights, depth, minbucket, (center, right)) or None
        self.center = center
 
 
    def find(self, start, stop):
        """find all elements between (or overlapping) start and stop"""
        overlapping = [i for i in self.intervals if i.stop >= start 
                                              and i.start <= stop]

        if self.left and start <= self.center:
            overlapping += self.left.find(start, stop)

        if self.right and stop >= self.center:
            overlapping += self.right.find(start, stop)

        return overlapping

Only 45 lines of code. I had added a couple extra attributes so that searching could do fewer checks, but it only improved performance by ~15% and I liked the simplicity. One way to improve the search speed, and the distribution on skewed data would be to sort the intervals at the top node, so they'd then be sorted for all other nodes. Then instead of using center = (left + right)/2, It'd could use the center point of the center interval at each node. That would also allow short-circuiting the brute-force search at the top of the find method with something like:


if not (start > self.intervals[-1].stop and stop < self.intervals[0].start):
    overlapping = [ .. list comprehension ]

But all told, that adds 5 or so lines of code. Oh, and depending on how it's used, it's between 15 and 25 times faster than brute-force search.

EDIT: I added the above check, but it can only do the 2nd comparison "stop < self.intervals.start as the first is invalid given a very long interval. Regarding speed, the smaller the search window, the better the performance improvement. The code is now > 20 times as fast brute force for a very (speaking in terms of looking for genomic features) large swath of 100K. with a search space of 50K, it's 50+ times as fast as linear search.

The full code (including a docstring with homer simpson quote) is in my google code repo. If I've made obvious mistakes or you have improvements, I'd be glad to know them.

twill with XHTML (not viewing HTML)

2008-10-25T17:09:00.000-07:00

Since I couldn't find this anywhere, I'll add it here for those who have the same problem:

I was trying to test a website with twill and got this at the end of my traceback:


    raise BrowserStateError("not viewing HTML")
BrowserStateError: not viewing HTML

After spending a bunch of time making sure that, yes, it was spitting out HTML, I figured out that it specifically means that twill (actually mechanize) doesnt like XHTML.

You can likely fix it by adding this at the top of the script:


b = twill.get_browser()
b._browser._factory.is_html = True
twill.browser = b

Presumably, there's a real reason that check is in place, but works-4-me...

appengine memcache memoize decorator

2008-10-24T16:47:00.000-07:00

[NOTE: see the 2nd comment below about using a tuple as a key. better to just use pickle.dumps]
I've been playing with google appengine lately. I'm working on a fun, pointless side project. Here's what I came up with for a cache decorator that pulls from memcache based on the args, kwargs and function name if no explicit key is given. The code for creating a key from those is from the recipe linked in the docstring.

"""
a decorator to use memcache on google appengine.
optional arguments:
  `key`: the key to use for the memcache store
  `time`: the time to expiry sent to memcache

if no key is given, the function name, args, and kwargs are
used to create a unique key so that the same function can return
different results when called with different arguments (as
expected).

usage:
NOTE: actual usage is simpler as:
@gaecache()
def some_function():
...

but doctest doesnt seem to like that.

    >>> import time

    >>> def slow_fn():
    ...    time.sleep(1.1)
    ...    return 2 * 2
    >>> slow_fn = gaecache()(slow_fn)

this run take over a second.
    >>> t = time.time()
    >>> slow_fn(), time.time() - t > 1
    (4, True)

this grab from cache in under .01 seconds
    >>> t = time.time()
    >>> slow_fn(), time.time() - t < .01
    (4, True)

modified from
http://code.activestate.com/recipes/466320/
and
http://code.activestate.com/recipes/325905/
"""

from google.appengine.api import memcache
import logging
import pickle

class gaecache(object):
    """
    memoize decorator to use memcache with a timeout and an optional key.
    if no key is given, the func_name, args, kwargs are used to create a key. 
    """
    def __init__(self, time=3600, key=None):
        self.time = time
        self.key  = key

    def __call__(self, f):
        def func(*args, **kwargs):
            if self.key is None:
                t = (f.func_name, args, kwargs.items())
                try:
                    hash(t)
                    key = t
                except TypeError:
                    try:
                        key = pickle.dumps(t)
                    except pickle.PicklingError:
                        logging.warn("cache FAIL:%s, %s", args, kwargs)
                        return f(*args, **kwargs)
            else:
                key = self.key

            data = memcache.get(key) 
            if data is not None: 
                logging.info("cache HIT: key:%s, args:%s, kwargs:%s", key, args, kwargs)
                return data

            logging.warn("cache MISS: key:%s, args:%s, kwargs:%s", key, args, kwargs)
            data = f(*args, **kwargs)
            memcache.set(key, data, self.time) 
            return data

        func.func_name = f.func_name
        return func

genedex.fasta with numpy.memmap

2008-10-03T10:10:00.000-07:00

EDIT: added job posting to comments.
I've been working a bit on genedex, I'm still not happy with the way it stores features. Which is a huge pickle of dictionaries where every dictionary is a 'feature' that looks like: {'name':'At2g26540', 'start': 1234, 'stop': 3456, 'strand': 1, 'chr': 2}. So the only way to do a search is by location--and that is _very_ fast, thanks to rtree, but there's no way to search by name or any other attribute--and an entire organism is loaded into memory at once--that part actually works out ok, but it feels dirty. I quickly wrote an SQLAlchemy backed interface to a simple db schema do allow this sort of searching here: http://code.google.com/p/genedex/source/browse/trunk/genedex/models/sqla.py. That already supports Feature.upstream(), downstream(), etc. methods, but it will work nicely once python supports sqlite rtree without any extra work--for now, it just uses BTree indexes on the start and stop. I could use rtree to index the sqlite database, but I'd like to move away from the LGPL. Maybe this KDTree that's already in a scipy branch with a more permissive license. Then it could do both spatial, and attribute queries...
That's all tinkering...

I also cleaned up the genedex.fasta module. The useage is nice, if not entirely the implementation. A fasta file can look like:


>chr1
ATGTCGTCGGCCGC
GGGCCAAGA
CAACGGAGA

>chr3
ATGGAGGAGGCTGGCGAGCGG

>chr2
ATGGCGTGC
ACGGCGGCG
CGCATGTT
CGCCT

where a line starting with > is the name of the sequence ('chr1') and everything up until the next > is the sequence. The problem is the newlines, so every time you want to look at chr1 basepairs 10 to 20, you have to find where the sequence starts, and account for newlines. -- Acutally one should never do that, as Biopython, and pretty much any library will take care of that for you. Pygr for example, creates a new file something.fasta.pureseq which removes all newlines and labels and indexes where the starts and stop are. genedex.fasta.Fasta now does something similar, here's example useage on the file above ('123.fasta').


from genedex import Fasta
f = Fasta('123.fasta')
print f.keys()
print f['chr1'][9:20].tostring()
print f['chr1'][9:20]

after, the fasta file looks like this:


>chr1
ATGTCGTCGGCCGCGGGCCAAGACAACGGAGA
>chr3
ATGGAGGAGGCTGGCGAGCGG
>chr2
ATGGCGTGCACGGCGGCGCGCATGTTCGCCT

with spaces removed--so it's still a valid fasta file (you can also send an argument to the constructor and it will create a new file) and there's a new file 123.fasta.gdx that is a python pickle containing:
{'chr3': (45L, 67L), 'chr2': (73L, 105L), 'chr1': (6L, 39L)}
which indicate the start and stop positions of the sequence in the file.
So the file remains a valid fasta file, but now it can be efficiently sliced. For now, it actually uses a numpy mmmap (numpy.memmap), to take advantage of the broadcasting, other than that, it'd be simpler to just use python mmap. So, when it sees f['chr1'][10:20] it acts just like it's indexing a numpy array, but it's accessing the data directly from the disk (well, not actually, but mmap works it's magic and I dont have to think about that). I like that--I can keep my fasta files as valid, add only a small python pickle file, and get simple, fast, pythonic indexing into them. It does take about 12 seconds to index and flatten the entire sorghum genome (629MB, ~660 million basepairs) on first view, after that first time, it's instantaneous.
Anyway, the source is available and easy_install-able:


svn checkout http://genedex.googlecode.com/svn/trunk/ genedex

As always, I'll gladly take any improvements, bugs, enhancements.

Also, our lab at UCB is looking for another (full-time or close to it, on-site) programmer who knows some biology, perl, and hopefully some python. If you're interested, or know anyone, send me an email. I have no real authority in the matter (or any matter) but I will have some say in this. I'd like to work with someone I can learn from. I'll add a link to the job posting in the comments below once it's posted.

choosing django

2008-08-06T12:16:00.000-07:00

I prefer sqlalchemy and genshi (or mako) and was therefore looking at using turbogears, but I saw a demo of the django admin, and that sold me. Certainly the templating language did not. Before this, I'd only used web.py in my projects. These are the things I've liked/noted:

The first and most important: community. (Oddly enough, as I write this there are 666 projects tagged as django. 'turbogears', 'tg2', 'tg' give less than 50 projects combined. Think someone might have already written what you need? yep.
Also, a great site: http://www.djangosnippets.org/, where I've learned a lot just by reading, and saved myself a lot of time, by extending ideas there.
And the development is active.

Second, django.contrib.*

User authentication is simple, and check google-code for various alterations on the theme.

admin. This was what first made me decide on django. And now, new-forms admin is in trunk. This gives you a pretty nice CRUD interface for models in your app. In the app I've been working on, we have row _and_ field level permissions. We also need to let users edit the fields of other users stored in the database--but only certain fields, with which fields depending on the user viewing and the user being edited. This bit is more hacky than it should be, but quite simple.
My biggest gripe about the admin is that it's too complicated to use custom widgets or validation. Hopefully, this will change.
Oh, and it's too hard to have read_only privileges. (Yeah, I know the admin mantra).

GIS: nuf said.

i18n, t9n. Anywhere there's some text to be displayed in the app, I wrap it in ugettext_lazy (aliased to _), and later, dump it and send it to someone who knows Chinese. When it returns, I make messages, and the app will show in English or Chinese depending on browser preferences. Simple.

docs/help: yeah, the docs may have trouble keeping up, especially with recent rate of change, but it's easy enough to find what you need, and the django official docs are pretty nice. And if you can handle the fact that most responses you'll get on #django will begin with "of course, ...", then it's a great place to get help.

ModelForms: I've just started using these, but, they've already saved me a lot of code.

There's a lot of other nice things, and a lot of django that I don't even know. I'd still consider myself a newb, but it's still possible to get sh*t done.

pylab matplotlib imagemap

2008-06-20T16:51:00.000-07:00

UPDATE 7-10-08:
+ add example for scatter plot
+ link to ken-ichi
===
Figuring how to make a client side image map from a matplotlib image has stumped me more than once. Andrew Dalke does have a working example. Below, I have the minimal example.

It's simple once you get the steps right: just use mpl's transform() to convert the data into the image's coordinate system. Then flip the y-axis as required by the imagemap, then do the normal imagemap stuff and save the html. The only real gotcha, is to make sure to put the dpi in the call to savefig().


import pylab
import sys
import random

name = 'imap'

# make some fake data
xs = range(15)
ys = [random.choice(xs) for i in range(len(xs))]

xys = zip(xs, ys)

# can also use : f = pylab.subplot(121)
f, = pylab.plot(xs, ys, 'ro')
dpi = f.figure.get_dpi()
height = f.figure.get_figheight() * dpi

# convert the x,y coords into image coords.
transform = f.get_transform()
icoords = transform.transform(xys)


# the minimal 'template' to generate an image map.
tmpl = """
<html><head></head><body>
<img src="%s.png" usemap="#points" border="0">
<map name="points">%s</map>
</body></html>"""

# change this as needed, e.g. if not plotting points.
fmt = "<area shape='circle' coords='%i,%i,2' href='http://example.com/%i/%i' >"

# need to do height - y for the image-map
fmts = [fmt % (ix, height - iy, x, y) for (ix, iy), (x, y) in zip(icoords, xys) ]

# NOTE, this dpi is needed!
pylab.savefig('imap' + '.png', dpi=dpi)
print >> open(name + ".html", 'w'), tmpl % (name, "\n".join(fmts))

UPDATE:
When trying to figure how to do this for a pylab.scatter plot, I found Ken-ichi had also done this for a scatter plot.
As of a matplotlib trunk revision 5711, the transform does not get set when the scatter plot is drawn. To set it, I had to set the transform explicitly:


   s = pylab.scatter(xs, ys)
   s.set_transform(s.axes.transData)
   transformed_xys = s.get_transform().transform(zip(xs,ys))

binary search over intervals

2008-06-04T22:42:00.000-07:00

[EDIT: update location of code repo]

This isn't particularly advanced or clever, it's just a simple implementation--less code than anything else I could come up with.

Binary search is easy. Just look at the python std library implementation (and the C API version). When you play the game with a friend of guessing a number between 0 and 100, you guess 50, your friend tells you "higher", you guess 75. That's pretty much binary search. It takes about 7 guesses max to guess a number between 0 and 100. It just requires that the numbers be in order.
Interval search is more difficult. It's not just looking for a single value, but rather for a range of values that overlap a given interval. Also, you can't sort on just start, because some intervals are longer than others, so when sorting by start, the stops may be out of order. So, you have to arrange some protocol. In all the examples I've seen, including the explanation here, that means storing not only the start, but the stop and often the center of every interval--and using them to do the search. That makes things considerably more complicated than binary search. I started with an implementation of interval search here, but couldn't figure out how to customize.

Binary search is kind of a special case of binary search where intervals are of exactly length 0, e.g. start == stop. So, if all intervals are of exactly length 2? Well, then you can sort by start, find the left-most index by looking for start - 2 and find the right most index by searching for the (highest) index of start. That will give you the indicies in the sorted array of all intervals that overlap the query. The highest and lowest correspond to python's bisect.bisect_left, and bisect.bisect_right respectively.
This carries to any length. But, if all the intervals are different lengths. Well, then you can save the longest length, and then for any search, it's:
p_overlaps = intervals[bisect_left(start - max_len):bisect_right(stop)]
but that only gives the putative overlapping. Since we extended the left by max_len, and we may have found an interval whose length was < max_len (meaning its stop is before the start of the query) we have to explicitly test for distance:
real_overlaps = [iv for iv in p_overlaps if distance(query_interval, iv) == 0]
which gives only the intervals that overlap. So, that part is the price to pay for the simplified search. Another way, as suggested in the wikipedia article is to store the length of each interval as part of the interval.

In this setup, the worst case scenario is when a single looooong interval covers the entire range of the list of intervals. Then every search is linear, brute-force search. However, my use for this is genomic data. There, I'll have an entire range of say... 50 million, and the intervals (genes) are rarely longer than 4000 in length (basepairs). So, it's useful to optimize simplicity.

My cython version of this is in my googlecode repo here. It has all the stuff I use, methods for left(), right(), upstream(), downstream(), nearest_neighbors(). Most of the searching work is in the binsearch_* functions -- I couldn't use the python ones. There are a couple of hacks in there:
1) because pyrex/cython don't support closures
2) the left() method is confusing because the intervals are sorted by start, and the left() has to find the nearest intervals by stop(). That's where complexity is increased because of this setup.

On my machine, it creates a tree with 6815 intervals in .016 seconds and does 50000 searches on that tree in .14 seconds. It seems to scale well with the number of features as 50K searches on 68150 features takes .50 seconds. Adding an interval that covers the entire range of all other features (results in worst-case linear search) makes the 50K searches (on 6185 intervals) take 1.54 seconds--which is only so good because the brute-force in method is pretty close to c-speed. This could be optimized by saving the stops in a separate array, or by saving long intervals in a separate array, but it's rare enough, and the code is simple enough as-is, that I'll probably leave it for now.

The "proper" way to do this is with an interval/segment tree, of which there's a very readable version in bx-python. If I'd found that earlier, I probably wouldn't have coded this... The tree is faster for larger number of intervals, but that's rarely going to be an issue, it does take much less memory...

wherecamp

2008-05-21T18:44:00.000-07:00

I agree with every report I've seen. Wherecamp was awesome. I've been telling people, and I'm still sure it's true, that I met exactly zero people who I'd physically seen before. Ordinarily, I avoid meetings, but this is a good format and seems to attract good people. It's fun to meet and work with people who are really into what they do. The talks are less "talky" and more like chat sessions--which is possible when the groups are small.

There was also plenty of time to hack, which was the original reason I went. During and after, I learned some simple things which I'm trying to incorporate into my usual workflow:
In the shell, background a job with "ctrl + z" then get back to it with %i where i is the number shown in the output from "jobs". That's a trick from jlivni.

From crschmidt, I added:
alias doctest="nosetests --with-doctest --doctest-extension=.txt"
to my .bash_aliases. Which let's me do:
doctest tests/
or
doctest tests/test_somefile.txt
to run my doctests instead of python -c "import doctest;doctest.testfile('...')"

And springmeyer showed me a ton of django and geodjango. The admin stuff is just ... nice -- it's how making a db front end should be. I still don't know how to learn that stuff on my own, it seems a lot of it, you just have to know which modules to import and the django book doesn't cover newforms or the new admin stuff as far as I can tell.

Quotes

On friday night, we met up in SF to do some hacking, the never went down, as I couldnt get wireless and it turned into more of a real bar trip. We, were however, talking about python. At one point, it was sorta quiet and out of the silence, comes:
"Python sucks"
from a true lisp hacker in the next booth--complete with curly grey beard and spectacles. He actually turned out to be a cool guy, I think maybe he even admitted that if he couldn't use lisp, python was a reasonable choice--I think that's about as much as you can expect from a lisper.

From crschmidt:
"I don't really know python that well"
Then who was it that basically rewrote featureserver between the hours of 2AM and 9AM when everyone else was sleeping?

slicehost, trac, wherecamp

2008-05-14T18:46:00.000-07:00

I have a development "server" here beside me. It's actually a budget laptop that sold for $799 2 years ago. It's a xubuntu machine, hosting a trac instance, a development server for mapping stuff, postgresql/postgis, mapserver, mysql, and couple svn repos, anything I do for contracting, etc. Oh, and it's also hosting a couple of sites for the multi-national company that my gf works for! all of their servers are windows machines (long rant suppressed).
It used to get warm, so propped it up on 4 tuna cans, 1 for each corner, now it stays cooler. Yep, it's a sweet setup.
Anyway, I pay AT&T or SBC --or whatever they are now called-- for static IP's and a supposedly faster internet connection. My 1 year contract for that is nearly up, so I'm switching to slicehost.

I'm not a sys-admin, I sorta do that for 4 gentoo (not my choice) machines at $work, and my strategy is to set up, rsnapshot and never, ever emerge -u world. ever. So far, it's mostly working. When I have the choice, I use (x/k)ubuntu, I don't care if they do magical stuff (or even if I had to redo my ssh keys today), it just works.

Anyway, I want something easy and idiot proof. I'd heard the hype about slicehost when it came out, and figured it was just that, hype. It's not. I've never used shared hosting before, but this was pretty simple. From entering my payment to ssh'ing as root into my slice took ~ 2 minutes. /proc/cpuinfo shows it's a machine with 4 opteron dual-cores. I started with a 256 slice with back-up. You can start new slices and restore from the backups. That's cool, and less $$ than I pay AT&T for static IPs and faster uploads.

I ran a script to apt-get all the packages I use, and had a base, working system in under 1 hour. They have a lot of articles about how to set stuff up, mostly basic (even for me), but I followed their info on setting up iptables. Predictably, I forgot to leave open port 22 and locked myself out of ssh, but they have a web-based console, so, not a problem. It seems to be idiot proof...

Trac

Also, in the theme of things that just work as they should, upgrading to Trac 0.11 (still in rc).


sudo easy_install -UZ Trac==0.11rc1
cd /path/to/trac/project/
sudo trac-admin . upgrade
sudo trac-admin . upgrade wiki 
sudo apache2ctl restart

WhereCamp

I'll likely be at wherecamp, it'll be good to learn some stuff, and meet people I only know from IRC. If anyone needs a ride from the east-bay, let me know.

openlayers, genomes and image-maps

2008-05-08T21:46:00.000-07:00

In response to Titus' post on using imagemaps for genomic visualization:
Why are imagemaps so popular in genomics? As an extreme and unfair comparison, just imagine if http://maps.google.com was an image map.
Given a CGI script that can accept a url like

&start=1024&stop=2048&chr=3

and return an appropriate image, you can provide a substantial set of tools using openlayers, which is developed by what must be one of the largest and active developer communities in GIS. (Yes, I am an openlayers fan-boy.)
You can do that with a small addition to openlayers which I updated a couple weeks ago to OL version 2.6. In that update, I removed > 140 lines of code. So, it's now even less of a change to OL. Maybe when 2.7 comes out, I'll figure out how to provide a patch that allows an extra argument to the OpenLayers.Map constructor that limits panning to the horizontal direction -- in which case genome-browser will cease to exist and only the single file containing OpenLayers.Layer.Genomic would be needed.

Maybe I'd need to make a real example of using openlayers for genomics, making more obvious use of layers, the vector stuff, fractional zoom, markers, geocoding, projections, power steering, etc to make it more obvious how badass OL can be. As far as I know, I'm the only one using it for genomics, and 98% of the time I spent developing it was in my spare time--I only have a good excuse to hack on it at $work when something breaks. If more people were using it, I might be more motivated to figure out how to do stacking of images vertically, but still restricting scrolling to the horizontal--to essentially allow the same thing as "tracks" in other genome browsers. Maybe that's the killer feature.

And

An observation on the difference between the bio and geo programming environments as I interpret them:
There's some tools that (IMHO) are better in the geo world than in the bio. Perhaps that's because of GDAL, the keystone for geo-data formats and projections. Since any other software (in any SWIG-able language) that makes use of gdal can access pretty much any format and I/O to any projection, then geo developers can do things like make nice renderers, or web-based map browsers rather than figuring how to convert that mrSID image in epsg:4326 to a tif in epsg:900913. (I'm also a GDAL fan-boy.)

Contrast this to the bio world where there's a bio for nearly every common language, bio-java, bio-perl, bio-python, bio-ruby, each with it's own blast parser, gen-bank parser, sequence objects, alignment objects. There is no keystone, so there's more duplication of effort, and there's no parser that's used across languagues as is the case for GIS. -- I'm not suggesting that shouldn't be the case, simply making an observation.

seqfind: levenshtein + bktree

2008-05-05T17:15:00.000-07:00

I've copied this recipe that I modified before, and added the BK Tree structure in cython. It's in my repo here.
Check it out with:

svn checkout http://bpbio.googlecode.com/svn/trunk/seqfind

or easy_install with

sudo easy_install http://bpbio.googlecode.com/svn/trunk/seqfind

It's now using the Damerau-Levenshtein distance which is more sensible for bioinformatics where transpositions are frequent.

Bearophile's original implementation used a tuple, which made sense, but in Cython, it's more efficient to use an object where the properties can be typed--as a class is converted to a c-struct--so there is no conversion when appending to a python array -- if i understand the generated c code correctly.

Using an object also allows arbitrary info to be passed along with the word when creating the tree, again, this is important for bio-informatics when the string is something like "actgcc ... acgtc" and it's useful to attach some annotation to it like:


words = [Word("actgc ... acgtc", {'name': 'At2g26540'}), ...]
#and then create a tree in the same fashion as with raw strings:
tree = BKTree(words)
# and search returns a word object:
[(w.word, w.info['name']) for w in tree.find("atc", 2)]

It'll suffice to say (the obvious) it's a lot faster searching with the tree, than comparing every word, on every search.

flash (n)back

2008-04-30T23:26:00.000-07:00

The PNAS article linked here found subjects could improve IQ with training. I've written a simple flash version of their protocol over a couple long evenings in haxe/flash.

The article methods list 2 stimuli, a moving box, and spoken letters. The test subject is to respond (click in my case) when the box position or the spoken letter is the same as it was 2 time steps ago. Where 2 is increased as the subject gets better. I didn't do sound, I just show a big letter. Clearly, the logical thing to do is use it for a couple weeks and then implement the sound when I'm smarter. (I've never really used swfmill, but I think that'd be useful here...)

The article is ambiguous about when the letter is to sound, I've made both the letter and the box appear at the same time. The default, as in the article is to have 3 seconds between events, and to show the box for 0.5 seconds. I also add some indication of whether the answer was correct (green +) or not (red -). That actually makes things more difficult as its distracting to see something new flash on the screen. It also keeps a running total of misses (didn't click when should have), correct, and incorrect (clicked when shouldn't have). The grid size is set at 3 * 3 as that's more than difficult enough for me, and it appears to be what they used in the article. It actually only has 8 positions as I use the center to display the letter. Another ambiguity from the article is the number of letters. I use 3. That's easily changeable in the code.

The length of the time step (time_step, 3000ms default), the amount of time to show the box and text (show_time, 500ms default) and the number of steps back (nback, 2 default) are settable via the url so the default equates to:
?time_step=3000&nback=2&show_time=500

Observations

1. It's freakin hard. I'm suck at it.
2. Haxe is nice, but jeez, I write ugly actionscript.
3. I like quick, pointless evening projects like this where I have a clear idea of the outcome. It's good for learning.
4. Whatever points there are against flash, it's easy to uh, "deploy".

Here's a live version of the app (untested in IE, but swfobject should do it's thing). Just click in the flash movie when there's something that's the same as 2 time-steps ago.
You can make it arbitrarily hard by sending in parameters on the url, see the links in the page for examples.

And here's the code. Get it from svn via:
svn co http://bpgeo.googlecode.com/svn/trunk/nback

levenshtein in cython

2008-04-27T22:43:00.000-07:00

EDIT(2): fix markup for <char *> casts ... fix malloc (see comments. thanks Bao).
NOTE: using a kwarg for limit slows things down. setting that to a required arg and using calloc for m2 speed things up to nearly as fast as the pylevenshtein.

Well, it seems to be popular to code up the levenshtein. I actually have a use for this and wanted to practice some Cython, so I've written a version. I used bearophile's recipe, wikibooks, this (from k4st on reddit) and this for reference. It follows bearophile's code closely, using only O(m) space instead of O(mn).


cdef extern from "stdlib.h":
    ctypedef unsigned int size_t
    size_t strlen(char *s)
    void *malloc(size_t size)
    void free(void *ptr)
    int strcmp(char *a, char *b)


cdef inline size_t imin(int a, int b, int c):
    if a < b:
        if c < a:
            return c
        return a
    if c < b:
        return c
    return b


cpdef int levenshtein(char *a, char *b, int limit=100):
    cdef int m = strlen(a), n = strlen(b)
    cdef char *ctmp
    cdef int i = 0, j = 0, retval
    cdef int achr, bchr
     
    if strcmp(a, b) == 0:
        return 0

    if m > n:
        ctmp = a;
        a = b;
        b = ctmp;
        #a, b = b, a
        m, n = n, m
        
    # short circuit.
    if n - m >= limit:
        return n - m

    cdef char *m1 = <char *>malloc((n + 2) * sizeof(char))
    cdef char *m2 = <char *>malloc((n + 2) * sizeof(char))
    
    for i from 0 <= i <= n:
        m1[i] = i
        m2[i] = 0

    for i from 0 <= i <= m:
        m2[0] = i + 1
        achr = a[i]
        for j from 0 <= j <= n:
            bchr = b[j]
            if achr == bchr:
                m2[j + 1] = m1[j]
            else:
                m2[j + 1] = 1 + imin(m2[j], m1[j], m1[j + 1])

        m1, m2 = m2, m1

    retval = m1[n + 1]
    free(m2)
    free(m1)
    return retval

I believe that's correct (if not, let me know!), it matches output from other versions (which i used fort testing) and it doesn't leak memory, so I must have done the malloc/free correctly despite my lack of C-fu.
Again, I'm not sure, but I think that's a pretty good example of how to mix python and C with cython as it's not much longer than the other python versions and pretty readable. And, it's very fast, it does 500K iterations of:
levenshtein('i ehm a gude spehlar', 'i am a good speller')
in 2.5 seconds.
bearophile's pysco-ed editDistanceFast took 3 minutes, 31 seconds to do the same (which is why I started this project...).

Literally, just before I was to post this, I found pylevenshtein which is even faster, does the 500K in 1.5 seconds. ho hum. It's doing more, handling unicode, and apparently doing some optimizations.... So, use that instead! -- and contribute a test-suite.
Next, I think I'll try a variant that allows transpositions (Damerau) and/or implement the BK-Tree, again cribbing off bearophile's recipe.

for reference, here are the contents of the setup.py to build the module.


from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

setup( name = 'levenshtein',
  ext_modules=[ Extension("levenshtein", sources=["levenshtein.pyx"] ), ],
  cmdclass = {'build_ext': build_ext})

numpy to tiff via gdal

2008-04-23T18:57:00.000-07:00

EDIT: 7 months later I came back to this and found an error. update in the code below with old line commented out.

Rather than venting about a project I've recently decoupled myself from, I'll try to do something constructive... I also posted this to the gispython mailing list, but I've had to figure it out a couple times, so I'll put it here for the record. Given an N * N numpy array, and a bounding box, it's actually fairly simple to make a georeferenced tiff:


from osgeo import gdal, gdal_array
import numpy
from osgeo.gdalconst import GDT_Float64

xsize, ysize = 10, 10
a = numpy.random.random((xsize, ysize)).astype(numpy.float64)

xmin, xmax = -121., -119.
ymin, ymax = 41., 43.


driver = gdal.GetDriverByName('GTiff')
# bad: out = driver.Create('a.tiff', a.shape[0], a.shape[1], 1, GDT_Float64)
# the args to Create are 'name', xsize, ysize. and .shape[0] is rows, which is y.
driver.Create('a.tiff', a.shape[1], a.shape[0], 1, GDT_Float64)

out.SetGeoTransform([xmin
                  , (xmax - xmin)/a.shape[0]
                  , 0
                  , ymin
                  , 0
                  , (ymax - ymin)/a.shape[1]])

gdal_array.BandWriteArray(out.GetRasterBand(1), a)

where the SetGeoTransform bit is (I believe) the same stuff you'd stick in a world file.

and plottable in pylab:


import pylab
tif = gdal.Open('a.tiff')
a = tif.ReadAsArray()
pylab.imshow(a)
pylab.show()

Anyone checked out matplotlib's basemap recently? I just wish they didn't rely on geos < 3...

python script as wsgi, cgi, or standalone

2008-04-11T19:00:00.000-07:00

EDIT:
See below for original, I realized this could be done cleanly with a decorator.
The decorator wrapplication takes the number of the port to use when called as a standalone server. The EMiddle class is unnecessary, it's just used as middleware to update the environ to show it came via wsgi. If there's a cleaner way, let me know.


#!/usr/bin/python
import os

class EMiddle(object):
    def __init__(self, app):
        self.app = app
    def __call__(self, env, start_response):
        env['hello'] = 'wsgi'
        return self.app(env, start_response)

def wrapplication(port):
    def wrapper(wsgi_app):
        if 'TERM' in os.environ:
            print "serving on port: %i" % port
            os.environ['hello'] = 'standalone'
            from wsgiref.simple_server import make_server
            make_server('', port, wsgi_app).serve_forever()

        elif 'CGI' in os.environ.get('GATEWAY_INTERFACE',''):
            os.environ['hello'] = 'cgi'
            import wsgiref.handlers
            wsgiref.handlers.CGIHandler().run(wsgi_app)
        else:
            return EMiddle(wsgi_app)
    return wrapper

@wrapplication(3000)
def application(environ, start_response):
    start_response("200 OK", [('Content-Type', 'text/plain')])
    yield "How do you like the teaches of peaches?\n"
    yield "from " + environ['hello']

ORIGINAL VERSION:

If you write a script with the application entry point that fits the wsgi spec, it's simple to make it run via mod_wsgi, cgi, or via standalone server depending on the context. I believe this is common knowledge, but for my own reference here's an example with the extra setup (which will work for any script) to do this:


#!/usr/bin/python

def application(environ, start_response):
    start_response("200 OK", [('Content-Type', 'text/plain')])
    yield environ['QUERY_STRING']

if __name__ == "__main__":
    try:
        from wsgiref.simple_server import make_server
        import sys
        port = int(sys.argv[1])
        print "server on port: %i" % port
        make_server('', port, application).serve_forever()
    except Exception, e:

        import wsgiref.handlers
        wsgiref.handlers.CGIHandler().run(application)

To change between wsgi and cgi toggle between

AddHandler cgi-script .py
AddHandler wsgi-script .py

For the stand alone server, just run it with an argument indicating the port:

python app.py 3000

to use the cherrypy server instead of the wsgiref, replace the make_server() line with:


from cherrypy import wsgiserver
server = wsgiserver.CherryPyWSGIServer(('0.0.0.0', port), [('/', application)], server_name='')
try:
    server.start()
except KeyboardInterrupt:
    server.stop()

That also handles the server shutdown more politely.

Genedex: query genomic features and sequence

2008-04-07T19:14:00.000-07:00

Normally, I don't write libraries, I figure smarter people than I should do such things, and I should just use them. But, I got tired enough of writing one-off scripts for genomic feature manipulation-- find the upstream, downstream neighbors and get the sequence -- and I saw enough of the pieces coming together that I decided to build it. I'd complained before about how rtree didn't support 1D indicies. Not only is this not a problem, it's beneficial. Genomic features should have strand information, so that's the 2nd dimension. Then rtree does containment queries, so it's simple to find only the features on a given strand. I realized this about the same time that the docstring for numpy's memmap went from 0 lines to about 100, and it was enhanced to take a filehandle, not just a filename. This means you can send in a start position and a shape to the numpy.memmap constuctor and it can create a numpy array of only that chunk. This means that it's possible to slice an unaltered fasta file using the numpy array syntax. That's very good.

So, if you put those 2 simple things together, you have the start of something powerful. That's what I did. Then I gave it a crappy name: Genedex (Gendex was taken) and slapped it into googlecode. Check it out: http://code.google.com/p/genedex/. My only design goal was to keep it as simple as possible. If the amount of features is under-whelming, that's good.

TDD

Also, I generally do TDD very half-ass, with asserts and maybe a couple doctests. However, I recently made fairly substantial changes to the SQLite datasource in featureserver, and wrote this set of doctests while doing so. It works! and I've been using it. So, I did what featuresever (presumably crschmidt) devs did and copied the setup for the shapely doctests. It's pretty useful for design, i'd just write out the code for how I wanted the API to look and then implement. The only thing is, for doctests, the way they're used (at least by me) is to copy the output from executing the code into the doctest. So, if your code is wrong to start with, you just copy the wrong answer into the doctest and it's broken but the tests pass. But, at least it's good for regressions, and I just had to remember not to blindly trust the output. That's true for all testing, but especially so for doctests.

So, there's now more tests than code. But, since it's mostly just tie-ing together pieces that do the real work, it's not much code. Doc-tests are also nice because (as the name suggests) they double as documentation. So, here's the genedex documentation:
http://genedex.googlecode.com/svn/trunk/doc/readme.html
It's pretty! It gets colored by pygments, using this script. The only major thing I'd like to add to the library is a plotting class using matplotlib. Then other smaller tasks like a method that takes 2 features and returns the sequence between them.
Any fixes, enhancements, ridicule, etc. will be greeted with commit access.

comparative genomics with openlayers

2008-04-06T21:37:00.000-07:00

Traditional genome browsers, look like this. In fact, I think that's the most popular genome-browser used--gbrowse. They display information in tracks, so any layer of annotation you just add on to the bottom of the image (after making the image taller). This doesnt work for genome-browser, the hack of openlayers to support only horizontal scrolling, because you if you have 2 adjacent tiles, if one has more features than the next, there's not guarantee that they'll be the same height, and no guarantee that a feature that's on both images will align correctly.

I was just hacking around, trying to test some work I'd done and realized that you can have annotation layers with OpenLayers, just add another map, and tie them together!

So that's 2 OpenLayers.Map() instances. What makes this easy is the new Map.panTo() methods in OpenLayers 2.6 (which is in release candidate 1). So, the top map registers for 'move' and 'zoomend' events with callbacks that update the bottom map with the position/zoom of the top map.
That's it! And layers of annotation are available, along with the slippy map. OpenLayers continues to amaze.
That site with the linked maps is here.

featureserver authentication

2008-03-30T21:39:00.000-07:00

As the name implies, featureserver serves vector features to various formats from a number of datasources, including OGR -- which means pretty much any vector format. That's extremely powerful. Really. That means, for instance, that when you're working on a really cool project and all anyone wants to know is if they can see it in KML/Google Earth, it's no extra work. Just point them to the REST-ful url like "http://example.com/featureserver/all.kml", and continue working on the cool project. Likewise for all.gml, .atom, etc. And, if you have a project with spatial data, if you put it in a format that featureserver understands, it's displayable, and editable in openlayers.

The next thing people want in a web application is some sort of user restrictions. In featureserver, by default, anyone can do any of the CRUD operations on any feature. I've been playing with a soon-to-be-open-sourced PPGIS (apparently the trendy acronym for that is now VGI) project where people can go to report sudden oak death. I want anyone to be able to report and view cases and add notes about existing cases, but only admins to be able to edit and delete existing reported sudden oak death cases. The simplest way is to use basic authentication in apache, but then anyone who goes to the site has to enter a user/password, and I think that really limits the public participation bit. If there's a way to do only authenticate sometimes with apache authentication, let me know.

unrested development

Since it's possible to use featureserver as an API, you can make your own server and add authentication in python. You can do this with any framework that supports sessions, I've done it with the development version (0.3) of web.py, using the nice sessions support. That it supports intuitive syntax for GET, PUT, DELETE, POST makes it a good fit as well. In the code below, the authentication related urls are /login and /logout, but those are never needed unless updating or deleting an existing point. Anyone can create new features. All featureserver related requests are made as with the original. Here's the wsgi script:


#!/usr/bin/python
import web
from FeatureServer.Server import Server
from FeatureServer.DataSource.SQLite import SQLite

urls = ( '/logout',      'logout'
        ,'/login',       'login'
        ,'/(.*)',        'features')

app = web.application(urls, globals())
session = web.session.Session(app
                , web.session.DiskStore('/tmp/sessions')
                , initializer={'authorized': False})

datasource    = SQLite('fsauth', file="/tmp/fsauth.sqlite")
featureserver = Server({'fsauth': datasource })

application = app.wsgifunc()

class login(object):
    """for a real app, save usernames, hashed pws in the db"""
    def POST(self):
        pw = web.input(password=None).password
        user = web.input(user=None).user
        if (user == 'abc' and pw == '123'):
            session.authorized = True
            return '[authorized]'
        return '[NOT-authorized]'

class logout(object):
    def GET(self): session.kill()

class features(object):
    """all the featureserver routing"""
    path = "/" + datasource.name + "/" # fsauth
    format = "geojson"
    def GET(self, feature_id=''):
        if "." in feature_id:
            feature_id, self.format = feature_id.split(".")

        # get web.py parsed url
        path = self.path + feature_id
        data = dict(web.input().items())
        data['format'] = self.format

        format, rsp = featureserver.dispatchRequest(data, path, "", request_method="GET")
        web.header('Content-type', format)
        return rsp

    def PUT(self, feature_id=None):
        return self.POST(feature_id, "PUT")

    def DELETE(self, feature_id=None):
        if "." in feature_id:
            feature_id, self.format = feature_id.split(".")
        # cant delete unless authorized.
        if not session.authorized: 
            web.header('Content-type', "text/plain")
            return "not logged in"
        path = self.path + feature_id
        data = dict(web.input().items())
        data['format'] = self.format
        format, rsp = featureserver.dispatchRequest(data, path, "", request_method="DELETE")
        web.header('Content-type', format)
        return rsp


    def POST(self, feature_id=None, method="POST"):
        if feature_id is None: return []
        if "." in feature_id:
            feature_id, self.format = feature_id.split(".")
        # must be an admin to do something with an existing feature.
        if not session.authorized:
            if not feature_id in ('new', 'create'):
                return 'not logged in'
        e = web.ctx.environ
        post_data = e['wsgi.input'].read(int(e['CONTENT_LENGTH']))
        path = self.path + feature_id
        format, rsp = featureserver.dispatchRequest({'format':self.format}, path, "", post_data=post_data, request_method=method)
        web.header('Content-type', format)
        return rsp

That'll be called from OpenLayers, but to demo, since everyone else is using curl which supports cookies:


FS="http://localhost/fsauth/"

echo "\n\nfirst borrow a feature. "
curl --url "http://featureserver.org/featureserver.cgi/scribble/35.geojson" > 35.geojson

echo "\nsee the features ... "
curl --url $FS

echo "\n\nadd a feature (no auth required.)"
curl -d @35.geojson --url "$FS/create.geojson"

echo "\n\nsee the features ... "
curl --url $FS

echo "\n\ntry to delete ... but cant"
curl -s -X DELETE $FS/1

echo "\n\nlogin ... \n"
curl -s --cookie-jar "cookies.txt" -d "password=123&user=abc" --url $FS/login > /dev/null

echo "\n\nthen delete ... "
curl -s -X DELETE -b "cookies.txt" $FS/1 > /dev/null

echo "\n\nsee the empty features ... \n"
curl --url $FS

echo "\n\nreturn the borrowed feature. thanks.  :-) "
curl -s -X PUT -d @35.geojson --url "http://featureserver.org/featureserver.cgi/scribble/35.geojson"

the important point there being that the DELETE fails until the user is logged in. I'm pretty sure adding authentication makes it un-REST-ful. but ???. Anyway, I won't lose any sleep over it.

reading code

2008-03-25T22:25:00.000-07:00

A fasta file of rice genomic sequence is 355MB. It's not easy to understand how large that is. This is an attempt to come up with a quick metric.
So, I downloaded Ulysses.
wc shows it to have 267235 words. Some googling says the average person can read 250 words per - minute. So that's 267,235 / 250 / 60 = 17.8 hours. Well, it's hard to believe anyone can really read Ulysses in 18 hours but... good enough.
So on the rice fasta file i ran:

grep -v ">" rice.fasta | wc -c

to get rid of the 12 header lines (1 per chromosome) and only count sequence (should be within 12 characters counting the extra new-lines). That gives 372,077,765 characters. The average word-size in ulysses is 5. I rounded up to 6. So, the rice sequence has the equivalent of 372,077,765 / 6 = 62,012,960 words
So, at 250 words per minute, it'd take:
62012960 / 250 / 60 = 4,134 hours to read the rice genome. That's 172 days. Also, from what I know, the plot is hard to follow.
Genome size varies widely among plants. I have a couple ideas for pointless visualizations of this...

point partitioning

2008-03-24T19:12:00.000-07:00

After you spend all day banging your head against your own problems, sometimes it's just nice to bang it on something else for a bit.
This question came through on the postgis mailing list and it seemed like a good diversion. I think it's a very clear description of the problem. To quote:

I have around 300,000 points, each with a lat, lon and altitude (also converted to geometry). I need to get a subset of those points, where none of them are within 5m (or some other arbitrary distance) of each other. It doesnt matter which points get picked over another, as long as whatever the data set it creates, that none of the points are within that 5m radius and that relatively most of the points are used

So, I hacked up a quick solution. It's probably inefficient -- deleting keys from a dict, and removing entries from an rtree index. But, it's easy to understand, and (without the plotting) it runs in about 2 minutes for the requested 300000 points.
When plotting, the image looks like this:

and the code...


import rtree
import random
import pylab


dist = 500
index = rtree.Rtree()
points = {}  # xy coords. use dict to delete without changing keys.
groups = {}  # store the grouped points. all keys are > dist apart

# create some random points and put them in an index.
for i in range(3000):
    x = random.random() * 10000
    y = random.random() * 10000
    pt = (x, y)
    points[i] =  pt
    index.add(i, pt)

print "index created..."

while len(points.values()):
    pt = random.choice(points.values())
    print pt
    bbox = (pt[0] - dist, pt[1] - dist, pt[0] + dist, pt[1] + dist)

    idxs = index.intersection(bbox)
    # add actual distance here, to get those within dist.

    groups[pt] = []
    for idx in sorted(idxs, reverse=True):
        delpt = points[idx]
        groups[pt].append(delpt)
        index.delete(idx, delpt)
        del points[idx]

# groups contains keys where no key is within dist of any other pt
# the values for a given key are all points with dist of that point.

for pt, subpts in groups.iteritems():
    subpts = pylab.array(subpts)
    pylab.plot(subpts[:,0], subpts[:,1], 'k.')
    pylab.plot([pt[0]], [pt[1]], 'ro')

pylab.show()

OGR python projection

2008-03-18T18:22:00.000-07:00

OGR Projection

If you're using shapely and you need to do projections, you'll either have a lot of boilerplate or a function like this one. Actually, even in OGR, there's a lot of bioler plate involved in transforming....


from osgeo import ogr
from shapely.wkb import loads

def project(geom, to_epsg=900913, from_epsg=4326):
    """utility function to do quick projection with ogr,
    to and from shapely objects
    >>> from shapely.geometry import LineString
    >>> l = LineString([[-121, 43], [-122, 42]])
    >>> lp = project(l, from_epsg=4326, to_epsg=26910)
    >>> lp.wkt
    'LINESTRING (663019.0700828594854102 4762755.6415722491219640, 582818.0692490270594135 4650259.8474613213911653)'
    """

    to_srs = ogr.osr.SpatialReference()
    to_srs.ImportFromEPSG(to_epsg)

    from_srs = ogr.osr.SpatialReference()
    from_srs.ImportFromEPSG(from_epsg)

    ogr_geom = ogr.CreateGeometryFromWkb(geom.wkb)
    ogr_geom.AssignSpatialReference(from_srs)

    ogr_geom.TransformTo(to_srs)
    return loads(ogr_geom.ExportToWkb())

import doctest
doctest.testmod()

spatially explicit metapopulation models in scipy.

2008-03-16T11:07:00.000-07:00

Making Pretty Pictures

I started to learn to program about 5 years ago running population ecology models in mathematica. Yesterday, I found an old mma notebook with a model modified to include differential parasitoid dispersal to adjacent host cells depending on the host density in those cells. It's bascially a discrete-time Nicholson-Bailey model. But in a grid of cells, where each cell contains a population of hosts (H) and parasitoids (P) that give birth, die, eat, and get eaten according to the NB model. Each generation, following birth/reproduction/predation, the hosts and parasitoids disperse. The hosts disperse equally to the 8 surrounding cells in their neigbhorhood. The parasitoids can move irrespective of host densities when the aggregation parameter (eta) is 0. When aggregation is 1 (eta == 1), the parasitoids move to adjacent cells in exact proportion to host densities in each of the surrounding cells. muH and muP (i'm too lazy to figure out how to write the symbols for mu) determine the proportion of individuals in each population that disperse. I no longer have mathematica, but I downloaded all 96 megabytes of the reader to open the model:

I'm no longer a big fan of mathematica, but, it does make it easy to use funky characters and show equations nicely.
There, the ListCorrelate does the dispersal. It's a convolution to "smear" out populations according to the kernels, which are defined by the muH, muP and for the parasitoid, eta.
It took literally 10 minutes to translate that to numpy/scipy.
I'll paste that code at the end.
The cool thing is that you can see complex spatial patterns arising within the gridded metapopulation, even when the sum of the densities across cells is constant. So, running the model, and plotting the time series (top), and the phase portrait (bottom) are shown.

The time-series is the average density of Hosts, Parasitoids across the metaopulation. It takes a while to stabilize, but they oscillate down to an equilibrium. The phase portrait is just plotting the mean density of hosts vs parasitoids for each generation. The generations are discrete, so the lines are merely to view the trajectory. It's hard to see from the scale of the time series, but the parasitoid cycles lag those of the host by about a generation.
But the model is running in a spatially explicit manner, and averaging loses all of that spatial data. So, starting at 200 generations, the program then saves a snapshot of the grid, with higher host densities in white, and lower values in black. This is done every third generation. After the model run, there's a directory with a bunch of images. It'd be nice to play them as a movie... Imagemagick's convert is perfect for this:


convert -delay 40 -loop 0 images/*.png metapop.gif

Where that command creates a nice little (4MB) animated gif:

That shows the spatiotemporal dynamics of the host population.
That's a pretty simple way to incorporate time into an visualization when you're otherwise limited to 2 dimensions... Maybe I'm too easily amused, but I could watch that all day.
Anyway, here's the hastily translated code:


import numpy
import pylab
from scipy.signal import convolve2d

def doplot(hm, pm):
    pylab.plot(range(gens), hm, 'b')
    pylab.plot(range(gens), pm, 'g')
    pylab.ylim(0, pm.mean() + hm.mean())
    pylab.legend(('hosts', 'parasitoids'))
    pylab.subplot(212)
    pylab.plot(hm, pm, 'r')
    pylab.plot(hm, pm, 'k.')
    pylab.xlim(0, 1.3 * numpy.mean(hm) + 5)
    pylab.ylim(0, 1.3 * numpy.mean(pm) + 5)
    pylab.xlabel('hosts')
    pylab.ylabel('parasitoids')
    pylab.show()

def metapop(H0, P0, a=0.05, l=3., K=1000, muH=0.2, muP=0.8
              , Hrange=1, Prange=1, eta=3, size=32, gens=1000
              , mode="same", boundary="wrap"):
    E = numpy.e
    Hmeta = []
    Pmeta = []
    Pkern = numpy.ones((2 * Prange + 1, 2 * Prange + 1)
            , dtype=numpy.double)
    Hkern = numpy.ones((2 * Hrange + 1, 2 * Hrange + 1)
            , dtype=numpy.double)
    Hkern *= muH/((2 * Hrange + 1)**2 - 1.)
    Hkern[Hrange + 1, Hrange + 1] = 1.0 - muH

    for gen in range(1, gens + 1):
        Hmeta.append(H0.mean())
        Pmeta.append(P0.mean())
        
        # poisson search.
        f = E**(-a * P0)

        # predation, births
        H1 = l * H0 * f * numpy.exp(-numpy.log(l) * H0*f/K)
        P1 = H0 * 1 - f

        # simple movement between adjacent cells.
        H0 = 0.001 + convolve2d(H1, Hkern, mode=mode, boundary=boundary)

        # biased movement by parasitoid according to density of hosts in
        # adjacent cells.
        heta = H0**eta
        B = heta / convolve2d(heta, Pkern, mode=mode
               ,  boundary=boundary)

        P0 = (1. - muP) * P1 + muP * B * convolve2d(P1, Pkern, mode=mode
                , boundary=boundary)
        P0 *= P1.sum()/P0.sum()
        if gen > 200 and not gen % 3:
            pylab.figure(figsize=(4,4))
            pylab.axes([0, 0, 1, 1])
            pylab.imshow(H0, cmap=pylab.cm.gray)
            pylab.xticks([])
            pylab.yticks([])
            pylab.savefig('images/%03i.png' % gen)
            pylab.close()

    return numpy.array(Hmeta), numpy.array(Pmeta)

if __name__ == "__main__":
    gens = 300
    size = 64
    H0 = 1000 * numpy.abs(numpy.random.randn(size * size).reshape(size, size))
    P0 = numpy.zeros_like(H0)
    P0[0, 0] = 1.
    P0[32, 32] = 1.

    pylab.close()
    pylab.subplot(211)
    hm, pm = metapop(H0, P0, gens=gens, a = 0.02)
    doplot(hm, pm)

The variable names are ugly, but it's a start for something to play with. It's interesting to see the effects of parasitoid aggregation (eta). And the parasitoid initialization--here, it introduces a single parasitoid at cells 0,0 and 32, 32. Should anyone actually use this, I have a derivation of the equations, and better explanation of the parameters available.

rtree: know your nearest neigbhors

2008-03-11T08:12:00.000-07:00

My computer spends a lot of time looking for neighbors of a given location-- even more so for bio, than for geo. This is what I've learned about the options for doing smarter search so far.


SELECT * from (
 (SELECT * FROM ((SELECT * FROM feature WHERE start<= ? ORDER BY start DESC  LIMIT 1) UNION (SELECT * FROM feature where start>= ?
ORDER BY start LIMIT 1)) as u)
 UNION
 (SELECT * FROM ((SELECT * FROM feature where stop<= ? ORDER BY stop   DESC  LIMIT 1) UNION (SELECT * FROM feature where stop>= ?
ORDER BY stop LIMIT 1)) as v)
) as w
ORDER BY ABS((start + stop)/2 - ?) LIMIT 1

if you fill in ? with an integer location, that query will return the closest feature most of the time. It's verbose if not ugly, and that's only for 1 dimension. It can return the wrong feature in certain cases.... You have to write it like that in MySQL, because it doesnt support functional indexes, so as soon as you do something like:
ORDER BY ABS((start + stop)/2 - ?)
it's no longer an indexed search
It's a hard problem, even if you're using postgis. And even if you're a postGIS badass.
Other than postGIS, there postgres's builtin geometric types. Even so, for most things, using SQL makes me feel far from the data.
Generally, I pull my data into a numpy array and use slicing to get to a region of interest:


feats = all_feats[(all_feats['bpmax'] > bpmin) & (all_feats['bpmin'] < bpmax)]

that's usually pretty fast. It can be even faster with numexpr. I suppose it could be done using indexed searching in pytables pro. Still doing a windowed search like that doesn't guarantee you get the nearest neighbor and you still have to pare it down and do some arithmetic to figure out the nearest gene. And, it can be the bottleneck.
For bioinformatics libraries, there's fjoin which does indexing, but doesnt seem much of a library (plus the indentation is messed up in the file). There's also pygr, which, presumably would be perfect for this, if I could understand it. I keep looking at pygr occassionally, but it just doesn't sink in yet.

This mention of rtree and nearest neighbor reminded me of this (see Sean's comment). Originally, I thought between shapely, rtree, and pylab, you could have a pretty useful python-specific plot-table, spatial data structure. It's a wrapper over spatialindex which is pretty comprehensive, supporting Rtree, MVR tree, TP Rtree (pdf). It's pretty low-level, you have to add each item to the index, but that lets you do powerful things, like add an index to uhhh, anything. As with most geo libraries, you have to hack them a bit to hand1e 1D shapes like genomics data, which only have a start and a stop, or maybe just a midpoint. But, I now have the simplest possible proof of concept to do stuff like this:


gstore = Gendex('genes')
for start, stop in features:
   gene = (start, stop)
   gstore.append(gene)
gstore.save()

and then to do nearest neighbor query:


gstore = Gendex('genes')
for loc in locations:
   neighbors = gstore.nearby(loc, 2)
   do_something(neighbors)

ideally, it'd allow for genomics upstream, downstream searches, with a syntax like:
gstore.nearby(loc, 2, direction=-1)
with -1 for upstream, and 1 for downstream locations.
Gendex is a simple class that automatically creates a store of locations and indexes them. Rtree can already create a persistent index, this just creates a persistent pickle of the start, stop locations. It has a bit of syntax sugar for 1D, automatically setting the y-values to 0.1 or -0.1 and returns the coordinates via .nearby, instead of their indicies from rtree's .nearest() search. Eventually, instead of just using a tuple location, it'll use take a class with name, strand, chromosome, etc. attributes but this is good enough for tinkering:


class Gendex(rtree.Rtree):
    _store = []

    def __init__(self, *args, **kwargs):
        self.filename = args[0]
        super(Gendex, self).__init__(*args, **kwargs)

        if not 'overwrite' in kwargs and os.path.exists(args[0] + '.pkl'):
            self.load()

    def append(self, bds):
        self.__setitem__(len(self._store), bds)

    def __setitem__(self, i, bds):
        # if < 2 y values, search performs poorly.
        y = (i % 2) and 0.01 or -0.01
        if isinstance(bds, (list, tuple)):
            assert len(bds) == 2
            self.add(i, (bds[0], y, bds[1], y))
        else:
            self.add(i, (bds, y))
        self._store.append(bds)

    def __getitem__(self, i):
        if isinstance(i, int):
            return self._store[i]
        assert(isinstance(i, (list, tuple)))
        return [self._store[ii] for ii in i]

    def save(self):
        f = open(self.filename + '.pkl', 'wb')
        pickle.dump(self._store, f, -1)
        f.close()

    def load(self):
        self._store = pickle.load(open(self.filename + '.pkl', 'rb'))

    def nearby(self, bds, n):
        if isinstance(bds, (int, float)):
            bds = (bds, 0.1)
        #return self.nearest(bds, n)
        return self.__getitem__(self.nearest(bds, n))

According to hobu, the spatial index library supports 1-d, it just has to be added to the rtree python wrapper. I've started to look at this, but given my limited c/c++ skillz, that may require a lot of hand-holding.

A nice lightweight bio-informatics package for feature management could be the combination of:

Pytables or a memmaped numpy recarray

rtree

pylab

that's it.
Pytables could work especially nicely if this gets implemented.

What a Shapely genome you have!

2008-03-04T07:00:00.000-08:00

This might be a case of if you have a really cool hammer, everything looks like a nail, but it was fun mixing tools from different disciplines.

After finding synteny, there's a bunch of paired genes whose neighbors are also pairs. Paired (homologous) genes have similar sequence because they have some function and can't change without loss of function. Non-gene sequence between the paired genes is mostly randomized via mutation, deletion, etc. But, there is non-gene sequence that is conserved between the genes. These CNS's-- conserved non-coding sequences--are usually sites that bind stuff that regulates the expression of a gene.
That looks like this.

With one gene on the top, and its pair below, both yellow. Pink lines in the foreground connect putative CNSs (similar sequences) between these genes. That the lines cross is bad. CNSs occur right at the level of noise. So even though a similar sequence occurs near both genes, it could be by chance. It is possible to reduce this noise using local synteny. In the figure, that means both ends of the line should be about equi-distant from the end of either yellow gene. And lines that are going diagonally across the image should be removed. The goal is to find the parallel lines and remove those that cross.

Luckily, there's a "bioinformatics" library that lets me write code like this:

for aline in pink_lines:
 for bline in pink_lines:
     if aline == bline: continue
     if aline.crosses(bline):
         aline.has_crossed.update(b)
         bline.has_crossed.update(a)

The .has_crossed is a set() to keep track which and how many lines any given line crosses. Then it's simple to find lines that crossed a lot of others and remove them.

# sort with most crosses first
pink_lines = sorted(pink_lines, cmp=lambda a, b: cmp(len(b.has_crossed), len(a.has_crossed)))

for cline in pink_lines:
 if len(cline.has_crossed) > THRESHOLD:
 cline.remove = True
 for dline in cline.has_crossed:
     # remove cline from all lines it crossed
     dline.has_crossed.difference_update(cline)

pinklines = [pl for pl in pinklines if not pl.remove]

Then there will still be some crosses remaining, so the code then loops through pink_lines, find intersections and removes based on the score of the hit, and how many crosses it has.

That is the simplest case. It also has to remove lines (CNSs) that are in the intron of one gene, and not in the intron of another:


if acns.within(agene) and bcns.within(bgene):
 do_something_with_intronic_cns(acns, bcns)
elif acns.within(agene) != bcns.within(bgene):
 # bad! in only 1 intron
 remove(acns, bcns)

There's also the cases where the lines of the CNSs don't cross, but the actual CNS's touch. So on one homeolog there might be a CNS from basepairs 920 to 936 and another from 922 to 937. We also want to get rid of those. So that'd be:

if cnsa.overlaps(cnsb):
 remove_either(cnsa, cnsb)

Finally, the CNSs themselves have to be syntenic relative to the gene--that is, if a CNS cns_a is 9000 basepairs upstream from gene a, and its homeolog, cns_b is 12 basepairs upstream from gene b, that's likely not a syntenic pair. As with synteny dot plots, sometimes it's good to flip the subject homeolog to the y-axis (in the images above the subject is the bottom gene) and keep the query (the top gene in the images above) along the x-axis. Then any hit between the x and the y gene can be drawn in the x-y space. For the genes above, that looks like this:

Starting with the top image. The thick blue line in the center of the image is the gene itself. Any large green dots without a black dot over them are what the program called as CNSs. Red CNSs are on the wrong strand and are not considered. CNSs that overlap the grey--either along the x-axis or the y-axis occur in the introns of other genes and are not considered. Green dots covering the blue gene are intronic CNSs. The thin blue line just extends the y=x line. The purple bowtie is what I use to enforce synteny (bowtie.contains(cns) ). Anything outside of the bowtie is either close to the query homeolog and far from the subject, or vice-versa. Real CNSs should fall along the diagonal. And we chose an arbitrary maximum distance from the gene of 12,000 basepairs. It's easy to see the syntenic genes line up along the y=x blue line.

The bottom image filled with lines of the bilious hue is to visualize the crosses (as in the code above). The line with the black stripes has been removed because it overlapped > THRESHOLD lines. After that removal, no other lines crossed and after removing non-syntenic genes outside the bowtie, only 20 CNSs remained.

After all that noise is removed, the original image looks like this in our real viewer:

and those are real CNSs with lovely parallel lines. This runs for thousands of homologous gene pairs and creates a list of CNSs. There are other complexities, i.e. when pairs that are on opposite strands or when there are other syntenic genes in the up/downstream regions.

Other than the extra dimension, MultiLineString() works perfectly for a single gene with many CNSs. The fact that it integrates so well with numpy, and therefore matplotlib is great for visualization. Compare that to making a map file for just a quick look some data. Anyway, a fun project. Cheers Sean and the makers of geos.

flash, vi, fcsh

2008-02-28T07:17:00.000-08:00

All my flash tinkering has been in VIM-- no IDE, no XML, just actionscript. It's a little tough to deal with the adobe compiler as it takes about 11 seconds to compile a large project like modestmaps on my machine. That's not good for a guess-and-check programmer. The typing does catch some errors.

The worldkit project compiles instantaneously with mtasc (the predecessor to haxe)--likewise for the as2 branch of modestmaps. The flash compiler shell drops the compile time for as3 modestmaps to under 3 seconds, so I've added this to my .bash_aliases:


alias fcsh="/usr/bin/rlwrap /opt/src/flex2/bin/fcsh"

the rlwrap is to use readline in the flash shell--meaning I can just press up-arrow to get the previous compile command. By default, one has to paste or type the entire command again.
With that, it's close to a reasonable workflow.

synteny mapping

2008-02-27T09:12:00.000-08:00

Living in Synteny
I've been working on automating synteny mapping between any pairs of genomes. Synteny is where there's a stretch of DNA or genes in some order on chromosomeA of organismX and due to a shared evolutionary history, you can find a similar stretch of genes in order on chromosomeB of organismY. Often there are small losses and inversions, but between closely related organisms like man and mouse, there's still a lot of synteny.
Plants can undergo polyploidy, following which, a species can have 2 entire copies of its genome. Over time, much of the duplicated cruft is lost, and the homologous chromosomes diverge, but if the divergence is not too great, it's still possible (actually common) to find synteny within the genome of a single organism--as well as between organisms.
I've written my own algorithm to find synteny which uses python sets, and numpy array slicing to do the heavy lifting. It is quite clever [wink]. And it _almost_ works but it's ... sorta non-deterministic.
The output looks like this for Arabidopsis against itself:

That's compressed from a 3000px*3000px image, but it's 25 boxes for the 5 vs. 5 chromsomes so above and below the diagonal should be mirror images. and with:
+ the diagonal dark line in the center being where each gene matches itself
+ the dark patches in the center of each box being all the crap in the centromere of each chromosome.
+ the speckled dots all over being because genes have many relatives. and there's lots of frequently appearing transposons.
+ the red lines being what my algorithm finds as diagonals.
The problem is it either extends too far into the sea of dots, or it doesn't seed on what should be a diagonal, depending on the parameters. The human eye is very good at this, but it's difficult to make a computer do it. Anyway, I'd previously tried published synteny finding programs and found them no better than mine, but just found dagchainer, where DAG is directed-acyclic graph. It's a but fussy in that given too many points, it will find spurious diagonals, but it's predictable, and fast, and it's a simple matter to cull the points before sending them in. DAGChainer is published and public domain, and the code is readable (made me think maybe C++ ain't so bad). Also, it doesnt rely on protein sequence as many synteny programs do. This is important when dealing with poorly annotated genomes. I had some trouble with it originally, and emailed the author, Brian Haas. He asked me to send my data set, so I sent the worst parts. He did some preprocessing, found some good parameters , ran it himself and sent me the results in under a day, including a couple of explanatory emails back and forth. So, now, I'm building my processing scripts around dagchainer, and it's working out great. Brian is extremely enthusiastic and helpful. Must be another one of those crazy folks who enjoy what they do.

open source gis and flash maps part two

2008-02-26T17:45:00.000-08:00

Mash Flap
I started looking into flash mapping stuff lately. For the patch I submitted to worldkit, Mikel gave me commit access! So, now the svn version of worldkit can be compiled with mtasc by typing "make". I feel unreasonably proud of that, given that much of what I did was some global search and replace stuff in VI, and then read and fixed mtasc compiler errors until they went away. It was good fun.
Michal Migurski saw in that same post that I mentioned modestmaps and gave me some good ideas on getting WMS going. I just figured out how to get that working and posted a message to their overly web2.0 forums. Hopefully someone with some real actionscript skillz will clean it up. A mapping library without a good WMS interface is much less useful for most of the stuff I do.

I haven't decided whether to use modestmaps or worldkit, or both. The time stuff the Mikel has done in worldkit is very cool and I haven't really looked at that yet. But I have a time based project starting soon. But, most of my GIS project are reliant on good, hi-resolution imagery and quality roads data, and it's nice if that's just "magically" included via one of the commercial map providers. That's a plus for modestmaps. Though I _still_ do not understand if the modestmaps usage of google, yahoo, microsoft images falls within the terms of use... Especially since the google starts sending "X" tiles instead of imagery after browsing with modesmaps for a while. Meanwhile, I keep using OpenLayers, because that's what sane people do.

fast python with shedskin

2008-02-24T09:36:00.000-08:00

There's a new release of the shedskin compiler. It is able to generate fast shared libraries that can be run from CPython. It can also create binaries so I thought I'd see how it did on some code from this BMC bioinformatcs article compared to psyco and CPython.
I took this iterative, brute-force (Needleman-Wunsch?) alignment code and modified it slightly. That's pasted here. (Notice the first line! that's how it appears in the original code). The modifications allow shedskin to infer the function and variable types. Plus, there's a couple changes I made that improve the run-time for all cases. The max() function is also in the original, but unnecessary because of python's builtin max(), however, pysco does run much faster using their hand-coded max(). For the shedskin run, I removed that extra code and used shedskin's builtin 'cause it made me feel better.

The python code was run as
$ time python -c "import alignment; alignment.imain()"

To run with psyco. The only change from the pasted script is to add
'import psyco; psyco.full()' at the top.

To run with shedskin built executable (after removing the max()):
$ shedskin alignment.py ### generates Makefile
$ make
$ time ./alignment

To run with shedskin as a shared (.so) library
$ shedskin -e -n alignment.py # generates Makefile
$ make
$ time python -c "import alignment; print alignment.imain()"
# python finds and imports the shared library (.so) before the .py file.

Timing:
time reported is real.
Python 2.5.1: 19.030s
Psyco: 1.1336s
Shedskin shared: 0.921
Shedskin binary: 0.818

CPython is much slower at 19seconds compared to ~1 for pysco and shedskin. The output format for the shedskin binary is slightly different because it calls the stuff in __main__, but they all generate the same alignments. It's interesting to look at the alignment.ss.py that shedskin generates, as you can see all the inferred types. The alignment.cpp contains the generated cpp code, which is also quite readable. It's also nice to be able to get a binary executable without jumping through any extra hoops.

Shedskin was very easy to use and faster than psyco for this case, I just pulled from SVN and it worked out of the box. It now has support for sets and regular expressions, and seems quite active. I can see using shedskin for purely brute force stuff where numpy is no help and I might otherwise have to resort to cython. I kinda like cython, but it's nice just to be able to get fast code with little to no modifications.

EDIT:
jython 2.3a0 runs this unaltered in 42 seconds.

Python Mapscript Tricks

2008-02-23T10:21:00.001-08:00

I mentioned previously a site which uses google maps with a WMS. The problem with using points (or labels) in a tiled application such as google maps, or OpenLayers is that each tile can only draw its own contents. So if you draw a point with a radius of 10 pixels whose center is 2 px away from the edge of the tile, then 8px of the entire 20px will be chopped. By default, those lost 8 pixels will not be drawn in the adjacent tile because the center of the point does not fall in that tile. so that gives something that looks like this for 2 adjacent tiles:

It's hard to even tell that those tiles belong together! Using some python mapscript and PIL, it's pretty simple to make those look like this( except I dont know how to tell blogger not to add the spacin...) :

That's actually the same WMS request(s), just changing the call from a simple WMS CGI script that calls mapserver to a WSGI script that uses python mapscript and PIL.
The script:
1. takes the current bounding box, width and height, and expands them all proportionally.
2. has mapscript draw the expanded image which is
3. saved into a stringIO object which is
4. sent that to PIL to be cropped down to size and extent.
5. returned to the browser.


#!/usr/bin/python

import mapscript
from cgi import parse_qsl
import os
import Image
from cStringIO import StringIO

BUFFER_PCT = 0.1
MAPFILE = 'soil.map'

########################################################
MAPFILE = os.path.join(os.path.dirname(__file__), MAPFILE)

def application(environ, start_response):
  winput = dict(parse_qsl(environ['QUERY_STRING']))
  format = winput.get('FORMAT', 'image/png')
  start_response('200 OK', [('Content-Type', format)])
  extension = format[format.find('/') + 1:]


  wms = mapscript.mapObj(MAPFILE)
  req = mapscript.OWSRequest()

  for k,v in winput.items():
      req.setParameter(k, v)

  buffer = BUFFER_PCT/2

  bbox = map(float, winput['BBOX'].split(","))
  rangex = bbox[2] - bbox[0]
  rangey = bbox[3] - bbox[1]
  w = int(winput['WIDTH'])
  h = int(winput['HEIGHT'])
  xdelta = int(round(w * buffer)) # e.g add 13px on each side
  ydelta = int(round(h * buffer)) # for 256x256 image

  bbox[0] -= rangex * buffer # extend the
  bbox[1] -= rangey * buffer # bbox in
  bbox[2] += rangex * buffer # all
  bbox[3] += rangey * buffer # directions

  # http://trac.osgeo.org/mapserver/ticket/2299
  req.setParameter('WIDTH',  str(w + 2 * xdelta))   # and adjust the
  req.setParameter('HEIGHT', str(h + 2 * ydelta))   # h, w by the same
  req.setParameter('BBOX', ",".join(map(str,bbox))) # amount
  req.setParameter("STYLES", "")
  req.setParameter("REQUEST", "GetMap")

  # PIL doesnt like interlace.
  wms.outputformat.setOption('FORMATOPTIONS', 'INTERLACE=OFF')
  wms.loadOWSParameters(req)

  im = Image.open(StringIO(wms.draw().getBytes()))
  if im is None: return ['']

  # crop the image back to the requested w, h
  im = im.crop((xdelta, ydelta, w + xdelta, h + ydelta))

  buffer = StringIO()
  im.save(buffer, extension)
  buffer.seek(0)
  data = buffer.read()
  return [ data ]

That's pretty simple, but pretty cool. And it actually runs slightly faster than the CGI version. It runs even more quickly if "wms = mapscript.mapObj(MAPFILE)" is moved into the global scope, so the mapfile doesnt have to be re-parsed each time, but I didnt test enough to make sure that's ok.

Open Source GIS

2008-02-21T16:59:00.000-08:00

Just saw this linked from Sean's post. The quotes in there are absurd. But reminded me there's an important point in favor of OSS that I haven't seen.

You can still get help directly via the email list or IRC from F Warmerdam, it's primary author. Likewise for the developers of Mapserver and OpenLayers and PostGIS. I wonder if the lead developers of ESRI products spend their off-time perusing the forum's or mailing lists and answering questions? ( I don't know, they may. But I suspect not. )
There's something to be said for enjoying what you do. And I think that's very true in the case of those in the open-source community. Happy coders make better software. Any programmer that denies that will leave me flabbergasted.

I just don't understand how I could be effective only clicking the menus that were provided to me if I chose a black-box solution. But, call me crazy, I like linux.

Also, in these parts, if you lock your single-speed bike to a wooden chair, it'll get stolen regardless of the cost of the lock or the weight of the bike.

Flash-y Map-plication vs. 500 marker limit

2008-02-12T17:21:00.000-08:00

I've been pushing capabilities of javascript mapping frameworks. There's a real limit on the number of markers a browser can display without getting bogged down. 500 is a good limit, and that's probably too high for internet explorer. You can put as many markers as you like via WMS tiles. You can even take map clicks, send them back as AJAX queries and open an info window. I recently helped my friend do this, and we have a pretty snappy map displaying 5,000+ clickable "markers" no problem. It's actually not markers, just tiles, but it works quite well:
That's a snapshot of the google maps application, with an info window that appears when a marker is clicked. It sends a GetFeatureInfo request back to a mapserver WMS.

Vector drawing in the browser is limited. The amazing efforts of the OpenLayers, featureserver projects make this difficult to assert, as they abstract away all of the browser incompatibilities and give a nice platform to do real vector editing in a map. Still, there's a limit, and even if browsers become twice as fast in the next 2-3 years (and that's a big if, because the slowest commonly used browser will be the weakest link), the limit will be at 1000 markers.
Flash is built to do vector drawing. It's cross-browser. I've got my toes wet with a flash-based project and now that I am passed (part of) the learning curve, I figured I'd see what's available for GIS applications in flash...

Yahoo fairly recently released their Actionscript 3.0 version of their maps API. The examples look very straight-forward, if a bit heavy on the under-scores--even for a pythonista. But it's not open source, and there is no indication on how one would request tiles from any source--say WMS.
There is a 3D earth viewer with a nice javascript API, but it's not open source and I dont need the 3D stuff.

I'm still looking at worldkit, and even provided a patch to allow it to compile with mtasc. The worldkit approach as I understand it is to just provide the functionality in the flash movie, and only offer customization through an inutitive config.xml file, so no actionscript programming is necessary. Of course, this isn't strictly the case, but I think it's the way that it's most used. It's licensed GPL...

Then there's modest maps who are truly open source, with a trac bug tracker and BSD licensed code, multiple developers, clean design, and libraries for both actionscript 2.0 and 3.0. And they can display imagery from any of the major providers, I don't quite understand how they get by without violating the licensing restrictions, but they are making javascript calls--presumably to get the copyright info associated with the region of interest. They have a coordinate conversion system that seems very clean, though I also dont understand how to relate it to proj4 definitions--or the EPSG code in a WMS request. If I can understand how to use this, and even if it's possible to have layers that overlay base imagery I'd go with modest maps.

Anyone know of any others? I guess my criteria are:

Active, currently developed
Open source
WMS friendly

EDIT:
I was able to get a modestmaps actionscript 2.0 movie to show 2500 markers, with the map still very responsive. The as3 version should perform even better because a marker can be a subclass of the more lightweight Sprite, rather than the full movieclip.
this is a shot of the randomly placed markers.

python bioinformatics

2008-02-10T09:58:00.000-08:00

There's a new article out in BMC Bioinformatics with a comparison of the speed and length of programs from various languages. This article was sent to the biology in python (BIP) mailing list. Looking at the code, it's not that bad, but it is clear that the authors are not pythonistas, and that the reviewers have done a great job on the actual paper, but likely there was no thorough review of the code.
The authors define their own max() function that needlessly overrides python's built-in, and they use the code:

line.rstrip('/n')

that indicates there was not a thorough understanding of python, or a complete code review. Even a non-python programmer should have seen the intent was to strip a newline '\n', but the operation is not inplace, so the desired behavior could be achieved by:

line = line.rstrip('\n')

Syntactical mistakes aside, python was given a poor review on speed. Andrew Dalke, of wide-finder (and general python-bio) fame ran the alignment program (which seems to be Needleman-Wunsch) and found the psyco JIT'ed version to run in 1.7 seconds instead of the original 18+.

Given the lack of polish on the programs, even from other languages, the question then is what's the recourse? On the BIP list, some even suggested requesting that the authors retract the article. I think that's going too far as (to my knowledge) the code runs, and does more/less what it should. An alternate approach may be to propose more thorough code review for all articles, not just those that are benchmarks. As an example of informal code review, had the authors sent their code to the python list, they would undoubtedly have received numerous suggestions for making the code faster and more idiomatic. Likewise for the other languages. Clearly, that's not a solution for all cases, but given reasonably active communities for Bio-Python, Perl, Java, Ruby, a journal or author should be able to find a competent reviewer--especially in the most common scenarios where a single language was used. Journals demand a very specific format of the text of an article, should have similar standards for any code that is used in the article? Should they require automated tests? This presents a problem for those that use proprietary software. But one of the points of a scientific paper is to document the methods sufficiently to reproduce the results, is the actual code required to do so?
It's an interesting question, I don't know the answer. I notice 2 things:
1) A reviewer is not expected to know how to program, but she is expected to understand the science.
2) We do not have the same standards for code review as for reviewing the text.

Do the ends justify the means?

parallel blasts using python's pp module

2008-02-02T16:15:00.000-08:00

BLAST handles utilizes multiple cores for some scenarios using the -a flag. however, i often do full genome blasts -- blasting all chromosomes of one organism against all others. For the case of rice (O.s.) against itself, this is 124 jobs. there are simple tools in python to run a queued blast in which on my 8 core machine, each core will run 1 of those blasts, as a job finishes, the pp or parallel python module starts the next job, based on the number of cpus it has detected for your machine.
the syntax for the script is:

python pblast.py rice_rice_10kmers

where rice_rice_10kmers is the section in a config/.ini file to get the parameters.
My fasta directory looks like:


$ ls /tmp/rice/fasta/
ricetenkmers_chr01.fasta  ricetenkmers_chr04.fasta  ricetenkmers_chr07.fasta  ricetenkmers_chr10.fasta  ricetenkmers.order
ricetenkmers_chr02.fasta  ricetenkmers_chr05.fasta  ricetenkmers_chr08.fasta  ricetenkmers_chr11.fasta
ricetenkmers_chr03.fasta  ricetenkmers_chr06.fasta  ricetenkmers_chr09.fasta  ricetenkmers_chr12.fasta

the script will create (144 - repeats ) output files with names that look like:

ricetenkmers_chr06_vs_ricetenkmers_chr09.blast

The code below is what does the work:


"""
this will parallelize the blast on a directory of fasta files. keeps
the number of CPUs on the given machine full until all jobs are done.
"""

import commands
import os
import sys
import pp
import glob

import ConfigParser


filen = lambda fullpath: fullpath[fullpath.rfind("/") + 1:fullpath.rfind(".")]
"""
>>> filen = lambda path: fullpath[path.rfind("/") + 1:path.rfind(".")]
>>> filen("/tmp/rice/fasta/ricetenkmers_chr06.fasta")
'/ricetenkmers_chr06'
"""

def gen_command(q_fastas, s_fastas, format_db, blast, out_dir):
   """
   generator of blast commands given arguments:
   'q_fastas' : a list/iterable of strings indicating the full path
                to the set of query fastas.
   's_fastas' : a list/iterable of strings ... subject ...
   'format_db': a formatdb command (see pblast.ini for example)
   'blast'    :  a blast command string (see pblast.ini) 
   'out_dir'  :  directory to write the blast output files
   yields the full blast commands.
   """

   for q_fasta in q_fastas:
       q_name = filen(q_fasta)
       for s_fasta in s_fastas:
           s_name = filen(s_fasta)

           out_name = os.path.join(out_dir, "%s_vs_%s.blast" % (q_name, s_name))
          
           if not os.path.exists("%s.nin" % s_fasta):
               print >>sys.stderr, "formatting %s" % s_fasta
               # formatting db is not parallelized here ...
               commands.getoutput(format_db % s_fasta)

           # tell it where to send the output.
           command = (blast % (q_fasta, s_fasta)) + " -o " + out_name
           yield command



def consume(command):
   return command, commands.getoutput(command)

def save_blast_info(out_dir, fomatdb, blast):
   """ keep a record of the blast and format_db params used"""
   bp = open(os.path.join(out_dir, "00blast.params"), "w")
   print >> bp, format_db
   print >> bp, blast
   bp.close()



def main(q_fastas, s_fastas, format_db, blast, out_dir):
   # create the server.
   s = pp.Server()
   if not os.path.exists(out_dir): os.makedirs(out_dir)
   save_blast_info(out_dir, format_db, blast)

   jobs = [] # this will hold all the pp jobs.

   for c in gen_command(q_fastas, s_fastas, format_db, blast, out_dir):
       # tell pp that we'll call the consume function
       # (c,) is the argument list, where c is a string for hte blast
       # command
       # ("commands",) is a tuple of modules that will be needed for the job
       jobs.append(s.submit(consume, (c,), (), ("commands",)))

   # loop through and collect the jobs and any output as they return
   for j in jobs:
       command, output = j()
       print command
       if output: print output


if __name__ == "__main__":
   config = ConfigParser.ConfigParser()
   config.read("pblast.ini")
   params = dict(config.items(sys.argv[1]))

   format_db = params.get("format_db", "/usr/bin/format_db -p F -i %s")
   blast     = params.get("blast", "/usr/bin/blastall -p blastn -K 80 -i %s -d %s -e 0.001 -m 8 ")
  
   query     = params.get("query_files")
   subject   = params.get("subject_files")
   out_dir   = os.path.join(params.get("out_dir", "/tmp") , sys.argv[1].strip())

   fasta_list = {"q": [], "s": []}
   patterns = ("fa", "fasta", "faa", "fas")

   for qs, fileset in (("q", query), ("s", subject)):
       # it's either a directory, in which case get all fasta files.
       if os.path.isdir(fileset):
           for pat in patterns:
               fasta_list[qs].extend(glob.glob(os.path.join(fileset, "*." + pat)))
       # or it's a glob pattern
       else:
           fasta_list[qs] = glob.glob(fileset)
       assert len(fasta_list[qs]) > 0, "didn't find any files for %s" % fileset
          

   main(fasta_list["q"], fasta_list["s"], format_db, blast, out_dir)

which reads this config file:


[rice_rice_10kmers]

# is it protein? leave the %s to be filled programmatically
format_db=/usr/bin/formatdb -p F -i %s

# just use correct path and use all parameters here:
# the %s's get filled with query and subject files programmatically.
# the output file is chosen and appended programmatically.
blast=/usr/bin/blastall -p blastn -K 80 -i %s -d %s -e 0.001 -m 8

# where to send the blast output becomes:
# /tmp/rice_rice_10kmers/
out_dir=/tmp/

# a glob pattern or a directory
query_files=/tmp/rice/fasta/
# it's a self-self blast. so query and subject are same.
subject_files=/tmp/rice/fasta/

This will run a blast of all rice chromosomes (split into 10,000mers) against all other 10,000 mers. This is what top looks like on my 8-core machine.

tinycc

2007-11-26T15:21:00.000-08:00

i've been _trying_ to learn C. tinycc, beside being tiny, it compiles very quickly, allowing you to do cool things like script in C

#!/usr/bin/tcc -run

#include 

int main(int argc, char *argv[]) {
 printf("Hello World %s, %s", argv[0], argv[1]);
 return 0;
}

and then run as
./file.c arg_1

which makes it easier for those c-fu to guess and check.
it also allows such nice things as c in python
which is like pyinline, but uses ctypes and doesn't need write access.

Sorting by proximity to a date in PostgreSQL

2007-10-11T20:53:00.001-07:00

postgreSQL has great support for dates,

=> SELECT '2007-08-23'::date - '2006-09-14'::date as days;
days
------
343

given a date column and a date, to find the nearest date, you can "extract the epoch", here, i used ABS as i just want the nearest date, before or after.

SELECT *, ABS(EXTRACT(EPOCH FROM(date - '2006-08-23'))::BIGINT) as date_order FROM record WHERE well_id = 1234 ORDER BY date_order limit 1

i suppose this could make a nice PL/PGSQL function...

k-means clustering in scipy

2007-09-29T10:22:00.000-07:00

it's fairly simple to do clustering of points with similar z-values in scipy:


import numpy
import matplotlib
matplotlib.use('Agg')
from scipy.cluster.vq import *
import pylab
pylab.close()

# generate some random xy points and
# give them some striation so there will be "real" groups.
xy = numpy.random.rand(30,2)
xy[3:8,1] -= .9
xy[22:28,1] += .9

# make some z vlues
z = numpy.sin(xy[:,1]-0.2*xy[:,1])

# whiten them
z = whiten(z)

# let scipy do its magic (k==3 groups)
res, idx = kmeans2(numpy.array(zip(xy[:,0],xy[:,1],z)),3)

# convert groups to rbg 3-tuples.
colors = ([([0,0,0],[1,0,0],[0,0,1])[i] for i in idx])

# show sizes and colors. each color belongs in diff cluster.
pylab.scatter(xy[:,0],xy[:,1],s=20*z+9, c=colors)
pylab.savefig('/var/www/tmp/clust.png')

using python mapscript to create a shapefile and dbf

2007-07-17T15:50:00.000-07:00

i always have trouble remembering how to use mapscript. it's pretty simple, but the docs are hard to find and the test cases (though excellent!) have a lot of abstraction.

heres some code that creates a shapefile and dbf (using another module). and does a quick projection at the start.


import mapscript as M
import random
from dbfpy import dbf

#########################################
# do some projection
#########################################

p = 'POINT(466666 466000)'
shape = M.shapeObj.fromWKT(p)
projInObj  = M.projectionObj("init=epsg:32619")
projOutObj = M.projectionObj("init=epsg:4326")
shape.project(projInObj, projOutObj)
print shape.toWKT()


#########################################
# create a shapefile from scractch
#########################################
ms_dbf = dbf.Dbf("/tmp/t.dbf", new=True)
ms_dbf.addField(('some_field', "C", 10))

ms_shapefile = M.shapefileObj('/tmp/t.shp', M.MS_SHAPEFILE_POLYGON)

for i in xrange(10):
    ms_shape = M.shapeObj(M.MS_SHAPE_POLYGON)
    ms_line  = M.lineObj()


    for j in xrange(10):
        ms_line.add(M.pointObj(random.randint(0,99), -random.randint(0,99)))

    ms_shape.add(ms_line)
    ms_shapefile.add(ms_shape)


    rec = ms_dbf.newRecord()
    rec['some_field'] = 'hi' + str(i)
    rec.store()

ms_dbf.close()

Note to self: using python logging module

2007-06-27T15:46:00.000-07:00


import logging

logging.basicConfig(level=logging.DEBUG
        ,format='%(asctime)s [[%(levelname)s]] %(message)s'
        ,datefmt='%d %b %y %H:%M'
        ,filename='/tmp/app.log'
        ,filemode='a')

logging.debug('A debug message')
logging.info('Some information')
logging.warning('A shot across the bows')

Fix indentation in VIM

2007-04-26T08:54:00.000-07:00

Often times, i get which has the indentation completely messed up, not just mixing tab/spaces, but really "whack"

these commands seem to magically fix for at least 2 test cases:


:set filetype=xml
:filetype indent on
:e
gg=G

Using Python MiddleWare

2007-04-25T15:59:00.000-07:00

just trying to figure this stuff out. it's pretty simple, but there's one level of abstraction through web.py. you can use middleware to add keys to the environ for example.
http://groovie.org/files/WSGI_Presentation.pdf


#!/usr/bin/python

import web
import random

class hi(object):
   def GET(self,who='world'):
       web.header('Content-type','text/html')
       print "hello %s" % who

class bye(object):
   def GET(self,who='world'):
       web.header('Content-type','text/plain')
       print "bye %s" % who

       for c in web.ctx.env:
           print c, web.ctx.env[c]

class other(object):
   def GET(self):
       web.header('Content-type','text/plain')
       for c in web.ctx:
           print c, web.ctx[c]

urls = ( '/bye/(.*)', 'bye'
       ,'/hi/(.*)' , 'hi'
       , '/.*'     , 'other')

class RandomWare(object):
   def __init__(self, app):
       self.your_app = app;

   def __call__(self,environ,start):
       environ['hello'] = random.random()
       return self.your_app(environ,start)

def random_mw(app):
   return RandomWare(app)

if __name__ == "__main__":

   web.run(urls,globals(),random_mw)

Install run, and benchmark mod_wsgi in < 10 minutes

2007-04-24T16:43:00.000-07:00

svn checkout http://modwsgi.googlecode.com/svn/trunk/ modwsgi
cd mod_wsgi
./configure
make
sudo make install
# note where mod_wsgi.so went on your system
echo "LoadModule wsgi_module /path/to/mod_wsgi.so" >> /path/to/apache2.conf

mkdir /var/www/wsgitest/
cd /var/www/wsgitest/

vi .htaccess
# [in .htaccess]
Options +ExecCGI
< Files hi.py >
SetHandler wsgi-script
</Files>
# [ end .htaccess]

vi hi.py
# [in hi.py]


#!/usr/bin/python

import web

class hi(object):
def GET(self,who='world'):
  web.header('Content-type','text/html')
  print "hello %s" % who

class bye(object):
def GET(self,who='world'):
  web.header('Content-type','text/html')
  print "bye %s" % who

urls = ( '/bye/?(.*)', 'bye'
  ,'/hi/?(.*)' , 'hi' )

application = web.wsgifunc(web.webpyfunc(urls, globals()))

#[end hi.py ]

you can then browse to
http://localhost/wsgitest/hi.py/hi/there
# see "hello there"
http://localhost/wsgitest/hi.py/bye/bye%20bye
# see "bye bye bye"

(meaningless) Benchmarking:
change last line in hi.py to:
if __name__ == "__main__": web.run(urls,globals())
and save as cgi.py

cgi
$ ab -n 1000 -c 30 http://localhost/wsgitest/cgi.py/hi/there | grep 'Requests per second'

Requests per second: 4.08 [#/sec] (mean)

wsgi
$ ab -n 1000 -c 30 http://localhost/wsgitest/hi.py/hi/there | grep 'Requests per second'

Requests per second: 351.05 [#/sec] (mean)

vim tricks

2007-03-21T09:45:00.000-07:00

i've been trying to learn new stuff in vim, instead of doing same old. recently, i've been using :tabe to edit in tabs. lately, i've been trying the :sp to edit in splits.
this set of tricks makes it even nicer:
http://www.vim.org/tips/tip.php?tip_id=173
now i can type ctrl+j to move down or ctrl+k to move up a split and
have that split maximized.
both tabs and split make it simple to yank and paste between files. something for which i had been using the mouse.

postgresql and mysql: benchmark? how?

2007-01-31T12:28:00.000-08:00

so somehow, my previous post on postgres / mysql made reddit , which i happened to be reading yesterday afternoon. i didnt even realize it was my post until following the link.
there were a couple harsh comments stating that i found what i wanted to find. ... which were merited given the sensationalist way i presented the results (50%) and the careless use of the term "benchmark". and yes, the config for mySQL was the default. still, i just presented what i found.
i was surprised noone commented on the hackish way that i checked to see if it was a protein sequence in perl, rather than mysql--or the coolness of pre-fetching in DBIx (which is available as eager loading or setting lazy=False in the mapper in python's sqlalchemy).

re the comments on things to change in the postgresql.conf... i'll try at some point. are there any suggestions for mysql?
the machine has 12G ram, 4CPUs. likely, the raid configuration (i dont know how it's set up) is not optimal, but that is out of my hands.

for a "real benchmark" it'd be nice to do this in sqlalchemy with the schema written in python and then just change the database engine between mysql/postgresql/sqlite. any pointers on how one would go about creating a "real benchmark"?

actually, that would be a good ask reddit topic: "How to design a 'real database benchmark'?"

real-world postgresql vs mysql benchmark

2007-01-14T21:48:00.000-08:00

At my work, we have a large MySQL database (15 MyISAM tables, 21 million rows, 10Gigs size). After seeing the benchmarks showing that Postgres out-performs MySQL on multi-core machines (our new db server has 4 CPU's), I ported the database to PostgreSQL.
We have begun using the DBIx perl module since Class::DBI is too sloooow. The DBIx module allows closer access to the generated SQL, and it allows "prefetch"ing which eliminates extra back-and-forth (and object creation) between the server and client. In addition, the connection string is in the script, not in the generated API. This makes it easy to benchmark as all that is required to change between db engines is to change the connection string. Using this script:


use CogeX; use strict;
# mysql
#my $connstr = 'dbi:mysql:genomes:host:3306';
# postgresql
my $connstr = 'dbi:Pg:dbname=genomes;host=host;port=5432';
my $s = CoGeX->connect($connstr, 'user', 'pass' );

my $rs = $s->resultset('Feature')->search({
          'feature_type.name' =>   'CDS' ,
          'feature_names.name' => {like => 'At1g%' }
      },
      {
          join => ['feature_names','feature_type'],
          prefetch => ['feature_names','feature_type']
      }
);

while (my $feat =$rs->next()){
  my $fn = $feat->feature_names;
  my $type = $feat->feature_type->name;

  map { print  $_->name . ":". $type . "\t" } $fn->next();
  print  "\n";

  # this prefetch avoids n calls where n is number of sequences #
  foreach my $seq ($feat->sequences({},{prefetch=>"sequence_type"})){
      print  $seq->sequence_data if $seq->sequence_type->name eq 'protein';
  }
  print  "\n\n";
}

fetches the protein sequence of any coding sequence (CDS) that has a feature name that starts with 'At1g' which would be any CDS on chromosome 1 of arabidopsis in our database.
The script consistently runs in 45 seconds on MySQL and in 29-31 seconds in Postgres. Other scripts seem to have about that difference--PostgreSQL finishes in about 60-70% of the time that the MySQL scripts do. Or, more dramatically: MySQL is 50% slower. That's pretty good for no change in the API, and all tables have the same indexing and structure.
Postgres was set up as default except for these values in postgresql.conf
shared_buffers 40000
max_connections 200
work_mem 4096
effective_cache_size 10000
This was the most concise indication of values to set, though likely closer tuning could improve performance (ahem, suggestions welcome).
An added benefit, is now, we can push more work into the database server using PL/Perl or another PL language, which can further reduce network back and forth, and reduce the creation of perl objects when not necessary.

python jumble solver

2006-11-26T17:14:00.000-08:00

over thanksgiving, i figured it'd be a good hack to solve the jumble . Originally, i tried to do all permutation of letter orders in the word but got confused and just did letter frequency:

import sys

word = len(sys.argv) > 1 and (sys.argv[1]).strip("\n").lower() or "egaugnal"
words = [w.strip("\n").lower() for w in open('/usr/share/dict/web2') if len(w) == len(word)+1]

def lfreq(w):
wfreq = {}
for letter in w: wfreq[letter] = letter in wfreq and wfreq[letter]+1 or 1
return wfreq

wfreq = lfreq(word)
match = [w for w in words if lfreq(w) == wfreq]

print match

that will print all matches for the word send in on the command line, assuming your dictionary file is in /usr/share/dict/web2 :-( . I tried for a while to do a recursive permute function which would take a word or list of letters and return all possible permutations: permute('abc') -> ['abc','acb','bca',bac',cba',cab'], which would clearly be a faster solution. but, as mentioned, i got confused and lazy. and the script above runs in <>

simple AJAX

2006-11-18T11:36:00.000-08:00

this is the AJAX implementation i use:

function jfetch(url,t,o) {
  var req = jfetch.xhr();
  req.open("GET",url,true);
  req.onreadystatechange = function() {
    if(req.readyState == 4){
      var rsp = req.responseText;
      if(t.constructor == Function) return t.apply(o,[rsp]);
      t = document.getElementById(t);
      t[t.value ==undefined ? 'innerHTML': 'value'] = rsp;
      req = null;
    }
  };
  req.send(null);
}
jfetch.xhr =
    (window.ActiveXObject)
   ? function(){ return new ActiveXObject("Microsoft.XMLHTTP"); }
   : function(){ return new XMLHttpRequest()};

it's short, it only checks for the transport (ActiveX or XHR) once, and it takes either an element id or a call back function. and, i can understand it.