(from a tea-room discussion on current selective pressures on Homo sapiens)
Friday, 4 July 2014
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:
Or to be mean:
(Ed: for those who don't know Perl this looks like)
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.
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:
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:
(& 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?
Aaaarrrrgggghhhh!
OK smug bioinformatician, how did the computer geeks go?
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.
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.
Monday, 7 April 2014
Genomics is a 4 letter word
Having trouble coming up with names for your genomics or genetics project? Perhaps you could consider an acronym using all of the DNA nucleotides? Act fast, as supplies are running out!
| Acronym | Name |
|---|---|
| ACGT | ACGT-a comparative genomics tool |
| ACTG | Now available! |
| AGCT | Applied Genomics & Cancer Therapeutics |
| AGTC | Advanced Genetic Technologies Center |
| ATCG | Now available! |
| ATGC | Alignable Tight Genomic Clusters |
| CAGT | Center for Applied Genetics and Technology |
| CGTA | Now available! |
| CTAG | Canadian Triticum Advancement Through Genomics |
| CTGA | Catalogue of Transmission Genetics in Arabs |
| GACT | Genetic Architecture of Complex Traits |
| GATC | GATC Biotech |
| GCAT | Genomics, Computational Biology and Technology Study Section |
| GCTA | Genome-wide Complex Trait Analysis |
| GTAC | Gene Technology Access Centre |
| GTCA | Now available! |
| TACG | Now available! |
| TAGC | The Applied Genomics Core |
| TCAG | The Centre for Applied Genomics |
| TCGA | The Cancer Genome Atlas |
| TGAC | The Genome Analysis Centre |
| TGCA | Now available! |
Feel free to make bets below for when they're all gone.
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?
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).
Create a new column based on calculations from other columns
Pull out a subset of columns:
Save to CSV:
Full diff and example project.
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.
Subscribe to:
Posts (Atom)