Howto

I’ve just finished writing up an example for the Biopython Tutorial of the new GenomeDiagram functionality added in Biopython 1.59. You can now control the start and end points of individual tracks, and you can add cross-links between regions of different tracks, as shown here:

This example attempts a simplified reproduction of Figure 6 in Proux et al. (2002), and shows three related phage genomes one above the other. Different classes of genes have been given different colors, while the strength of the red shaded cross-links indicates the percentage identity of the linked genes. Note there are some minor differences in the GenBank annotation we’ve used and the genes shown in the original figure.

One of the new things coming in Biopython 1.59 is improved chromosome diagrams, something you may have seen via Twitter. I’ve just been updating the Biopython Tutorial (current version here, PDF) to include an example drawing this:

Here’s a PDF version too. This example just parses the Arabidopsis thaliana GenBank files to get the chromosome lengths and the tRNA gene placements. There are so many tRNA on the forward strand of Chr I that their labels are forced to overlap. Here the figure just uses a different color for each chromosome, but you can color each feature individually.

In another quirk to the FASTQ story, recent Illumina FASTQ files don’t actually use the full range of PHRED scores - and a score of 2 has a special meaning, The Read Segment Quality Control Indicator (RSQCI, encoded as ‘B’).

Hats off to Dr Torsten Seemann for raising awareness of this issue in his post on the seqanswers.com forum, referring to a presentation by Tobias Mann of Illumina which says:

The Read Segment Quality Control Indicator:

[Read More]

This is another blog post to highlight one of the neat tricks you’ll be able to do with Biopython 1.54 (which you can help test with the Biopython 1.54 beta release).

It is often useful to be able to extract a few records from a larger sequence file - for example, some sequences of interest from a full UniProt or GenBank dump. One obvious way to try to do this is to parse the file into an object representation (i.e. SeqRecord objects using Bio.SeqIO.parse(...)), filter to pick out the entries you want, and then write them back to disk (using Bio.SeqIO.write(...)). However, for complex file formats like GenBank this can be lossy ( Bio.SeqIO does not support a 100% identical round trip), and Biopython don’t currently support writing out the SwissProt plain text format used by UniProt. So, that approach won’t work.

Biopython 1.51 onward includes support for Sanger, Solexa and Illumina 1.3+ FASTQ files in Bio.SeqIO, which allows a lot of neat tricks very concisely. For example, the tutorial ( PDF) has examples finding and removing primer or adaptor sequences.

However, because the Bio.SeqIO interface revolves around SeqRecord objects there is often a speed penalty. For example for FASTQ files, the quality string gets turned into a list of integers on parsing, and then re-encoded back to ASCII on writing.

As per Peter’s recent post we are using this space to show of a couple of the new features in Biopython 1.52 before it is released. In this post we’ll look at the new convert() function that both Bio.SeqIO and Bio.AlignIO will get in Biopython 1.52.

No one has ever complained that bioinformatics just doesn’t have enough file formats - you probably frequently find yourself converting sequence files to suit particular applications with Bio.SeqIO. At the moment this is usually a two step process, something like this:

The forthcoming release of Biopython 1.52 will include a couple of nice improvements to the Bio.SeqIO module, and here we’re going to introduce the new index function. This will of course be covered in the Biopython Tutorial & Cookbook ( PDF) once this code is released.

Suppose you have a large sequence file with many many individual sequences in it. This could be next generation sequence data for example, maybe a FASTQ, FASTA or QUAL file. Or, it might be a big annotation rich file, such as the whole of UniProt, or a chunk of GenBank.

Constructing complicated NCBI Entrez searches can be tricky, but it turns out one of the Entrez Programming Utilities called Entrez EInfo can help.

For example, suppose you want to search for mitochondrial genomes from a given taxa - either just in the Entrez web interface, for use in a script with EFetch.

I knew from past experience about using name[ORGN] in Entrez to search for an organism name - but how would you specify just mitochondria? I actually worked this out from the NCBI help and exploring the Entrez website’s advanced search - but it took a while.

Hi all,

You may have noticed we’re trying out using the wiki for Biopython cookbook entries. It’s a new idea so at the moment there are only a few ‘recipes’ on offer. If you have some tricks you find yourself using time and again to solve a problem why not share them? Similarly, if you find yourself coming up against a problem you can’t seem to solve easily with Biopython’s tools send a message to one of the mailing lists proposing it as a cookbook example and someone just might solve it for you!