Viewing online figures in an Nucleic Acids Research article…an exercise in frustration

I found a Nucleic Acids Research article that I wanted to read online. The article in question had — like most articles — some figures. Here's how my web browser displays one of the figures:

View fullsize

I've blurred out the surrounding text so as not to risk any copyright issues (and also to let you focus on just the figure). If your eyesight is like mine, you may feel that the three subpanels of this figure are too small to be of much use.

So I clicked the image…

View fullsize

The white rectangle enclosing the figure increases in size from about 250 x 130 pixels to about 460 x 210. If I squint, I can just about make out some of the figure text, but it remains far from clear.

So I clicked the image…

View fullsize

Now I see exactly the same image as before, but it's in a new webpage which has a wider figure legend. The new legend is acually a little bit harder to read than previously as there are more words per line. I'm really not sure what the point of this page is.

So I clicked the image…

View fullsize

Hooray! Now I can finally see the figure at a decent size where all of the figure text is legible. It only took me three clicks to get there. To make sense of the figure I turn to the legend, only to find that there is no legend!

So you can have your figure at a reasonable size without the legend, or you can have the legend but only with a small version of the figure. It is obviously beyond the journals ability to give you both.

Great post by Lex Nederbragt about why we need graph-based representations of genome sequences →

In a recent blog post, Lex Nederbragt explains why we all need to be moving to graph-based representations of genome sequences (and getting away from linear representations).

In this post I will provide some background, explain the reasons for moving towards graph-based representations, and indicate some challenges associated with this development.

After listing the many challenges involved in moving towards a graph-based future, he refers to the fact that current efforts have not been widely adopted:

Two file formats to represent the graph been developed: FASTG and GFA. FASTG has limited uptake, only two assembly programs (ALLPATHS_LG and SPAdes) will output in that format. GFA parsing is currently only experimentally in the ABYSS assembler, and [the VG program] is able to output it.

The lack of a widely recognized and supported standard for representing variation inherent in a genome sequence is, in my opinion, a major barrier to moving forward. Almost all bioinformatics software that works with genome sequence data expects a single sequence with no variation. It will require a whole new generation of tools to work with a variant-based format, but tool developers will be reluctant to write new code if there is no clear agreement on what is the new de facto file format.

Duplicate names in bioinformatics: #10 MAGMA

I've written about this problem enough times now, that I thought I'd promote it to its own series. Check out the previous posts that deal with this subject.

Compare and contrast:

• 2002: MAGMA: Multiobjective Analyzer for Genetic Marker Acquisition
• 2007: MAGMA: analysis of two-channel microarrays made easy
• 2007: MAGMA: An Algorithm for Mining Multi-level Patterns in Genomic Data
• 2008: MAGMA (molecular animation, graphics and modeling application framework)
• 2015: MAGMA: Generalized Gene-Set Analysis of GWAS Data

Filtering a SAM file generated by TopHat to find uniquely mapped, concordant read pairs: AWK vs SAMtools

Software suites like SAMtools (or should that be [SAMsam]tools?) offer a powerful way to slice and dice files in the SAM, BAM, and CRAM file formats. But sometimes other approaches work just as well.

If you have aligned paired RNA-Seq read data to a genome or transcriptome using TopHat you may be interested in filtering the resulting SAM/BAM file to keep reads that are:

• a) uniquely aligned (only match one place in target sequence)
• b) mapped as concordant read pairs (both pairs map to same sequence in correct orientation with suitable distance between them)

TopHat has a --no-discordant command-line option which only report read alignments if both reads in a pair can be mapped, but the name of option is somewhat misleading, as you still end up with discordantly mapped read pairs in the final output (always good to check what difference a command-line option actually makes to your output!.

So if you have a TopHat SAM file, and you wanted to filter it to only keep uniqueley mapped, concordant read pairs, you could use two of the options that the samtools view command provides:

• -q 50 — This filters on the MAPQ field (5th column of SAM file). TopHat uses a value of 50 in this field to denote unique mappings (this important piece of information is not in the TopHat manual).
• -f 0x2 — This option filters on the bitwise FLAG score (column 2 of SAM file), and will extract only those mappings where the 2nd bit is set. From the SAM documentation, this is described as each segment properly aligned according to the aligner. In practice this means looking for values that are either 83, 89, 147, or 163 (see this helpful Biobeat blog post for more information about this).

So if you have a SAM file, in essence you just need to filter that file based on matching certain numbers in two different columns. This is something that the Unix AWK tool excels at, and unlike SAMtools, AWK is installed on just about every Unix/Linux system by default. So do both tools give ou the same result? Only one way to find out:

Using SAMtools

The 'unfiltered.sam' file is the result of a TopHat run that used the --no-discordant and --no-mixed options. The SAM file contained 34,340,754 lines of mapping information:

time samtools  view -q 50 -f 0x2 unfiltered.sam > filtered_by_samtools.sam

real    1m57.068s
user    1m18.705s
sys     0m13.712s

Using AWK

time awk '$5 == 50 && ($2 == 163 || $2 == 147 || $2 == 83 || $2 == 99) {print}' unfiltered.sam  > filtered_by_awk.sam

real    1m31.734s
user    0m46.855s
sys     0m15.775s

Does it make a difference?

wc -l filtered_by_*.sam
31710476 filtered_by_awk.sam
31710476 filtered_by_samtools.sam

diff filtered_by_samtools.sam filtered_by_awk.sam

No difference in the final output files, with AWK running quite a bit quicker than SAMtools. In this situation, filtering throws away about 8% of the mapped reads.

Correcting needless jargon in bioinformatics literature: part 1

A new series in which I attempt to streamline the existing corpus of bioinformatics literature by removing some of the largely pointless jargon that can occur in the titles and abstracts of papers.

Compare my before and after attempts at reworking this 2010 paper published in Genome Biology:

Before:

After:

More duplicate names for bioinformatics software: a tale of two HIPPIES

Image from Flickr user ochre_jelly

Thanks to Sara Gosline (@sargoshoe) for bringing this to my attention. Compare and contrast the following:

• HIPPIE: Integrating Protein Interaction Networks with Experiment Based Quality Scores
• HIPPIE: a high-throughput identification pipeline for promoter interacting enhancer elements

The former tool, published in 2012 in PLOS ONE, takes its name from 'Human Integrated Protein-Protein Interaction rEference' (it was doing so well until it reached the last letter). The latter tool ('High-throughput Identification Pipeline for Promoter Interacting Enhancer elements') was published in 2014 in the journal Bioinformatics.

Leaving aside the issue of whether these names are worthy of a JABBA award, the issue here is that we have yet another duplicate set of software names for two different bioinformatics tools. The authors of the 2nd paper could, and should, have checked for 'prior art'.

If you are planning to develop a new bioinformatics tool and have thought of a possible name, please take the time to do the following:

1. Visit http://google.com (or your preferred web search engine of choice)
2. In the search box type the proposed name of your tool followed by a space
3. Then add the word 'bioinformatics'
4. Click search
5. That's it

Inconsistent bioinformatics branding: SAMtools vs Samtools vs samtools

View fullsize

The popular Sequence Alignment Map format, SAM, has given rise to an equally popular toolkit for working with SAM files (and BAM, CRAM too). But what is the name of this tool?

SAMtools?

If we read the official publication, then we see this software described as 'SAMtools' (also described by Wikipedia in this manner).

Samtools?

Head to the official website and we see consistent references to 'Samtools'.

samtools?

Head to the official GitHub repository and we see consistent references to 'samtools'.

This is not exactly a problem that is halting the important work of bioinformaticians around the world, but I find it surprising that all of these names are in use by the people that developed the software. Unix-based software is typically — but not always — implemented as a set of lower-case commands and this can add one level of confusion when comparing a tool's name to the actual commands that are run ('samtools' is what you type at the terminal). However, you can still be consistent in your documentation!

How do people choose a single isoform of a gene to use for bioinformatics analyses?

Update 2015-09-29: in addition to the comments at the end of the post below, also see the follow up post that I wrote which offers some more suggestions including the APPRIS database/webserver which looks very useful.

This post is somewhat of a follow-up to something that I wrote earlier this week. In bioinformatics, we often want to analyze all genes from an organism (or from multiple organisms). In many well-annotated genome databases, there is often a choice of isoforms available for each protein-coding gene, and the number of isoforms only ever seems to increase.

For example, in the latest set of human gene annotations (Ensembl 78), there are 406 protein-coding genes that have more than 25 transcripts. At one extreme, the human GPR56 gene has 77 transcripts, 61 of which are annotated as protein-coding! The length of these 61 putative protein products ranges from just 6 amino acids (!) all the way up to 693.

In Caenorhabditis elegans, sequence identifiers for genes were historically based on appending numbers to the identifier of the BAC/YAC/Cosmid clone containing that gene. E.g. B0348.1 would represent the first predicted gene on the B0348 clone, B0348.2 the second gene…and so on. When splice variants were discovered, curators appended letters for each isoform. E.g. B0348.2a and B0348.2b represent the two alternative isoforms of this gene. In the latest WS248 release of WormBase, one gene (egl-8) has 25 isoforms (all the way up to B0348.4y). I wonder what WormBase will do when a 27th isoform is discovered?

So how does one attempt to choose a single variant for use in a bioinformatics pipeline, and is this something that we should even be attempting? Historically, people have often opted for a quick-and-easy approach in order to get around this problem. Some examples from papers indexed by Google Scholar:

"In cases of alternative splicing, we chose the longest protein to represent a gene"

"In cases of multiple transcript isoforms, we chose the isoform with the longest CDS supported by transcript and protein homology in other mammalian species"

"Because of the redundancy of protein sequences, we chose only the longest isoform for every entry"

"In cases where a gene possesses more than one reference sequence, we chose the longest"

"When multiple protein entries are found for the same EntrezGene identifier, choose the longest sequence isoform"

This methodology is obviously not without problems (as others have reported on). So I'm genuinely curious as to what people do in order to choose a 'representative' isoform (whatever that means). The problem is further complicated when the reality might be that some genes consistently use different isoforms in different tissues or at different developmental time points.

Please comment below if you think you have found a good solution to this problem!

24 carat JABBA awards

Here is a new paper published in the journal PLOSBuzzFeed…sorry, I mean PLOS Computational Biology:

• A DIseAse MOdule Detection (DIAMOnD) Algorithm Derived from a Systematic Analysis of Connectivity Patterns of Disease Proteins in the Human Interactome

It's a good job that they mention the name of the algorithm ninety-one times in the paper, otherwise you might forget just how bogus the name is. At least DIAMOnD has that lower-case 'n' which means that no-one will confuse it with:

• Fast and sensitive protein alignment using DIAMOND

This second DIAMOND paper dates all the way back to November 2014. Where does this DIAMOND get its name?

Double Index AlignMent Of Next-generation sequencing Data

This DIAMOND gets a bonus point for having a website link in the paper which doesn't seem to work.

So DIAMOnD and DIAMOND are both the latest recipients of JABBA awards for giving us Just Another Bogus Bioinformatics Acronym.