Yes! Why should the evolution people have all the fun with their blog carnival? (After all, it is only a theory.) It’s time for bioinformaticians to show what we are made of, and to have a carnival of our own. Bio::blogs had a good run some time ago. I decided to reconnect what is hopefully now a larger and more networked community under the title of Bioinformatics Blog Carnival.

Northampton: monument to Francis Crick, born in Northampton. Copyright Ian Rob, licensed for reuse.

Programming

For many new bioinformatic programmers, there is a question of which Bio* package to choose from. The Bio* packages (Biopython, BioPerl, BioJava and BioRuby) are open source licensed packages which are used heavily in bioinformatics coding. They contain parsers for bioinformatic data file, file format converters, SQL interfaces and other goodies that make a biohacker’s life manageable. Walter Jessen presents a rapid-fire no-nonsense review of Open Source Programming with Bio* Libraries on Expressing Scientific Insight. Finally, if you would like to do everything wrong, Manuel Corpas presents 10 Sarcastic Rules on How to Be a Bioinformatician posted at Manuel Corpas’ Blog.

Databases

Morgan Langille talks about BioTorrents – a file sharing resource for scientists posted at his blog Beta Science. My take on BioTorrents is that it is a cool idea, but as most institutes block BitTorrent along with other peer-to-peer sharing utilities, I doubt there is a critical mass of feeders to make BioTorrents viable. Things can change though: institutes might decide to set up dedicated servers to Torrent scientific data, just like there are legitimate Torrent servers for Linux distros. George presents Bio-graphics, BioSQL and Rails part 1 posted at Biorelated. He talks about how to quickly add graphics support to a bioinformatics database rails application.

Brad Chapman’s Blue Collar Bioinformatics is a treasure-trove of useful bioinformatics methods. Brad is very thorough in his writing, and he covers a wide variety of topics. I particularly enjoyed reading about his adventures at the biohackathon 2010 in Tokyo, and the resulting Python query interface to BioGateway SPARQL endpoint and InterMine.

Genomics

Nick Loman from Pathogens: Genes and Genomes gives some excellent tips for de-novo genome assembly. He also gives the necessary scripts,which are great companions to Velvet, the popular short read assembler. Luke Jostins writes about AGBT: Speculating on Third Gen Tech posted at Genetic Inference, “An investigation into what data from various third generation sequencing technologies may look like.” SM presents Resources for Exome sequencing annotation posted at Organizing the Strands of Curiosity.

Structural biology

Sean Seaver talks about the retraction of several protein structure papers published by one group at the University of Alabama, Birmingham. Structuregate was reported both in the media and in the blogosphere. Sean walks us through how it might have been done, and how structural bioinformatics techniques found out the wrong structures in Origin and Orientation posted at P212121. Maria Hodges talks about how difficult it is to explain structural genomics to the man in the pub. Menachem Fromer presents Tradeoff between stability and multispecificity in the design of promiscuous proteins. Promiscuous proteins are proteins which bind different partners (ligands, or other proteins). However, the more partners they are able to bind, the less stable they are, as shown in a series of simulated evolution studies, summarized at Nir London’s Macromolecular Modeling Blog. OK, promiscuity and and a pub. There must be a joke somewhere there.

Credit: Iddo Friedberg & 3D chem for the Anthrax toxin lethal factor image

A common lament among bloggers and other enthusiastic adopters of Web 2.0 technology is the lack of mainstream uptake of these tools by active scientists. A recent report from University of California Berkeley confirmed this: “The advice given to pre-tenure scholars was consistent across all fields: focus on publishing in the right venues and avoid spending too much time on public engagement, committee work, writing op-ed pieces, developing websites, blogging, and other non-traditional forms of electronic dissemination (including online course activities)“. Maria Hodges argues that Web 2.0 thrives where journals don’t, and that the NMR community might be the first to reach the tipping point, where your career is harmed by not contributing. She talks of two wikis used by the NMR community. It is a small community, which has a need for sharing methods that are generally not publishable in peer-reviewed journals. The wiki venue makes for an ideal dissemination method for this community. On the subject of data sharing, and how it can backfire: here is an interesting connection between the history of the PDB, and the email affair known as ‘climategate’.

There are more microbial cells in our body than our own. Those microbes are not just passive hitchhikers or conversely, malicious agents of disease. They affect our well-being and health in a much broader spectrum than simply “bad” or “passive”. Among other things our gut microbes play an important role in digestion, have been linked to obesity, conditions as severe as certain inflammatory bowel syndromes or as relatively mild as traveler’s diarrhea. The microbes living on our skin also affect many things: from body odor and dandruff and acne to dermatitis and psoriasis. Also, being the most exposed microbial population means that they are themselves affected by our constant exposure to various agents,and their population varies by our own behavior –  such as how often, and with what, do we wash our hands: antibacterial soap has been named as a major culprit in the development of resistant bacterial strains. In all organs the native flora, relatively benign, protects us against colonization by more virulent bacteria.

Indeed, our gut, skin, mouth and genital microbiomes can be viewed as  additional organs in the way that they affect us. If you take a long course of powerful systemic antibiotics, the ensuing diarrhea and sometimes mouth thrush are the result of these “organs” — in your mouth and gut — being removed temporarily from your body.

Credit: David Gregory&Debbie Marshall. Wellcome Images images@wellcome.ac.uk

Today, another big step has been taken towards understanding the role of our microbiomes play. Just how big, in what direction, or what will be the consequences of this step is unclear, and will remain unclear for quite some time. A group of Chinese and European researchers have published the largest sequencing effort yet of gut bacteria. Their current yield, 576 Gigabases of DNA from the feces of 124 European people is considerably larger than the previous large effort in the US: 3GB from sequenced from 33 US and Japanese adults. Also, Qin and colleagues looked at obese, lean, healthy and sick (inflamatory bowel syndrome) individuals. They identified 1,000 to 1,100 different bacterial species, with about 160 different species per individual. Healthy individual’s bacterial population was markedly different that those with inflammatory bowel syndromes, and the populations of those with IBS differed depending on disease type.

What does this all mean? Well, like in the human genome project, it will take a while before  we understand not only what these data contain, but what are the limits of our ability to interpret them, and how best they can serve us. What we have right now is the equivalent of the outline of a newly discovered continent. It is up to many individual explorers to discover and chart the myriad living things in terra excreta. Which genes are most associated with obesity or with ulcerative colitis? How prevalent is gene transfer between gut bacteria, and how much of a role does this play in antibiotic resistance? Are there microbial species more prone to changes in their genomes than others? Are there metabolic pathways that are shared between different microbial species? Are some genes faster evolving than others, what would be the ecological role different species play in the gut? And how do different microbial populations ultimately affect their human hosts? Do different bacterial species have a preference with whom they share our guts? These broad questions can be broken down into individual questions relating to a lab’s pet genes, metabolic pathways, microbial species and metabolic conditions. There is a lot to sift through here and these data will keep us busy for years to come.

Pathogenic E. coli on the intestinal lining. Credit: S. Schuller.

In other blogs Carl Zimmer has a great post on how he, for one welcomes our bacterial overlords, while Ed Yong talks about the science of things to come.

Oh, and for those of you who, like me, can’t wait to plug the data into your favorite analysis pipeline, you can get the gene and assembled and annotated sequence data data from Peer Bork’s lab at the European Molecular Biology Laboratory, or the Beijing Genome Institute. The raw sequence data are available from the European Bioinformatics Institute under the accession ERA000116. (I can’t find the latter myself right now, hopefully they will show up in a couple of days). Update March 7, 2010: the sequences are now deposited at EBI. Also, on BioTorrents (hat tip to Morgan Langille).

Some interesting events cropped up recently. The Marco Island Advances in Genome Biology and Technology meeting was heavily tweeted and blogged about.  Pacific Biosciences unveiled their third generation sequencer. Ostensibly, it can sequence reads of 20,000 length, but the fraction of actual long reads in a run, and their quality is still a bit hazy. The most interested to me is the Ion Torrent. Being rather low on budget, this seems like the family budget car of high throughput sequencing: cheap, reliable, and does not offer more than I really need. $50,000 for a sequencer with $500 runs with 160MB/hr? Nice. Genetic Inference has a great summary of the various technologies presented.

Overall, we are starting to see a divergence in sequencing technologies, as each tech concentrates on having clearly defined advantages and potential applications that differ from all others. This means that the scientists themselves can more closely tailor their choice of tech to fit their situation. Are you a small lab that needs 10 high-quality genomes on a budget? Go to Complete. Want a cheap, fast machine for library validation? Use Ion Torrent. Setting up a pipeline for sequencing thousands of genomes? Go Illumina.

The review article on metagenomics I recently published in PLoS Computational Biology (yeah, yeah, shameless plug) already starts to feels somewhat outdated on the sequencing technology front.

Carnival of Evolution #21 the superstar edition is up: check it out. It’s a nice and detailed one,. Some posts I liked included talking about how human fingers evolved, and why it is important to consult evolutionary biologists when making decision about conservation.

An interesting email I got yesterday: PubGet, a search engine for PDFs of scientific articles, is no linked to PLoS. PubGet is a very useful service that gets  you the article PDF immediately, without going through he usual clickeroo via Google,  pubmed, publisher’s gateway, journal gateway and then squinting along the sidebar to find the PDF link. Nice to see that these two are teaming up.

Finally, two reminders. First, Ada Lovelace day, a blogging day celebrating the achievements of women in science and technology is coming up, March 24. Go ahead, pledge and blog! Second, the Byte Size Biology will be hosting a Carnival of Bioinformatics. Quite a few posts have been submitted already, please submit yours, deadline: March 9.

Byte Size Biology will be hosting the first edition of the bioinformatics blog carnival. All you bioinformatics bloggers, submit your entries by Mar 9, 2010 23:59:03  EST. Note the 3 second extension I have already given. There will be no more deadline extensions, I’ve been generous enough as it is. The carnival will be posted here by March 15, 2010.

Any blog posts that have to do with the computational aspects of: genomics, nextgen sequencing, sequence analysis, gene expression, systems biology, ontologies, databases, structural biology, metagenomics, phylogenetics, function prediction and I probably forgot a few other categories so don’t hold it against me, just submit. Early and often. Your own posts and others that you liked. Nothing too old please, 1/1/2009 and later. I reserve the right to be a less-than-benevolent dictator and screen out posts. This applies especially to commercial plugs with no other merits.

Please retweet, reblog, rebuzz and remember:  submit to the bioinformatics blog carnival!

If you have a cool logo for this carnival, email bioinfo.blog.carnival AT gmail.com OhOne and OhToo will pick the winning logo to be displayed on the carnival. Likewise, if you would like to host the next edition, let me know. Do not submit posts by email, only via blogcarnival.com.

Happy carny-ing.

As part of the process of manufacturing  a new car,  the designers will take the blueprints to the factory floor. There they will set up an experimental assembly line, tinkering with the manufacturing process of the prototype until it is ready for mass-production. Can we do the same with the machinery of life – the assembly of proteins? Can we set up an alternative assembly line for a new protein prototype — and then actually set up a working assembly line for the whole new protein?  A proof-of-concept has been published this week in Nature by Jason Chin’s group at the Medical Research Council Laboratory of Molecular Biology, Cambridge UK.

If there is a single common denominator to all life, it is the genetic code. All life is built around DNA encoding information for proteins  nucleotide triplets or codons. Since there are four types of nucleotides (A,T,G,C)  that are read in words of thee, there are 4³ = 64 possible codons: more than enough to encode for the 22 amino acids that make up proteins. There is nothing more basic and fundamental to life on Earth than the three-letter based genetic code.

Until now.

Chin’s group has created a four-nucleotide codon system.  It is not that the DNA is different: it is the way the cellular machinery decoding  RNA transcripts interprets the nucleotide sequence. Ribosomes –large RNA and protein complexes  which are the platform upon which messenger RNA is read and decoded — are set to serve up messenger RNA three nucleotides at a time. (Messenger RNA or mRNA is a transcript of the DNA which is carried to the ribosome.)  Transfer RNA or tRNA is a short RNA molecule that shuttles the proper amino acid to the ribosome, but will only attach if the proper codon is served up by the ribosome. The whole protein synthesis “assembly line” looks something like this:

Protein synthesis. Credit: Wikimedia Commons.

To change the interpretation of the genetic code from three lettered words  to four, Chin and his colleagues had to make new ribosomes, and new tRNAs.  To create these new ribosomes, they designed orthogonal ribosomes, or o-ribosomes. O-ribosomes are genes inserted to produce extra ribosomes that operate in the cell alongside the regular ribosomes. The cell functions because it has the regular ribosomes to maintain its viability. The ribosomal RNA in the o-ribosomes is free to be mutated to create new unnatural traits: in this case, the ability to serve as a platform read four-letter codons. They selected for Escherichia coli bacterial cells that expressed a o-ribosomes which translated a four-letter codon in a gene, which would otherwise go untranslated by the regular ribosome. The gene gives the bacterial cells resistance to the antibiotic chloramphenicol. So cells that survive a dosage of chloramphenicol are those which have functioning o-ribosomes, as they have the chloramphenicol resistance gene that is being translated by the o-ribosomes.

They also needed to create new tRNAs that have an four-nucleotide anticodon (the part that complementarily binds to the messenger RNA –  see figure above.)  So the surviving E. coli cells have a population of working o-ribosomes, regular ribosomes, modified tRNA (with a  four-letter anticodon) and regular tRNA.

Then they took their work a step further. Each three-letter tRNA carries a specific amino-acid, depending on its anticodon. Thus tRNA_AAG will always have a phenylalanine attached, because CTT (the complement of AAG on the messenger RNA) codes for phenylalanine. If you start messing with that, the translation machinery will produce non-functional proteins, which will probably kill the cells pretty quick. But with the orthogonal 4-letter code machinery, that is not really a problem: the orthogonal machinery operates alongside the normal one. Also, there are no amino acids naturally assigned to any four letter code, because this code does not appear in nature in the first place! So Chin’s lab assigned an unnatural amino acids to a four-letter code. The non-naturally occurring p-azido-l-phenylalanine amino acid was assigned to tRNA_UCCU. They then showed that the whole alternative translational machinery worked by synthesizing a mutant of the protein calmodulin which used this amino-acid in its structure.

Why do it? Well, personally I don’t see the need for justification: just being able to do it is so cool!  But seriously: think of the ability to design proteins from up to 4⁴=256 different amino acids other than the 22 we have now.  The possibilities of tinkering with existing proteins using this orthogonal, four-letter based machinery are huge. The other benefit of this orthogonal synthesis setup is the ability to control this orthogonal translational machinery: because it does not use the three-letter vocabulary, this orthogonal machinery would be much easier to manipulate, tinker with and switch on and off without getting in the way of regular cellular translational machinery. The analogy to a car assembly line breaks here. It is as if two different models are being assembled on the same line just by using different robots. The better analogy is for a program source code to be read by two different compilers, each producing a different program. Awesome.

Seriously, this is what I first thought when I saw the cover of this week’s Nature, and the associated drawings  in the press.  The dude’s haircut seems like it was even bad in the ’80s… 2080 BCE that is, which is when his body is dated. Approximately.

Inuk, artist's impression

A large group of researchers were involved in analyzing DNA from a 4,000 year old man found in Greenland. There is an interesting story on how the sample was obtained. Eske Villerslev spent 2 months looking for human remains in archaeological sites in Greenland, to no avail. When complaining about it to a friend, he learned that that friend’s father had already given hair samples to the national museum of Denmark… back  in 1986. Serious slap on the forehead!

The preserved DNA was of great quality. Also, this is 2010: we have a plethora of genotypic data that lets us associate a genome with phenotypes, and geographical origin. Inuk, as he was dubbed, shares many traits with east Asians, and much less with contemporary North American Inuit or other Native Americans. This means that Inuk and his people migrated fairly recently to Greenland from Siberia and across North America. Inuk and his people, the (now extinct) Saqqaq have no relationship to contemporary Native Americans. This was also shown in an earlier study of the mitochondrial DNA.

What do we know about Inuk himself? He had blood type A+; thick black hair and very likely a receding hairline due to a baldness allele; brown eyes; shovel-graded front teeth (East  Asian characteristic); dry type earwax… before you laugh that is actually, another characteristic of Asian populations on the Siberian and Chinese east coast. He also had alleles associates with human adaptation to cold temperatures.

And a bad haircut; but great press nevertheless (just got to Google News and search for  “Greenland Genome” – see, for example, the NYTimes.)

This post has been Slashdotted. Exercise extreme caution.

Gilbert MT, Kivisild T, Grønnow B, Andersen PK, Metspalu E, Reidla M, Tamm E, Axelsson E, Götherström A, Campos PF, Rasmussen M, Metspalu M, Higham TF, Schwenninger JL, Nathan R, De Hoog CJ, Koch A, Møller LN, Andreasen C, Meldgaard M, Villems R, Bendixen C, & Willerslev E (2008). Paleo-Eskimo mtDNA genome reveals matrilineal discontinuity in Greenland. Science (New York, N.Y.), 320 (5884), 1787-9 PMID: 18511654

If the title of this post makes you cringe, then you belong to a minority of people who realize why the phrase “highly evolved” is so wrong. Unfortunately, “highly evolved” (as an absolute term) and “more evolved” (as a comparative term) seem  to be used all-too frequently.  They are uttered not only by non-scientists and non-biologists but even by scientists who should know better. Even when they catch themselves after blurting out “highly evolved” in a conversation (or, more embarrassingly, in a lecture), the damage is done. Yet another Freudian (Darwinian?) slip that tells of a fundamentally bad grasp on evolution. And, yes, I know, this topic has been written about by many of my betters, who are vastly more evolved better writers than I, with much better breadth and depth of knowledge of evolution, and a reach to a much wider audience.

More evolved? Source: Wikimedia commons, public domain.

So why am I writing about it? Well, this is my blog and ranting in it is my prerogative. And despite the Richard Dawkinses and Steven Jay Goulds of this world, the use of this phrase still persists. So it is up to us foot soldiers of the blogging community to do our own modest bit. If I prevent any of my six readers from being tempted to utter this phrase the next time it is (wrongly) deemed appropriate, then I have done my bit.

Why is this “highly evolved” used so much? And why is it wrong?

Consider the sponge, and then consider Albert Einstein. There are certain traits that Einstein had, that a sponge does not. We deem these traits to be of merit. Einstein developed a fundamental theory in physics. He  played the violin. He  ate with a knife and fork, had binocular color vision, opposable thumbs and he cultivated his facial hair in the form of a mustache.

A sponge… well, to be brief, does not have all those qualities we hold in such high merit. It kinda sits there at the bottom of the shallow ocean, flopping about, filter feeding, pooping and apparently not much else. Clearly, there are qualities to Einstein that make him more interesting than the sponge.

Less evolved? Source: Wikimedia Commons

Einstein seems, intuitively, to be more complex than a sponge, and that complexity can be quantified directly, in many ways. Actually, this is a pretty contentious point by itself: can we speak of organism complexity? Can we quantify the complexity of an organism and compare between different species? And what exactly would the complexity metric we choose tell us?

But let us assume, for argument’s sake, that our intuition that Einstein is more complex than a sponge is correct. For example, we can imagine a measure derived from the diversity and number of cells. Obviously there are more cell types in Einstein than in a sponge. Does that mean he is also more evolved? Are humans a more  evolved than sponges? Chimps? After all, did life not start 3.85 billion years ago as simple and over time became more complex? Progressing, as it were from simple unicellular bacteria through more complex sponges all the way to the crowning achievement of humans? Had humans not, in a sort of (alas, Pyrrhic) victory, mastered the Earth and competed with many of earth’s species to the latter’s extinction? Isn’t competition what evolution is all about? And isn’t human victory a direct result of human complexity making humans “more evolved”? So isn’t “complexity” an end product of evolution, the more complex you are the more successful you are, and the more evolved you are?

No, no, no, no, no, and no.

The reason for this series of compounding errors is the mistaken notion that evolution by natural selection is a progression resulting in a production of increasingly complex life.  Evolution is not goal oriented, and there is no teleology involved. The increasing complexity of organisms along time may seem to involve a  progressive process, but there is none. It is a “statistical illusion”.  What do I mean by that? Well, life did start out in less complex forms, that became more complex. But the less complex forms remained as well. Thus the distribution of complexity increased over time, but there is no directionality towards progress: the less complex life remained around as well. But over 3.85B years, complexity has had a chance to manifest itself in life, as natural selection favored some initial complexities, and those extended to become even more complex. Yes, we can trace a direct route from the first multicellular organisms, through sponges, invertebrates, vertebrates. But humans, chimps, sponges and bacteria living on Earth today are the result of exactly the same selective forces that have shaped life since  it crawled out of an underwater volcano, or wherever. Complexity emerged over time, and is still emerging. But complex organisms are being added to the pool of life, rather than replacing the simple organisms. The result is an increase of a distribution of complexity levels, not the moving of an entire curve of complexity rightwards.

Apparent progress due to a to a 'wall' restricting where random change can take things. Adapted from SJ Gould. Reproduced under CC from talkorigins.org

The point I am trying to make is that humans may be more complex than sponges, but we are not “more evolved” nor are we “highly evolved”.  There is no progressive process, and all of life on earth is the result of the same 3.85B years of selective pressures.

For a really good historical overview of teleological, or purpose-driven, thought in evolution, look to talkorigins.org.

Few know that Einstein was teaching evolution at Princeton. Physics was just a cover.

All of this does not mean that Highly Evolved by The Vines is not a kick-ass song. Listening to it is also  a good way to get the rage from hearing “highly evolved” out of your system.  Note the low complexity of the video:

In celebration of the biohackathon happening now in Tokyo, I am putting up a script that is oddly missing from many bioinformatic packages: extracting intergenic regions. This one was written together with my student, Ian. As for the biohackathon itself, I’m not there, but I am following the tweets and  Brad Chapman’s excellent posts:

• Day 1
• Day 2
• Day 3
• Day 4

About intergenic regions: intergenic regions are as interesting and sometimes even more interesting than the genes themselves: when you are interested in promoters, transcription factor binding sites or almost any other transcription regulation mechanism. Here’s a simple script to find intergenic regions. It reads a genbank formatted file and uses the information there to extract the intergenic regions. The sequences are written to a FASTA file.

#!/usr/bin/env python
import sys
import Bio
from Bio import SeqIO, SeqFeature
from Bio.SeqRecord import SeqRecord
import os

# Copyright(C) 2009 Iddo Friedberg & Ian MC Fleming
# Released under Biopython license. http://www.biopython.org/DIST/LICENSE
# Do not remove this comment
def get_interregions(genbank_path,intergene_length=1):
    seq_record = SeqIO.parse(open(genbank_path), "genbank").next()
    cds_list_plus = []
    cds_list_minus = []
    intergenic_records = []
    # Loop over the genome file, get the CDS features on each of the strands
    for feature in seq_record.features:
        if feature.type == 'CDS':
            mystart = feature.location._start.position
            myend = feature.location._end.position
            if feature.strand == -1:
                cds_list_minus.append((mystart,myend))
            elif feature.strand == 1:
                cds_list_plus.append((mystart,myend))
            else:
                sys.stderr.write("No strand indicated %d-%d. Assuming +\n" %
                                  (mystart, myend))
                cds_list_plus.append((mystart,myend,1))

    for i,pospair in enumerate(cds_list_plus[1:]):
        # Compare current start position to previous end position
        last_end = cds_list_plus[i][1]
        this_start = pospair[0]
        strand = pospair[2]
        if this_start - last_end >= intergene_length:
            intergene_seq = seq_record.seq[last_end:this_start]
            strand_string = "+"
            intergenic_records.append(
                  SeqRecord(intergene_seq,id="%s-ign-%d" % (seq_record.name,i),
                  description="%s %d-%d %s" % (seq_record.name, last_end+1,
                                                        this_start,strand_string)))
    for i,pospair in enumerate(cds_list_minus[1:]):
        last_end = cds_list_minus[i][1]
        this_start = pospair[0]
        strand = pospair[2]
        if this_start - last_end >= intergene_length:
            intergene_seq = seq_record.seq[last_end:this_start]
            strand_string = "-"
            intergenic_records.append(
                  SeqRecord(intergene_seq,id="%s-ign-%d" % (seq_record.name,i),
                  description="%s %d-%d %s" % (seq_record.name, last_end+1,
                                                        this_start,strand_string)))
    outpath = os.path.splitext(os.path.basename(genbank_path))[0] + "_ign.fasta"
    SeqIO.write(intergenic_records, open(outpath,"w"), "fasta")

if __name__ == '__main__':
    if len(sys.argv) == 2:
         get_interregions(sys.argv[1])
    elif len(sys.argv) == 3:
         get_interregions(sys.argv[1],int(sys.argv[2]))
    else:
         print "Usage: get_intergenic.py gb_file [intergenic_length]"
         sys.exit(0)

What are we seeing here?

Lines 11-16 are the preamble: we read the GenBank file using Biopython’s genbank parser in line 12.  Beacuse we expect a genome file, which contains a single record, this is a one-time read. Note that this is a rate limiting step, and can take a couple of seconds. Took me ~2secs to read the full E. coli genome on my Linux box.  We prepare one list for the + strand intergenic regions (13), another one for the minus strand intergenic regions (14) and one for all the records (line 15).

The rest of the code are three loop blocks: lines 16-28 I loop over the genbank features, extracting the coordinated of the genes themselves. Line 32-41 I find the intergenic regions on the + strand. Lines 42-52 I do the same for the “-” strand.

Now for a philosophical interlude: although there is a way to read all the intergenic regions in a single pass, I subscribe to the “code simple” doctrine of research software writing. Code performance optimization is a low priority for me. I’d much rather have something that is simple to write,read and modify. I also don’t want to spend too much time coding and elegant script for elegance’s sake, especially if I may not use it too much. Historically, scientific code written for research is mostly extinct: thrown away after a short lived hypothesis was tested and ended its days. Research coding is mostly throwaway glue code. Very rarely it matures into a product. Then, and only then, can you apply all those fine software engineering you learned in college. Before that, write fast and simple.

But I digress. Line 17 loops over the features in the genome file. Line 18 we identify if it is a coding sequence (CDS). If so, we identify the start position, and position and the strand the CDS is on. The list cds_list_minus is a list of 2-tuples. Each 2-tuple is the start and end positions of a CDS on the minus strand. (If you would like to go over the genes, as opposed to coding sequences, change line 18 to:

if feature.type == 'gene':

(or better yet, pass an argument that defines it.)

cds_list_plus, is, yes, the same as cds_list_minus, only for the plus strand (line 24).

Sometime, a CDS does not contain information on which strand it is. With genome files, that is usually the case with single stranded viral genomes. Therefore, we put in the default assumption that if there is no strand indication, then the feature is is on the plus strand. We generate a warning message nevertheless (lines 25-28).

Lines 30-41 we loop over the plus  strand list, and identify the coordinates between the genes. Python’s enumerate function is very useful here. The enumerate function allows us to iterate over a list, but at the same time keep track of which index we are in when looping over the list. So in line 30, pospair receives the start and end coordinates of a CDS as a 2-tuple, while i receives the actual number if the index in the plus strand CDS list. In that way, we can look back to the previous list member, find the coordinates where that CDS ends, and where the current CDS begins. The two coordinates make up the beginning and end of the intergenic regions between those two genes on that strand. In line 35 we check if the intergenic region length is equal to or larger than a threshold: suppose we are only interested in those intergenic regions that are longer than 100 bases? (The default value is 1, see line 11.) In lines 38-39 we build a biopython sequence object that contains an informative header, and the sequence of the intergenic region. The description which goes in the sequence header contains the start and end coordinates of the intergenic region, and the locus ID of the CDS directly downstream from it. The sequence object is appended to a list, which will eventually get written (lines 40-41).

Lines 42-52 are a repeat of lines 30-41, only for the minus strand.  Lines 53 & 54: the list that contains all the intergenic region sequence objects gets written to its own fasta file.

Finally, line 56-63 are boilerplate code, that make this script runnable from the command line. Have fun looking at intergenic regions. Let me know of you find something interesting.

A  small spike on my blog traffic yesterday led me to look for the source via Google Analytics. (If you are a blogger, you should really use this tool, lots of useful traffic information.) Seems like most of the traffic came from the page of a high school science teacher at Badin High School in Hamilton, OH. Apparently the students were to be quizzed today on two of my posts about endosymbiosis (and one from 80Beats; I’m in good company.) So they were very busy Sunday. It’s encouraging to know that some of my posts are accessible enough for high school science. Finally, quite a few Miami students come from Hamilton (we’re close). So I might see some of them next year.

Muahahaha!

Pharmaceutical companies are always on the lookout for secondary drug targets. After all, if you invest billions developing a single drug, you would be more than happy to sell it as a treatment for two, three, or more different ailments. Sildenafil citrate was developed to treat angina and hypertension, but during phase I clinical trials, it was found that Sildenafil induces penile erections. The drug was branded Viagra, and the rest is history. Eflornithine, an anti-cancer drug, is also effective against the agent of African sleeping sickness, Trypanosoma brucei. African Sleeping Sickness is known as a “neglected disease”, for which drug development is not profitable and therefore not a priority. However, having a drug already on hand makes it easier to distribute in affected areas, since the R&D costs are recovered elsewhere.

Another example of polypharmacology is a drug that binds to multiple targets in the human body. This could be used for overcoming drug resistance, a known problem with cancer. Cancer tumors often develop a resistance to anti-cancer drugs by simple natural selection: the protein that the drug binds to mutates, and no longer binds the drug. However, if the drug acts by binding redundantly to several proteins, it would be more effective, since several mutations would be required to effect drug resistance.

Another important polypharmacological consideration  is toxicity. If a drug binds to one protein, its drug target, it may also bind to another one which it should not bind to as it disrupts the normal functions and the health of the patient. If the side effects outweigh the cure, the drug is no good.

How one drug (cyan) can bind to two different proteins with different overall shapes (pink and green), but with similar binding sites

Because it can either increase  profits, or conversely derail a whole process of drug development, predicting polypharmacophoric effects is very much something drug developers want. A study published yesterday in PLoS Computational Biology by Jacob Durrant and his colleagues suggest a way  bioinformatics and theoretical biophysics can help in identifying multiple drug targets. Durrant’s goal was simple: given the molecular structure of a candidate drug, which proteins are expected to bind it? The strategy this group took is a combination of bioinformatic and experimental screening.

Finding additional drug targets in four  steps

Reproduced under CC license from doi:10.1371/journal.pcbi.1000648.g001. Click for original image.

Step 1 (A-C in the figure above): identify the known target protein. Now pick all the protein structures that are not similar to it. Why those that are not similar? Similar proteins could be potential drug targets, since they have a similar shape to the known target protein. But here they are interested in finding targets from proteins that are of a different shape, and have no homology to the known target protein: secondary targets.

Step 2 (D): take this set of dissimilar proteins, and look for binding site similarities. Binding sites are clefts in the protein that may bind drugs. If those clefts are similar in shape to the cleft in the known target proteins, they may bind the drug. Leave only those non-homologous proteins that have similar binding sites (D in the figure)

Step 3 (E): Now add all the proteins that are homologous to the set generated in step 2. This increases the number of possible targets to homologs of the proteins that were initially selected only for binding site similarity.

Step 4 (F): take the drug molecule, and try to dock it to the various protein structures. Rank the druggability of each protein according to the score provided by the drug docking software (Autodock).

Experimental verification

Now for the cool part. Durrant and colleagues tested  the computational prediction in the lab, using the compound NSC-45208. NSC-45208 inhibits a protein responsible for RNA processing in Trypanosoma brucei. So we know it is a potential drug against African sleeping sickness. What else is it good for?

(NSC-45208), 4,5-dihydroxy-3-(1-naphthyldiazenyl)-2,7 -naphthalenedisulfonic acid, a recently discovered inhibitor of T. brucei RNA editing ligase 1 (TbREL1)

“The predicted secondary targets that gave the best docking scores, H. sapiens mitochondrial 2-enoyl thioester reductase (HsETR1), T. brucei UDP-galactose 4′ epimerase (TbGalE), H. sapiens phosphodiesterase 9A (HsPDE9A2), and Streptococcus pneumoniae teichoic acid phosphorylcholine esterase (SpPce), were subsequently tested experimentally.”

Durrant and his colleagues tested their predictions that NSC-45208 also binds to two human proteins (HsETR1 and HsPDE9A2), one additional Trypansome protein (TbGalE), and a bacterial (Streptococcus pneumoniae) protein (SpPce). Not only binds, but also inhibits their enzymatic activity. Their predictions worked well on the top two predicted targets: NSC-45308 inhibited the enzymatic activity of HsETR1and TbGalE, but HsPDE9A2 and SpPce were not affected by the drug.

At first blush, this does not seem to be much of a batting average: two positives and two false positives. But we have to remember that the experimental verification of these predictions — even if you have a good enzymatic assay to check predictions– can be very time- and resource consuming. An exhaustive test of the predictions for several compounds and many target enzymes is still not possible. As initial proof-of-principle, this work goes much farther than most other joint experimental / computational works I have read. The authors also go into lengthy and detailed discussions on limitations and improvements, as well as on another form of non-specific inhibition of the secondary targets, makes for a really interesting read on polypharmacological considerations in drug screening.

Seems like I had a couple of taxonomic posts recently.  So BsB would like to know what is your favorite Kingdom of life. I used the 6 kingdom system just to make things slightly more interesting. Vote in the pollbox on the righthand sidebar. You only have yourself to blame if your favorite kingdom is not represented adequately.

[View with PicLens]

Despite our best attempts to remove species from the face of the Earth, there is still quite a bit of life out there and it is still quite diverse. Also, there are still quite a few people who want to document, describe and make the rest of us aware of the magnitude and diversity of life. In the time-honored tradition of natural history they take the tools of their trade to the field to sample, collect, identify, tag, photograph, film, mount and otherwise understand and advertise Nature’s progeny.  There are probably thousands of web sites dedicated to taxonomy of families, genuses genera and species. Trouble is, they are all “out there”, and require quite a bit of searching to find those that are of interest to you. And what if you just want to browse, the Internet equivalent of going through a natural history museum?

Museums were historically the place where natural historians brought their findings from far away places to be studied and displayed (if lucky) or simply archived. If the finding was still alive, it would go to the zoo or the botanical garden.  And people would go to the museum, zoo or garden to see the  latest wonders of life.

[View with PicLens]

(From Life Through a Lens: the Natural History museum (of London) 1880-1950)

Virtually every biology lab with some taxonomic interest in  a given species, genus or family would have such a site. That’s a lot of taxonomic information , but is not very accessible or organized! What if all these sites could be rounded up together, connected in some contextual fashion, some semblance of standardization? Then we could have an online Natural History museum, where the displays accumulate, are constantly updated, and we can walk through virtual halls looking for the plant or animal of our fancy.

Scratchpads is a cool way of getting to that online Natural History museum. Actually, it is hosted in the great-grandaddy of brick-and-mortar natural history museums, the one in London.

Scratchpads are an easy to use, social networking application that enable communities of researchers to manage, share and publish taxonomic data online. Sites are hosted at the Natural History Museum London, and offered free to any scientist that completes an online registration form. Key features of the Scratchpads include tools to manage:

	Classifications		Phylogenies
	Bibliographies		Documents
	Image galleries		Custom data
	Specimen records		Maps

Users control who has access to content, which is published on the site under Creative Commons (by-nc-sa) license.

Data added to a Scratchpad are automatically classified and grouped around a taxonomy that is supplied by the users or imported from EOL. This is optionally supplemented with information from high quality web accessible databases, to automatic construct content rich web pages about any documented taxon. Currently these sources include Genbank, GBIF, Biodiversity Heritage Library, Yahoo! Images, flickr, Google Scholar and Wikipedia.

When you start a Scratchpad site, you are given a template on which you can build your taxonomy site. the template includes a tree-building widget, space for image galleries, a blog area, a bibliography area etc. You insert the content, and your site links back to the scratchpad mother site.  The Scratchpads linked from the mother site are in various stages of construction, which is actually pretty cool as you can see a village of taxonomy sites being constructed. Some are fairly complete and interesting to browse through, like the Scratchpad about freeloader flies or the one about nanofossils.

Time will tell if Scratchpads will catch on with taxonomists. I certainly hope it does. It seems  quite a bit of thought has gone into balancing standardization with creativity, and making scratchpading easy and manageabe. The authors have chosen templates based on the Drupal content management system, which has a very gentle learning curve. I am waiting on my own scratchpad account, where I intend to… well, you’ll see it on this blog eventually, I hope.

So if you are involved in taxonomy in some way, consider getting your own Scratchpad. This could make for a great student project, by the way.  If you know someone who is a taxonomist, pass this information on to them.

I just spent the better part of a Saturday playing with Foldit. Foldit is an ongoing experiment in finding protein structures by harnessing the power of the mob – or gamers, as is the case here. The player is presented with a backbone & sidechain configuration, with the secondary structures mostly pre-determined. The problem is to get the protein to fold into the correct conformation. You can tweak the secondary structures, rubberband the beta strands together into sheets and rotate the sidechains. The residues are colored by hydrophobicity, so you know who should be facing where. The 23 tutorial cases walk you through the  simple yet powerful interface to folding the structures. You can rotate helices, rubberband strands, flip sheets, etc. The interface gives you feedback on sidechain clashes and voids in the structure among other things.  The examples also teach you the basic of protein structural considerations: maximize hydrogen bonds, hydrophobic side chains should be buried in the structure, strands should combine to sheets, and so on. When you feel you are ready, you can start solving the various folding puzzles presented online. You can work solo, or as part of a team. The “correct solution” is, of course, unknown. The best you can do is accumulate as many points as you can, which represent how stable is the conformation you are building. Foldit, like many other cool things in structural biology, is the product of David Baker’s lab at the University of Washington. Here is a video from the YouTube UWFoldit channel showing the coolness of it all. If you decide to get your fold on, make sure you can make the time. It’s flippin’ addictive.

The source file associated with this post can be downloaded here.

The last time I talked about how to read a GOA gene_associations file into a Python dictionary data structure.  Our goal was to find all genes that are annotated as hydrolases in the GOA gene_associations file. The tricky part is, most enzymes are not annotated with a general term such as “hydrolase_activity”, but with a more specific term, such “phosphatase_activity” or “tyrosine_phosphatase_activity”. “phosphatase_activity” can be viewed as a “private case” or in ontological terms a descendant of  hydrolase_activity. “tyrosine_phosphatase_activity” is even more specific, being a descendant of “phosphatase_activity” and, by transitivity of “hydrolase_activity”. The bottom line is,  we need to know all the descendants of  our query term to get the properly annotated genes form the GOA gene_associations file.

To do that, we would need to access the GO database. You can import it from the Gene Ontology FTP site as a MySQL dump, and install it locally. This would probably make your query faster. However, you would need to worry about maintaining and regularly updating your database from GO. Instead, it might make sense to connect to one of the mirrors of the GO MySQL database, and query it remotely. This is what I will do here.

To do that, I will create a class that does two things: (1) opens (and closes) a MySQL connection to a GO mirror site and (2) reads the GOA gene_associations file. We saw how to read the file in the previous part. What I am doing here is simply rolling that code into a method (a function called by a class).

import UserDict
import MySQLdb

class GOA_ga(UserDict.UserDict):
    def __init__(self):
        self.data = {}

    def godb_open(self):
        self.dbcon = MySQLdb.connect(host="mysql.ebi.ac.uk", user="go_select",
                   passwd="amigo", db="go_latest", port=4085)
        self.dbcursor = self.dbcon.cursor()

    def godb_close(self):
        self.dbcursor.close()

    def read_gene_assoc(self,inpath):
        for inline in file(inpath):
            db, db_object_id, db_object_symbol, qualifier, go_id, \
            db_reference, evidence, withit, aspect, \
            db_object_name, synonym, db_object_type, \
            taxon_id, date, assigned_by = inline.strip().split('\t')
            key = db+":"+db_object_id
            self.data.setdefault(key,[]).append({'db':db,'db_object_id':db_object_id,
                                            'db_object_symbol': db_object_symbol,
                                            'qualifier': qualifier,
                                            'go_id': go_id,
                                            'db_reference': db_reference,
                                            'evidence': evidence,
                                            'with': withit,
                                            'aspect': aspect,
                                            'db_object_name': db_object_name,
                                            'synonym': synonym,
                                            'db_object_type': db_object_type,
                                            'taxon_id': taxon_id,
                                            'date': date,
                                            'assigned_by': assigned_by})

What we have here is a class with four methods. The class itself is an instantiation of the Python UserDict class. This means it is primarily a dictionary, and can be used as such. But there are other methods too.

Lines 5-6: The __init__ method describes what happens when the class is instantiated, and is actually superfluous here. I just like to put it as good programming practice. Instantiating GOA_ga created an empty dictionary. That’s it.

Lines 8-11: The second method will open a MySQL connection to a remote server. dbcon and dbcursor will be used to access the database from Python, as we shall see soon. I import the Python MySQLdb module which contains functions to access MySQL through Python. You can get it from SourceForge, and it is also in the repositories of most Linux distributions. To install in Ubuntu Linux simply type at the command-line:

sudo aptitude install python-mysqldb

Lines 13-14: The go_dbclose method simply closes the database connection. It is good practice to do so after all queries are done.

Lines 16-36: this is basically the code to read the GOA gene_associations file that we saw in part 1.
The file go_associations.py contains this code (and more below). You can download it, import it, and run the following on your GOA gene association file of choice, using the Python shell. I chose arabidopsis in this example.

>>> import go_associations as ga
>>> my_goa = ga.GOA_ga() # create a class instantiation
>>> my_goa.read_gene_assoc("gene_association.goa_arabidopsis")  # read all GOA

So what we have is a class that lets us read GOA gene_associations file, and also query the GOA MySQL database.

Remember the original goal: get all the records associated with a certain function, or its children. The MySQL code for that is:

SELECT
  rchild.acc
FROM
  term AS rchild, term AS ancestor, graph_path
WHERE
  graph_path.term2_id = rchild.id and
  graph_path.term1_id = ancestor.id and
  ancestor.acc = 'GO:0016787'

This mysql query will return a list of all the GO accession codes that are the descendants of the go accession code “GO:0016187″. How do you know which GO accession code belongs to which term? You can use the keyword search feature in AmiGO. By the way, I pulled this query from the GO Example Queries page on the GO site. There are many useful queries there!

So how do we insert this into Python code? Here it is:

def find_go_descendants(goa_ga, parent_go_id):
   find_descendants_go = """
   SELECT
     rchild.acc
   FROM
     term AS rchild, term AS ancestor, graph_path
   WHERE
     graph_path.term2_id = rchild.id and
     graph_path.term1_id = ancestor.id and
     ancestor.acc = '%s'
   """
   go_children = {}
   found_genes = {}
   # query GO database for all the child terms of the parent_go_id
   goa_ga.dbcursor.execute(find_descendants_go % parent_go_id)
   rpl = goa_ga.dbcursor.fetchall()
   # Flatten the list. Put it in dictionary keys for faster searching.
   for i in rpl:
      go_children[i[0]] = None
   for gene_id in goa_ga.keys():
      for goa_rec in goa_ga[gene_id]:
         if goa_rec['go_id'] in go_children.keys():
            found_genes.setdefault(gene_id,[]).append(goa_rec)
   return found_genes

Lines 2-11 are the templated MySQL query to get all the descendants of the provided GO access number.
Line 12: go_children is a dictionary whose keys are all the GO accession numbers that are the descendants of our query. The values are null. It is faster to search a dictionary keys than a list.
Line 13: found_genes is a GOA_ga dictionary which will contain the found go_associations entries. The keys are the UniProt IDs, and the values are GOA_ga dictionary entries
Lines 15-16: Now we start working. First, we query the GO MySQL database for all the descendant GO terms of go_id. The returned value is a tuple of tuples which looks something like this:

(('GO:0012345',), ('GO:0034567',),('GO:003617',),....)

That is simply how MySQLdb’s fetchall method works: each MySQL record fetched is a tuple. Since we are fetching only one field with our query, each tuple element is  We now need to flatten this tuple — remove the superfluous parentheses. Lines 18-19 do that. They also put the descendant GO terms a dictionary keys structure, rather than a list structure. A dictionary keys structure is faster to search as we do repeatedly in line 23.

Line 20 is the outer loop of the search. We loop over all UniProt gene IDs in the goa_associations file we read into our goa_ga data structure.
Line 21 is the inner loop. Each Uniprot ID may appear in more than one record, as the protein it represents may be associated with more than one GO term.

Lines 22-23: if the associated GO term is in our list of descendant, add that GOA record to our found_genes dictionary, and key it with the proper UniProt ID.

Line 24: return the dictionary of all GOA records that are annotated as the descendants of the query term, or as the query term itself.

Putting it all together:

>>> import go_associations as ga
>>> my_goa = ga.GOA_ga() # create a class instantiation
>>> my_goa.read_gene_assoc("gene_association.goa_arabidopsis")  # read all GOA
>>> my_goa.godb_open()     # open a connection to the GO database
>>> all_hydrolases  = ga.find_go_descendants(my_goa,"GO:0016787") # find the descendants of "hydrolase_activity"

all_hydrolases now contains the GOA records of all the arabidopsis proteins annotated as hydrolases. How many are there?

>>> len(all_hydrolases)
3617
>>>

(This may vary when the latest version of this file changes. GOA is regularly updated). Some genes may be annotated as hydrolases more than once. This might either be an annotation error, as the is annotated twice with the same function, or the protein might be promiscuous with more than one function, or it may have two domains, each with a specific function. To find out all the hydrolase annotations:

>>> n=0
>>> for i in all_hydrolases:
...      n += len(all_hydrolases[i])
>>> print n
5753
>>>

How many of those annotations have a good evidence code?

>>> n=0
>>> for i in all_hydrolases:
...     for j in all_hydrolases[i]:
...             if  j['evidence'] in ('EXP','IDA','IPI','IMP','IGI','IEP'): n += 1
...
>>> print n
416
>>>

416 out of 5763: just over 7%. Keep that in mind next time you delve into an annotation database. (Exercise: find out how many enzymes have been annotated with good evidence in arabidopsis, or in another GOA-annotated organism).

Have fun GOing, and getting there!

Speaking of sampling bacteria, this ties in well with the previous post about GEBA. And by “well” I mean “in an alternate-universe/ altered-consciousness manner”.

The voices in the song are sampled from this KFC employee training tape. The video won a prize in machinima.com. So if you like World of Warcraft, bacteria, KFC, sampled music, or any combination of the above, you’re gonna love this.