Rationale

Visualization is an important component of genomic analysis, primarily because itfacilitates exploration and discovery, by revealing patterns of variation andrelationships between experimental data sets and annotations. Data on the genome fallinto two classes: annotations, such as gene models, and experimental measurements, suchas alignments of high-throughput sequencing data. The unique and unifying trait of allgenomic data is that they occupy ranges on the genome. Associated with the ranges isusually multivariate meta-information both at the feature level, such as a score orfunctional annotation, and at the sample level, such as gender, treatment, cancer orcell type. These data ranges can range in scale from hundreds to billions of datapoints, and the features are dispersed along genomes that might be many gigabases inlength. Visualization tools need to slice and dice and summarize the data in differentways to expose its different aspects and to focus on different resolutions, from asensible overview of the whole genome, to detailed information on a per base scale. Tohelp focus attention on interesting features, statistical summaries need to be viewed inconjunction with displays of raw data and annotations.

Various visualization tools have been developed, most of which are implemented in theform of a genome browser. Data are typically plotted along with annotations with genomiccoordinates on the horizontal axis with other information laid out in different panelscalled tracks. Examples of genome browsers include the desktop-based browsers IntegratedGenome Browser [1, 2] and Integrative Genomics Viewer [3, 4]. There are also web-based genome browsers, including Ensembl [5], UCSC Genome Browser [6], and GBrowse [7], and several new web-based browsers, like Dalliance, which rely ontechnologies like HTML5 and Scalable Vector Graphics [8], or Adobe Flash, like DNAnexus [9]. Other software, like Circos [10], provide specialist functionality. R also has some new tools for visualizinggenomic data, GenomeGraphs [11] and Gviz [12]. They all have advantages for different purposes: some are fast, while othershave easier user interfaces. Some are interactive, offer cross-platform support orsupport more file formats.

Data graphics benefit from being embedded in a statistical analysis environment, whichallows the integration of visualizations with analysis workflows. This integration ismade cohesive through the sharing of common data models [13]. In addition, recent work on a grammar of data graphics [14, 15] could be extended for biological data. The grammar of graphics is based onmodular components that when combined in different ways will produce different graphics.This enables the user to construct a combinatoric number of plots, including those thatwere not preconceived by the implementation of the grammar. Most existing tools lackthese capabilities.

A new R package, ggbio, has been developed and is available on Bioconductor [16]. The package provides the tools to create both typical and non-typicalbiological plots for genomic data, generated from core Bioconductor data structures byeither the high-level autoplot function, or the combination of low-level components ofthe grammar of graphics. Sharing data structures with the rest of Bioconductor enablesdirect integration with Bioconductor workflows.

Basic usage

In ggbio, most of the functionality is available through a single command, autoplot,which recognizes the data structure and makes a best guess of the appropriate plot.Additional file 1 Table S1 lists the type of plot produced foreach type of data structure. The plot in Figure 1 was renderedwith the autoplot function in ggbio, using the following code:

Figure 1
figure 1

Gene structure. Example of exons of gene SSX4 and SSX4B isoforms, annotated to illustrate the grammar of graphics extensions used.The filled rectangle represents exons and the chevron represents introns. They aregrouped by transcript ID and the y axis shows the stepping levels, which stackstranscripts to avoid overplotting. Color is mapped from strand direction.

Full size image
autoplot ( grl , aes ( color = strand ) )

The data, grl, is a GRangesList object, a Bioconductor data structure for representingcompound ranges, including a set of transcript structures. The autoplot functionrecognizes the GRangesList object and draws the intervals in the typical fashion for agene model or alignment. The y axis is generated by a layout algorithm to ensure thatthe transcripts are not overplotted. The x axis is automatically set as the genomiccoordinates. The call to aes maps the strand variable to the color aesthetic. Users canspecify labels, titles, layout, and so on by passing additional arguments to autoplot.Compared to the more general qplot API of ggplot2, autoplot facilitates the creation ofspecialized biological graphics and reacts to the specific class of object passed to it.Each type of object has a specific set of relevant graphical parameters, and furthercustomization is possible through the low-level API, introduced later.

Plotting tracks

In many genome visualizations, different datasets are typically plotted separately andthen stacked on top of the same x axis, the genomic coordinates. These plots are oftencalled tracks, because they are usually much wider than they are tall and run parallelto each other. Each track might contain a heatmap, or a histogram, or a density plot,for example. The data displayed in each track is related to the data displayed in theother tracks through the shared genomic axis. The goal is to be able to observedifferent patterns in these snapshots for regions of the genome. When displaying arelatively large region, the tracks might show a summary of the data, whereas moredetails are depicted for smaller regions.

The ggbio package provides a function called tracks, which stacks plots in a specifiedorder and creates a Tracks object. The object allows users to zoom and shift the plot,as well as modify parameters like the height and theme of individual tracks. In thefollowing example, p1, p2 and p3 are plot objects in the session:

tracks ( p 1 , p 2 , p 3 , heights = c ( 4 , 1 , 1 ) )

Figure 2 illustrates the creation of tracks. One tumor/normal pairof RNA-seq samples is shown with the goal of interpreting splicing changes. A key aspectof the plot is how the read alignment coverage and junction counts are composited andthen juxtaposed with the data from the other sample, as well as the annotated transcriptmodels. The viewer is then able to relate the changes in coverage to the correspondingtranscript structures, via the common x axis. At the top of the figure is an ideogramoverview of the chromosome, using colors corresponding to the Giemsa stain.

Figure 2
figure 2

Splicing summary with coordinate truncate gaps. An example of a plot madefrom multiple tracks. At the top, the relevant chromosome is drawn with thesubregion of interest marked in red. The middle track shows the slicing summaryplots for the gene ALDOA for normal and tumor samples. Splicing is shownas arches and size is used to represent junction counts and color representsnovelty: blue indicates known splicing events against the model and red indicatesnovel splicing events. The height of arches is proportional to the distancebetween the two ends of the arches, or the distance between the junction reads.Coverage is shown by position to address supporting evidence in the raw data. Thesplicing summary plots are aligned with a view of the gene structure in the bottomtrack. The thicker rectangle represents the Consensus Coding DNA Sequence (CCDS)transcripts. The plots in both tracks are made with truncate gaps coordinatetransformation. The space dedicated to introns has been significantly reduced, andthe exonic regions are shown in detail, even though the entire gene region is inview.

Full size image

Genomic overviews

The purpose of an overview plot is to give a grand view of the entire genome. Bydefinition, this means that the resolution will be poor and that only large featureswill be visible. An overview may reveal large features that might be missed if onefocused too narrowly. Different methods for mapping the genomic axis to the screen havebeen applied to address the space issues, and also to ease the drawing of connectionsbetween regions.

Grand linear view

Figure 3 shows an example of a grand linear view, which laysout the entire genome along a single linear axis. The plot shown in the figure is aspecial case of the grand linear plot called a Manhattan plot, due to its resemblanceto the Manhattan skyline. It emphasizes extreme events, which show up as high-valuedoutliers. This view has been used for many genome-wide association study reports [17]. The data here come from a genome-wide association study on Angus cattle,and the data are faceted by three different scoring and classification methods [18]. The horizontal axis shows the global genomic coordinates, and thevertical axis is generally mapped to some quantity of interest, in this case, thegenetic variance.

Figure 3
figure 3

Manhattan plot. Grand linear view applied to a Manhattan plot as part ofa genome-wide association study in Angus cattle. The y axis shows geneticvariance, calculated by sliding windows of five consecutive SNPs for theinfectious bovine keratoconjunctivitis (IBK; a type of pinkeye) score. The xaxis is the genomic coordinates with all the chromosomes side-by-side. Thehorizontal striping of color helps to indicate the end of one chromosome andbeginning of another. The plot is faceted by three different analysis methods.There is one extreme variance in the middle facet, in the region of chromosome23. There are also a few large values in other regions. According to theresults from the paper, three of these regions, 2, 13, and 23 are found to bepotentially indicative of a quantitative trait locus associated with IBK.

Full size image

The plot uses linear layout and employs the genome coordinate transformation, whichtransforms the chromosomal coordinates into global genomic coordinates as if all ofthe chromosomes were concatenated together. The transformation supports bothproportional and uniform scaling of chromosomes, so that the plot area consumed by achromosome is either proportional to its length or the same as the other chromosomes.It is also possible to add a buffer or break between chromosomes.

Karyogram overview

Figure 4 shows a (single copy) karyogram overview plot, withthe color indicating RNA-editing locations in human [19]. The karyogram layout represents chromosomes as rectangles and stacks themvertically or in a grid layout. Genomic position is relative to the chromosome,starting from the first position at the left for each rectangle. Associatedinformation is overlaid on or over the box. Applications like Genome Graphs in theUCSC genome browser [6] provide a similar plot. The advantage of this layout over the grand linearview is more efficient use of horizontal space, and hence finer resolution detail onthe positions of the features. The trade-off is that the second variable has lessspace - instead of a full vertical axis the information needs to be fit into eachrectangle. Thus, we obtain genomic position resolution at the cost of less data layerresolution. It is common to see this type of plot used for SNP density, varyinglevels of identity-by-descent (IBD) [20] and length of linkage disequilibrium spans [21].

Figure 4
figure 4

Stacked karyogram overview. Karyogram plot shows a subset of humanRNA-editing sites, and they are color coded for different regions as follows:red indicates exons, green indicates introns and blue indicates exons/intronsstatus is unknown.

Full size image

Circular overview

The primary purpose of the circular view is to show links between genomic regions.This is generally infeasible with the linear or karyogram layouts. In a circularlayout, features are organized into concentric rings. Figure 5illustrates the circular overview on the data from a gene fusion study conducted byBass and colleagues [22], who sequenced the genomes of nine individuals with colorectal cancer andidentified an average of 75 somatic rearrangements per tumor sample. This circularview shows only a single sample (colorectal tumor sample CRC-1), the structuralrearrangements are shown as links with intrachromosomal events in green andinterchromosomal translocations in orange. An ideogram of the autosomes is shown inthe outer ring, with somatic mutation and score tracks in the plot. There are manysoftware packages that provide a circular overview plot, including Circos [10], CGView [23] and DNAplotter [24].

Figure 5
figure 5

Single sample circular view. DNA structural rearrangements and somaticmutation in a single colorectal tumor sample (CRC-1). The outer ring shows theideogram of the human autosomes, labeled with chromosome numbers and scales.The segments represent the missense somatic mutations. The point tracks showscore and support for rearrangement. The size of the points indicates thenumber of supporting read pairs in the tumor and the y value indicates thescore for each rearrangement. The links represent the rearrangements, whereintrachromosomal events are colored green and interchromosomal events arecolored orange.

Full size image

Specialized plots

There are some typical types of plots used to examine specific biological questions.This section describes how ggbio builds two of these: a mismatch summary and anedge-linked interval plot.

Mismatch summary

Mismatch summary is one typical way to visualize alignments from sequencing data,especially in the context of variant calling. Other genome browsers, such asIntegrative Genomics Viewer [4], Savant [25] and Artemis [26] render similar plots from BAM and variant call format (VCF) files. Figure6 shows two different summaries of mismatches from a set ofRNA-seq read alignments. The top plot shows one DNA-seq sample from the first phaseof the 1000 Genomes Project [27], represented as a stacked barchart. It provides a detailed view of thecoverage, where the counts of bases that match the reference are indicated by graybars, and the counts of non-reference bases are indicated by a different color thatis specific to the base (A, C, G or T).

Figure 6
figure 6

Mismatch summary. An example of a mismatch summary plot, with associatedvariant calls. The top track shows a barchart of reference counts in gray andmismatched counts colored by the nucleotide. The middle track shows SNPs asletters, color coded also by nucleotide. There is one mismatch, 'T', that isdifferent for all of the reads from the 'A' in the reference genome (bottomletter plot).

Full size image

Edge-linked interval to data views

Interval data, like genes, regulatory sites, read alignments, and so on, aredifferent lengths. Differences in length can be distracting when looking atassociated numerical information. Thus, length is sometimes best ignored, and theinterval treated as an id or categorical variable. Figure 7shows an example. The top plot shows a profile display of expression levels for twosamples, GM12878 and K562, where the genomic position of the exons is treated as acategorical variable, forcing equal width in the plot. This allows us to see exonswhere the expression level is different, without being distracted by the relativeinterval size of the exons. We could also consider this display to be a parallelcoordinate plot [28, 29].

Figure 7
figure 7

Edge-linked interval to data view. Edge-linked interval to data view forthe expression of the exons of gene PDIA6. The top track shows theexpression level for each of the exons, and the color indicates the sample(GM12878 or K562). The second track shows the links between the even-spacedexpression track and the exons track, below. The package DEXseq, which producesa similar graphic, computes differential expression and significance, andsignificance is indicated by coloring the connecting lines red. The track atthe bottom shows the annotated transcripts.

Full size image

Biological extensions to the grammar of graphics

The base grammar

The work introduced in this paper builds upon the grammar of graphics conceived byWilkinson [14] and expanded upon by Wickham [15]. The grammar is composed of interchangeable components that are combinedaccording to a flexible set of rules to produce plots from a wide range of types.Table 1 explains the components of the grammar (data, geom,stat, scales, coord, facet), as utilized in ggplot2, and indicates how these are usedto create two graphics: Figure 8 and 1.

Table 1 Components of basic grammar of graphics
Full size table
Figure 8
figure 8

MA-plot. MA-plot for differential expression analysis in four RNA-seqsamples with two cell lines GM12878 and K562, annotated to illustrate the useof the grammar of graphics. Points is our geometric object, x axis indicatesthe normalized mean and the y axis indicates the log2 fold change.Aesthetics mapping took place between the groups and the color to use red toindicate the most significant differently expressed observation (gene). Thisplot uses Cartesian coordinates.

Full size image

Genomic data and abstractions

Data are the first component of the grammar, and data may be collected in differentways. Wilkinson makes a distinction between empirical data, abstract data andmetadata [14]. Empirical data are collected from observations of the real world, whileabstract data are defined by a formal mathematical model. Metadata are data aboutdata, which might be empirical, abstract or metadata themselves. We will use the termdata source to refer to concrete data in specific databases and file formats. This isroughly analogous to Wilkinson's empirical data.

The ggbio package attempts to automatically load files of specific formats intocommon Bioconductor data structures, using routines provided by Bioconductorpackages, according to Additional file 1 Table S2. Theloaded data are then considered by Wilkinson to be abstract, in that they are nolonger tied to a specific file format. Analogously, a data structure may be createdby any number of algorithms in R; all that matters is that every algorithm returns aresult of the same type. The type of data structure loaded from a file or returned byan algorithm depends on the intrinsic structure of the data. For example, BAM filesare loaded into a GappedAlignments, while FASTA and 2bit sequences result in aDNAStringSet. The ggbio package handles each type of data structure differently,according to Additional file 1 Table S1. In summary, thisabstraction mechanism allows ggbio to handle multiple file formats, withoutdiscarding any intrinsic properties that are critical for effective plotting.

Extension overview

Genomic data have some specific features that are different from those of moreconventional data types, and the basic grammar does not conveniently capture suchaspects. The grammar of graphics is extended by ggbio in several ways, which areillustrated in Figure 9 and described in Additional file1 Table S3. These extensions are specific to genomicdata, that is, genomic sequences and features, like genes, located on thosesequences.

Figure 9
figure 9

Diagram of the ggbio framework for processing sequence data. It startswith a mapping from different file types to different objects or data structurein R, using Bioconductor tools, followed by general and extended grammar ofgraphics mapping of data elements to graphical components. The final stagearranges the graphics in a designed layout to show annotation tracks ormultiple data sets. Orange boxes and dark brown arrows indicate the extensionsprovided by ggbio.

Full size image

Figure 1 illustrates how the components of the grammar arecombined to plot gene structures. A sample of the data is shown in Table 2. The data are passed to ggbio as a GRangesList object. Thechevron geom mimics the typical splice junction diagram found in textbooks and itdraws the introns in the example. The exons are drawn using the rectangle geom, andthe high-level alignment geom ensures that the introns and exons from the sametranscript are drawn connected, according to the tx_id column. The position on thechromosome is mapped to the horizontal axis and strand is mapped to color. Thevertical axis is mapped to a variable generated by the stepping statistic, whichavoids overplotting between the transcripts. We will explain these aspects in thefollowing sections.

Table 2 Example of GRanges object
Full size table

Biological geometries

A geom is responsible for translating data to a visual, geometric representationaccording to mappings between variables and aesthetic properties on the geom. Incomparison to regular data elements that might be mapped to the ggplot2 geoms ofpoints, lines and polygons, genomic data has the basic currency of a range. Rangesunderlie exons, introns and other features, and the genomic coordinate system formsthe reference frame for biological data. We have introduced or extended several geomsfor representing ranges and gaps between ranges. They are listed in Additional file1 Table S3. For example, the alignment geom delegates totwo other geoms for drawing the ranges and gaps. These default to rectangles andchevrons, respectively. Having specialized geoms for commonly encountered entities,like genes, relegates the tedious coding of primitives, and makes user code simplerand more maintainable.

Biological statistical transformations

A statistical transformation (stat) transforms or summarizes the data in a particularway. These statistics may be mapped to visual aesthetics in the same manner as theoriginal data. In this work we introduce several statistical transformationsspecifically for genomic data as shown in Additional file 1Table S3. For example, given a large number of read alignments, computation ofcoverage is useful, as shown in Figure 10. Thesetrans-formations were implemented with significant help from Bioconductor tools. Newstatistical transformations are readily incorporated.

Figure 10
figure 10

Coverage transformation. Statistical transformation, coverage andstepping, are used to summarize short reads data. Top: a set of (simulated)short reads, displayed using the stepping transformation, vertically, and thedefault geom 'rectangle'. Bottom: coverage is shown on the vertical axis, usingthe geom 'area'. This example applies the data model GRanges object.

Full size image

Biological coordinate transformations

Coordinate systems locate points in space, and we use coordinate transformations tomap from data coordinates to plot coordinates. The most common coordinate system instatistical graphics is cartesian. The transformation of data to cartesiancoordinates involves mapping points onto a plane specified by two perpendicular axes(x and y). Why would two plots transform the coordinates differently for the samedata? The first reason is to simplify, such as changing curvilinear graphics tolinear, and the second reason is to reshape a graphic so that the most importantinformation jumps out at the viewer or can be more accurately perceived [14].

Coordinate transformations are also important in genomic data visualization. Forinstance, features of interest are often small compared to the intervening gaps,especially in gene models. The exons are usually much smaller than the introns. Ifusers are generally interested in viewing exons and associated annotations, we couldsimply cut or shrink the intervening introns to use the plot space efficiently. Forexample, Figure 2 is able to show the entire gene region, withvirtually no loss in data resolution. In ggbio, we propose three sets of coordinatesystems, shown in Additional file 1 Table S3, which areuseful for genomic data.

Biological faceting

Almost all experimental outputs are associated with an experimental design and othermeta-data, for example, cancer types, gender and age. Faceting allows users to subsetthe data by a combination of factors and then lay out multiple plots in a grid, toexplore relationships between factors and other variables. The ggplot2 packagesupports various types of faceting by arbitrary factors. The ggbio package extendsthis notion to facet by a list of ranges of interest, for example, a list of generegions. There is always an implicit faceting by sequence (chromosome), because whenthe x axis is the chromosomal coordinate, it is not sensible to plot data fromdifferent chromosomes on the same plot. As an aside, generating a set of tracks mightresemble faceting, but it is easier to fit into the grammar framework if we think ofit as a post-processing step.

Biological layout

We have also extended the grammar of graphics with an additional component calledlayout, upon which the mapping from genomic coordinates to plot coordinates depends.The default layout simply maps the genomic coordinates to the x axis and facets bychromosome. The currently supported layouts are: linear (genomic coordinates mappedto the x axis), karyogram (each chromosome displayed separately, in an array), andcircular (like linear, except wrapped around in a circle). The high-level genomicoverview plots take advantage of these layout mechanisms.

Low-level grammar-oriented API

For custom use cases, ggbio provides a low-level API that maps more directly tocomponents of the grammar and thus expresses the plot more explicitly. Generallyspeaking, we strive to provide sensible, overridable defaults at the high-level entrypoints, such as autoplot, while still supporting customizability through the low-levelAPI.

All lower level functions have a special prefix to indicate their role in the grammar,like layout, geom, stat, coord, and theme. The objects returned by the low-level API maybe added together via the conventional + syntax. This facilitates the creation of newtypes of plots. A geom in ggplot2 may be extended to work with more biological datamodel, for example, geom rect will automatically figure out the boundary of rectangleswhen the data is a GRanges, as do geom bar, geom segment, and so on. As an example, thefollowing code produces the same plot as the code shown above, using the low-level APIinstead of autoplot:

ggplot ( ) + geom _ arrowrect ( unlist ( grl ) ) + geom _ chevron ( gaps ( unlist ( grl ) ) )

The reader will notice how the low-level code is more descriptive about the compositionof the plot. In this example, it says we start with an empty plot as created by ggplot.We then use geom arrowrect for exons and add a second layer for the gaps using geomchevron.

The low-level API may be used in conjunction with autoplot, via the + syntax. This makesit possible to save a plot as an object in a session and modify it in different wayswhile experimenting. For example, the following code applies a new theme to the existinggraphic object. The theme null function removes the background labels and legend.

p < - autoplot ( grl , aes ( color = strand ) ) p + theme _ null ( )

Materials and methods

The ggbio package is an extension for R, a free cross-platform programming environmentfor statistical analysis and graphics with more than 3, 000 contributed packages. Thepackage depends upon Bioconductor libraries for handling and processing data, includingthe implementation of the statistics in our extension of the grammar. The Bioconductorproject is a collaborative effort to develop software for computational biology andbioinformatics with high-quality packages and documentation [30]. The visualization methods in ggbio depend heavily on the package ggplot2 [15], which implements the grammar of graphics. The new geoms in ggbio areconstructed from primitives defined in ggplot2. We use ggplot2 as the foundation forggbio, due to its principled style, intelligent defaults and explicit orientationtowards the grammar of graphics model. The color schemes in ggbio were derived fromstandard palettes available in R [31–33].

The RNA-seq data used in this paper are from ENCODE [34]. Two cell lines, GM12878 (blood, normal, female) and K562 (blood, cancer,female), are selected, and there are two replicates for each sample. The data weremapped against hg19 using Spliced Transcript Alignment and Reconstruction (STAR) [35]. The Bioconductor packages Rsamtools [36] and GenomicRanges [37] were used to import the BAM files and count reads overlapping exons. Thepackage DEXSeq [38] was used to conduct the expression analysis and find the most differentlyexpressed exons. We used the rtracklayer package [39] to import BED format files and cast them into GRanges objects for ggbio. TheDNA-seq BAM files and VCF files used in Figure 6 were downloadedfrom the 1000 Genomes Project [27].

All figures, code and data links are available from the documentation section of theggbio website [40].

Discussion

We have demonstrated how ggbio supports both the convenient construction of typicalgenomic plots, while simultaneously supporting the invention of new types of plots fromlow-level building blocks. Use cases of ggbio range from generating reproducible,exploratory plots in the course of an analysis to the prototyping of new ways of lookingat these complex data. Lessons learned might be applied to the design of more complex,interactive systems. A new package, visnab, is being developed to make interactivegraphics for genomic data [41].

One such lesson learned is the importance of color choices, which are inconsistent inmany existing tools. Color is one of the primary visual clues in a data graphic andneeds to be handled with some intelligence. For example, the ggbio package builds onwell-specified color palettes used in ggplot2 and biovizBase, including one that isbased on the biologically inspired Giemsa stain colors, as shown at the top of Figure2.