Sunday, 29 June 2014

Random bases and Infinite Monkey Theorem

An email conversation from work....

A:  We needed a big file,  do you have some test data?

B: cat /dev/urandom > random_data.txt

A: You know about Infinite Monkey Theorem, right? A monkey typing forever will eventually produce the works of William Shakespeare? I think you've just automated that monkey.

But I was thinking of something with somewhat smaller entropy... The tragedy would be if we couldn’t notice a difference...

B:  OK what about:

perl -e "while(1) { print(('G','A','T','C')[rand 4]); }" > random_bases.fasta

Or to be mean:


perl -e "while(1) { print(('G','A','T','C')[rand 4] x rand 10); }" > ion_torrent.fasta
 
(Ed: for those who don't know Perl this looks like)
CCCCTTTTTTTTTTTTTGGGTTGGGGGGGGTTTTTTTTTT
TTTTCCCCCCCCTTTTTTTAAAAGGGGGGGAAGGGGGTTT
TTTTAAAAAAAGGGAAAAATTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTGGAAAAAAGGGGGGAAAAAAAAATTTTTTT
CCGGCCCCTTTTTCCCCCCCCCGGGGGGGGGTTTTTTAAA

A: Actually it looks like we need a corollary, Infinite Monkey Genome Theorem: if you removed all the keys from the monkey's keyboard except  'G', 'A', 'T' and 'C', they will eventually type out the complete genome of William Shakespeare.

Monday, 2 June 2014

My aching chromosomes

Before you can do anything in genomics, you need a way to keep track of DNA locations - for instance where a gene  or a sequence starts and stops.

Nature packed billions of bases of DNA into separate chromosomes, somewhat like this:


chromosomes = {'chromosome 1' : np.array(249250621),
               'chromosome 2' : np.array(243199373) } # etc

And so a pretty natural system drops out, allowing us to refer to a genomic location via: "chromosome, start, end, strand".

Bioinformaticians can download multi-billion dollar text files for free, and they look like this:

head iGenomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa
>chrM
GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCAT
TTGGTATTTTCGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTG
GAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCATT
CTATTATTTATCGCACCTACGTTCAATATTACAGGCGAACATACCTACTA

(& on & on & on)

So if you look for "chromosome 1, 1234567-1234599 +" it's "GCCGGAGCCAGCCCCTCACAGCCGCATCAGCA" in hg19. This is all pretty cool, but nature didn't name them chromosome 1 and chromosome 2, that was up to humans. The humans decided to name them by size, with chromosome 1 being the biggest.

The biologists had one job: arrange the chromosomes in order, and label them 1 to N

The bioinformaticians had one job: take the chromosome names, and put them into text files so that people can share information.

How did the biologists do?

chromosome_sizes.py iGenomes/Homo_sapiens/UCSC/hg19/Sequence/WholeGenomeFasta/genome.fa.fai

chr20 longer than chr19 by 3896537 bases!
chr22 longer than chr21 by 3174671 bases!

Aaaarrrrgggghhhh!



OK smug bioinformatician, how did the computer geeks go?

chromosome_sizes.py iGenomes/Homo_sapiens/Ensembl/GRCh37/Sequence/WholeGenomeFasta/genome.fa.fai

20 longer than 19 by 3896537 bases!
22 longer than 21 by 3174671 bases!

Where's the 'chr'? Is chromosome 1 called 'chr1' or '1'?

Should the chromosome labelled "20" really be called "chr20", "19" or "chr19"? Which one is it? Aaaarrrrgggghhhh!

Chr or no chr... it's like picking which side of the road to drive on. It's completely arbitrary... and doesn't matter, so long as you are CONSISTENT!

"I've always driven on this side of the road!"
"Well I've always driven on this side of the road!"

Then they build a bridge and all of a sudden there are car smashes and dead bodies everywhere.

Except in bioinformatics we have misses instead of collisions, and a huge waste of people's time.

Wednesday, 26 March 2014

A ChIPseq Rorschach Test

Two people can look at the same picture and see two different things. Some psychologists think this gives a deep insight into a person's psyche.

I see histone H3 Lys4 trimethylation peaks relative to transcription start sites, what do you see?

I mapped a sample to the wrong genome...



Thursday, 13 February 2014

Refactor to Pandas example

A lot of bioinformatics involves loading and manipulating data which often comes in spreadsheet or similar data files (csv).

Pandas and Numpy are Python libraries which allow you to efficiently manipulate data. I recently gave a talk giving an example of common data analysis tasks (filter rows, calculations, sort by column etc) written first in straight Python, and then refactored to use Pandas (which has similar Dataframes to R).
Benchmarking the initial version:
In [1]: %timeit do_experiment(SIZE)
1 loops, best of 3: 1.86 s per loop
The refactored version is over 14x faster:
In [1]: %timeit do_experiment(SIZE)
10 loops, best of 3: 129 ms per loop
Besides performance improvements, I find Pandas code much simpler to read. With less boiler plate, the important operations are more clear.
Create a table with columns "one", "two", with a million random rows.  

Create a new column based on calculations from other columns


Pull out a subset of columns:


 Save to CSV:


Full diff and example project.

Hardest to Google Bioinformatics Software Awards

In the old days, programmers walked miles through snow between keypresses, and so brevity was key. They had short names - "cal" for calendar, "ed" for editor, "word" for word processor. Life was simpler back then, but life was good.

You can usually drop vowels in English (strt tlkng lke ths) or the middle of words (SxxxT TxxxxxG LxxE TxxS) and so "copy" becomes "cp", "make directory" - "mkdir" and "change group" is "chrgrp".

This worked for a while, but coders kept cranking out new code, doing ever more clever, but hard to describe things. Long before I was born, most of the good names were gone.

After a certain amount of time, a name becomes obscure enough to reuse, but this is not clear cut. And so we needed to keep coming up with more and more names for things, heading out into ever more difficult and obscure territory.

There are many complicated rules about naming, and it is very hard to pick a good name. You can't be too similar to something else, but you need to try and strive for memorability, pronunciation, efficiency and googlability. If shortened words aren't available, perhaps you can use acronyms: "awk" and "gatk" work great, but not everything does. I guess you have to write it out and say it a few times, to try it out.

The next step is creating a word, or borrowing an obscure one from somewhere else. When you invent one it still has to "look like a word". I can't really think of a compact way to describe this, but you know it when you see it. This is language and culturally specific - "Häagen-Dazs" (ice cream) was an invented name for a trademark - and it really does sounds like it's from a foreign country that makes great ice cream.

Two big winners via this strategy are Kleenex and Google. They also won big by being first to market or dominating a market. Nobody had seen the word "Kleenex" before, nor I guess seen paper so thin and cheap you could blow snot into it, and the two new things became associated together. I think you are supposed to fight to protect your trademark, but for a company this is the ultimate win, as their name becomes generic for a whole product.

While new names may sound funny, the borrowing or re-approriate strategy has a danger - if the word you pick isn't obscure enough, or your tool popularity / word obscurity balance isn't right - the name will backfire and you'll be lost in un-googlable waters.

It might seem that a new field gets a blank slate, without historical cruft. But Bioinformatics came from a mish-mash of fields - molecular biology, Unix, statistics, programming, maths - and it seems we imported all of the names and so we were born saturated.

One of the oldest and most important pieces of Bioinformatics software was called BLAST (Basic Local Alignment Search). The name is cool, and there weren't many tools available then, so it worked really well. It's the top hit for BLAST.

The trouble is, we can't all be BLAST. Skim a few bioinformatics journals and you'll see a stream of new software, with everyone trying to capture a name just like BLAST. In no way am I disparaging the fine software - the programs may be cool, but the names are too much of a stretch - they'll never going to claim top spot.

But there is a place where they can claim top spot - in the annual bad naming awards. The runner's up are:
  • Macs - ChIPseq software, but for some reason I feel like listening to my iPod while eating a hamburger.
  • DAVID - (webapp for bioinformatics annotation) - this is really annoying when you are called David, and do bioinformatics annotation.
The winner was revealed a few weeks ago, as my co-worker and I sat with monitored univerity proxies between our lab and the internet. We were investigating whether a sequencing run could be saved by calibrating the basecalling on the 3' end of reads rather than the 5' (which had low complexity) and found someone who had written scripts for this very purpose.

We were impressed and looked more into it. After telling me about this tool, she added "bareback is a very difficult word to Google, by the way".