Software for Computing and Annotating Genomic Ranges — TEI-tagged view

Michael LawrencepersName, Wolfgang HuberpersName, Hervé PagèspersName, Patrick AboyounpersName, Marc CarlsonpersName, Robert GentlemanpersName, Martin T. MorganpersName, Vincent J. CareypersName (2013date). “Software for Computing and Annotating Genomic Ranges.” PLoS Computational Biologytitle·j.

Object of study — software infrastructure for representing, computing on, annotating, and integrating genomic rangesterm in R

Research question — how to make range-based genomic analysis efficient, interoperable, reproducible, and statistically integrated across packages and data types

Contribution — the IRangesterm, GenomicRangesterm, and GenomicFeaturesterm infrastructure with scalable rangeterm structures, operations, annotationterm access, and experiment containers

Full text — tagged

Software for Computing and Annotating Genomic Rangesterm Abstract We describe Bioconductor infrastructureterm for representing and computing on annotated genomic rangesterm and integrating genomic data with the statistical computing features of R and its extensions. At the core of the infrastructureterm are three packages: IRangesterm , GenomicRangesterm , and GenomicFeaturesterm . These packages provide scalable data structuresterm·mth for representing annotated rangesterm on the genome, with special support for transcript structuresterm, read alignmentsterm and coverage vectorsterm. Computational facilities include efficient algorithms for overlap and nearest neighbor detectionterm, coverageterm calculation and other range operationsterm. This infrastructure directly supports more than 80 other Bioconductor packages, including those for sequence analysis, differential expression analysis and visualization.

This is a PLOS Computational Biology Software Article.

Introduction The genome is typically represented as a linear sequenceterm, split over multiple chromosomes, and data are linked to the genome by occupying a rangeterm of positions on the sequence. These data fall into two broad categories. First, there are the annotationsterm, such as gene models, transcription factor binding site predictions, GC percentage, polymorphisms, and conservation scores. Such annotationsterm are highly processed and are often served by public databases such as NCBI or EBI. Second, there are primary experimental measurementsterm, such as read alignmentsterm from high-throughput sequencing. Data integration, within and between those two categories, is made possible by treating the data as rangesterm on the genome, which acts as a common scaffoldterm. Thus, rangesterm play a central role in genomic data analysisterm·dom, and statistical tools should consider rangesterm to be as fundamental as quantitative and categorical data types.

For example, rangesterm are integral to the manipulation of gene model annotationsterm. Examples include deriving candidate promoter regions, finding introns, calculating the total exonic length of a transcript or finding the exonic regions that are unique to a particular transcript in an alternatively spliced gene. Rangesterm also play a central role in the analysis of experimental data, where they are used to represent read alignmentsterm. In the analysis of ChIP-seq data, it is typical to calculate the depth of alignment coverageterm, which then serves as input to calling algorithms which output peaks as rangesterm. These rangesterm are then annotated according to their overlap with and proximity to other rangesterm, such as gene structures. Similarly, for RNA-seq data, analysts measure gene expression based on counting the alignments overlapping exons.

All these analyses depend on specialized, range-based algorithms and data structures. For example, computations on gene models involve set operations on rangesterm, including intersectionterm, unionterm and complementterm. Coverageterm calculation is important for detecting regions of enrichment and for producing visual summaries. Overlap and nearest neighbor detectionterm is fundamental to the annotationterm of ChIP-seq peaks, estimating expression from RNA-seq data and many other integrative analyses.

The primary argument for storing rangesterm in specialized, formal data structuresterm·mth is efficiency, in terms of both implementation and language. The notion of rangesterm can be made explicit in the application programming interface (API)term·mth, permitting the expression of algorithms in a succinct and readable language that illustrates concepts instead of exposing implementation details. Another goal is interoperabilityterm: by using the same data structures, multiple routines, spread across different packages, can operate on the data without cumbersome conversions. Also, a data structure can be accessed through an abstractionterm·mth that hides the details of the optimized implementation, and this results in looser coupling between components. Together, these benefits lead to more robust, maintainable softwareterm.

Data structures should support the storage of per-range metadataterm, because genomic data is multivariate and consists of much more than the rangesterm alone. This enables the storage of gene identifiers and other symbols with the gene rangesterm, and the peak heights or confidence scores with the peak rangesterm. Some metadataterm merit special treatment, such as the chromosome nameterm·var and the strandterm·var. Also necessary is a data structure for storing summaries and processing results for a common set of rangesterm across multiple samples. Such a structure would hold, for example, the RNA-seq per-exon counts or a set of variant calls. Finally, there should be support for storing hierarchies of rangesterm, at least for one level of nesting, to represent, for example, the nesting of exons into transcripts. Whether it is appropriate to treat the exons as individual rangesterm or the transcript as a compound rangeterm depends on the use case; both should be supported.

These data structures are represented as classes, through which we communicate the formal definition of each data structure to the programming language. One benefit is that we can defer the regulation of data access and the tracking of data integrity to the language. In the case of functional object-oriented languages, there is another benefit: we can implement behaviors as methodsterm·mth on generic functionsterm·mth. A generic functionterm·mth is one that dispatches to a particular implementation, termed a methodterm·mth, based on the classes of passed arguments. This means that the same API will exhibit specialized behavior depending on the input. For example calling startterm·var on a rangeterm data structure would return the starting positions for the rangesterm, while calling the same function on a base R time-series object would behave differently.

This paper describes the infrastructure in Bioconductor [1] for the integrative statistical analysis of range-based genomic data. Main features include scalable data structuresterm·mth for annotated genomic rangesterm and genome-length vectors, and efficient algorithms for overlap detectionterm and other range operationsterm. The packages that form the core of the infrastructureterm include IRangesterm , GenomicRangesterm and GenomicFeaturesterm . Source code for the packages is included in the supplement, under Software S1 , S2 , and S3 , respectively. The IRangesterm package provides the fundamental rangeterm data structures and operations, while GenomicRangesterm builds upon it to add biological semantics to the metadataterm, including explicit treatment of sequence name and strandterm·var. Finally, GenomicFeaturesterm enables access to and manipulation of gene models and other annotationsterm. Together, these packages support more than 80 other packages in Bioconductor.

Other software toolsterm provide facilities for working with genomic rangesterm, e.g., bedtools [2] and cisGenome [3] . Those provide UNIX command-line interfaces and rely on common file formats (which are often incompletely specified) to interoperate with other tools, leading to workflows embodied as: collections of heterogeneous scripts, system dependencies and data files. Such workflows can be difficult to maintain and challenging to reproduce. In contrast, the Bioconductor infrastructureterm is tightly integrated with other R packages through in-memory data structures, while still supporting interaction with external tools. The Bioconductor package genomeIntervals provides data structures for representing genomic rangesterm and utilities, such as overlap detectionterm, that have much in common with the tools described here, but our tools are more extensive and have been more widely adopted.

Design and Implementation Working with Simple Rangesterm We use the term “rangeterm” to denote an ordered set of consecutive integersterm. A rangeterm is represented by a pair of integers satisfying . In Figure 1 , and correspond to the startterm·var and endterm·var columns, respectively. The “width” of a rangeterm is given by , so a rangeterm for a single integer (modeling, for example, a single nucleotide position) has .

Tabular (top) and visual (bottom) representation of the exons for the human KRAS gene, derived from the UCSC known gene annotationterm. In the table, the columns seqnames , startterm·var and endterm·var locate the exons in the genome. The strandterm·var column indicates the direction of transcription. The exons are grouped into transcripts by tx_id , and the exon IDs are given by exon_id . Virtually all genomic data sets fit this pattern: genomic location, followed by a series of columns, often including strandterm·var and/or score, that annotate that location. In the plot, the rectangles represent exonic regions, and the arrows represent the introns, as well as the strandterm·var.

The IRangesterm package, which is designed to be general and thus avoids biology-specific considerations, introduces the IRanges classterm·mth to represent a vector of rangesterm. The GenomicRangesterm package builds on IRangesterm to include biologically relevant features such as strandterm·var and sequence (e.g., chromosome) name.

In Figure 1 , we show a table of the exons of the human gene KRAS . The tx_id column indicates the transcripts to which each exon belongs. A single IRangesterm object can store those exon rangesterm, and this model is appropriate for per-exon analyses.

The IRanges classterm·mth supports the basic R vector API, including the length accessor, extraction and subsetting functions like [[ and [ , concatenation with c, etc. This will hold true for all vector-like objects in the rangeterm infrastructure.

The available range operationsterm are listed in Table 1 . The IRangesterm object supports direct manipulation of the startterm·var, endterm·var and width of the contained rangesterm. In applications, many of these operations follow recurrent patterns, and manipulating startterm·var and endterm·var directly can be needlessly tedious and error-prone. For this reason, shift , resize and similar frequently useful range operationsterm are provided. Rangesterm can be simplified and summarized with several functions, including rangeterm , reduce and disjoin . Figure 2 illustrates the latter two. It is often appropriate to conceive of an IRangesterm object as a mathematical set of integers, or, in the biological context, a set of nucleotide positions. gaps (complementterm), unionterm , intersect and setdiff support this notion. For example, taking the unionterm of two transcripts would yield the rangesterm covered by any KRAS exon. The flank function could be used to demarcate putative promoter regions of transcripts.

Summary of the Rangesterm API.

Illustration of the reduce and disjoin operations on the last exon from each of the KRAS transcripts.

A recurrent operation is overlap detectionterm; various instances are illustrated in Figure 3 . In later sections, we apply overlap counting for finding the percentage of ChIP-seq peaks that overlap a promoter, counting the number of RNA-seq reads for each transcript, and other tasks. The findOverlaps function uses an efficient interval tree algorithm [4] to detect overlaps between two IRangesterm objects, as well as the more complex range-based data structures introduced later. The algorithm supports several types of overlap, including those defined by Allen's Interval Algebra [5] . The one-time cost of constructing the interval tree is , and queries are performed in logarithmic time. In accordance with the vectorized semantics of R, if multiple queries are submitted, they are efficiently processed in batch, without restarting at the root of the tree for each query. The language of implementation is C, which avoids the potentially expensive iteration over the tree in R.

Illustration of overlap (top) and adjacency (bottom) relationships. The any mode detects hits with partial or complete overlap, while within requires that the query rangeterm represents a subregion of the subject rangeterm.

Working with Genomic Rangesterm The IRanges classterm·mth encodes only the startterm·var and endterm·var of rangesterm but not the chromosome, strandterm·var nor other information that is important in genomic applications. The GenomicRangesterm package adds biological semantics on top of IRangesterm . At its core is the GRanges classterm·mth. Each element of a GRangesterm·mth instance includes a chromosome identifier and strandterm·var designation. Each data set is associated with a particular, versioned reference genome sequenceterm consisting of a discrete set of chromosomes or contigs, along with their lengths, if known. The GRanges classterm·mth thus fully represents the data in Figure 1 , and encourages best-practices (e.g., tracking genome build) to minimize book-keeping errors.

The GRanges classterm·mth supports many of the same range operationsterm as IRangesterm and specializes them for genomic data. We achieved API specialization by implementing methodsterm·mth for both classes on the same generic functionsterm·mth. In general, we believe methodterm·mth specialization is an effective practice for providing the same interface on top of different data structures. Abstracting implementation details leads to user code that is more robust and easier to write and maintain.

The GRangesterm·mth methodsterm·mth give special consideration to the chromosome and, when appropriate, the strandterm·var. For instance, the findOverlaps generic has methodsterm·mth for both IRangesterm and GRangesterm·mth , and the GRangesterm·mth methodterm·mth is specifically able to take advantage of the chromosome information when detecting overlaps. Operations that depend on a notion of direction optionally consider strandterm·var. For example, the resize function will resize from the startterm·var or endterm·var of the rangesterm in a IRangesterm object. For a GRangesterm·mth object, resize will take the startterm·var to be the leftmost position for positive strandterm·var features and the rightmost position for negative strandterm·var features.

Some types of genomic data, for example gene models or aligned paired-end reads, have a hierarchical structure. To represent this, multiple GRangesterm·mth objects may be combined into a GRangesListterm·mth , where each GRangesterm·mth is considered a compound featureterm. GRangesListterm·mth groups transcripts by gene, groups exons by transcript, and represents read alignmentsterm, where each alignment consists of multiple segments separated by gaps. For example, we group the KRAS exons by transcript using a GRangesListterm·mth . A note on performance: although the user interface presents each element of a GRangesListterm·mth as a GRangesterm·mth , internally there is only a single GRangesterm·mth , along with an assoicated partitioning that forms the list elements.

For a GRangesListterm·mth , overlap detectionterm reports a hit at the element level, i.e., when any rangeterm within an element overlaps a query rangeterm. See Figure 4 for an illustration. This semantic is convenient, for example, when counting the total number of RNA-seq read pairs overlapping the exonic regions of each transcript. In that case, both the reads and the transcripts are GRangesListterm·mth objects.

Illustration of overlap computations between two GRangesListterm·mth objects. Each set of rectangles linked by solid lines represents a compound rangeterm, i.e., an element of the list. Rangesterm in the query (top) are being matched against rangesterm in the subject (bottom). The labels between them indicate the type of overlap (any, within, none).

Accessing Gene Models Recalling our KRAS gene model example, there are multiple models for representing transcript structuresterm, and the applicability of each depends on the use case. To support the analyst in asking a broad rangeterm of questions, there is a need for a mechanism that draws from a data source of gene annotationsterm and returns them in a variety of different data models. The GenomicFeaturesterm package fills that role by distilling multiple data sources into a single database schema and wrapping that database in an API that returns, for example, the exons grouped by gene, or the bounds for every transcript. The databases are implemented in SQLite and are thus accessible from environments outside of R. For reproducibility, a database may be encapsulated in a redistributable R package.

The database is represented by the TranscriptDbterm·mth class and stores the rangeterm of each exon, the coding rangeterm, the transcript ID, the gene ID, and metadataterm about the source of the transcript information. The GenomicFeaturesterm package provides an automated mechanism for constructing a TranscriptDbterm·mth object from tracks defined in the UCSC genome browser, Biomart, or GTF/GFF files. Bioconductor provides pre-built packages for the most widely adopted gene models, like the UCSC known gene annotationsterm on hg19. These packages follow a standard naming convention, e.g., TxDb.Hsapiens.UCSC.hg19.knownGene .

There are functions for performing common queries that return the exons, coding regions, and transcript boundaries as a GRangesterm·mth object. Transcript and gene-level groupings are preserved by GRangesListterm·mth objects. The rangesterm in Figure 1 were derived from the TxDb.Hsapiens.UCSC.hg19.knownGene package using the following call to the function exons: > library(“TxDb.Hsapiens.UCSC.hg19.knownGene”) > library(“org.Hs.eg.db”) > kras_gene <- org.Hs.egSYMBOL2EG$KRAS > kras_exons <- exons(TxDb.Hsapiens.UCSC.hg19.knownGene, + vals = list(gene_id = kras_gene), + columns = c(“tx_id”, “exon_id”)) To retrieve the exons corresponding to a particular transcript, such as transcript 48666 of KRAS, we call exonsBy , which returns a GRangesListterm·mth of exons grouped by transcript, and extract the element corresponding to the desired transcript identifier: > exonsByTx <- exonsBy(TxDb.Hsapiens.UCSC.hg19.knownGene) > krasA <- exonsByTx[[“48666”]] The contents of krasA are shown in Table 2 .

Contents of the krasA object, representing the exons in isoform A of KRAS.

Associating Annotationsterm with Rangesterm Figure 1 demonstrates how genomic data consist of both rangesterm and uni- or multivariate annotationsterm on those rangesterm. In that table, the annotationsterm are the exon ID and a variable grouping the exons into transcripts. If we had read alignmentsterm from an RNA-seq experiment, we might use countOverlapsterm to generate a read count for each exon. Other examples of annotationterm would include the reference and alternate bases for a single nucleotide variant (SNV) or the position weight matrix (PWM) score for a putative transcription factor binding site.

Every multi-element data structure in the IRangesterm suite supports the storage of per-element metadataterm: data about data. In this case, the metadataterm are the annotationsterm, and the primary dataterm are the rangesterm being annotated. This includes all of the data structures for storing rangesterm, such as IRangesterm, GRangesterm·mth and GRangesListterm·mth . The metadataterm are stored in a DataFrameterm·mth with as many rows as there are elements in the object. We introduce a DataFrameterm·mth class that behaves similarly to the base R data.frame , but supports storage of complex vector-like objects (e.g., a DNAStringSet , representing DNA sequences, or a GRangesterm·mth ) in columns.

Working with Coverageterm and Similar Vectors A common methodterm·mth of summarizing a genomic data set is to calculate the coverageterm, i.e., the number of features in the data set overlapping each position in the genometerm. This is useful in ChIP-seq analysis, where many peak detection methodsterm·mth operate on the coverageterm.

For this example, mouse genomic DNA was cross-linked with DNA-binding proteins, fragmented and precipitated with an antibody for CTCF. An antibody for GFP was used for the control. The CTCF and GFP samples were each sequenced in a single lane on an Illumina sequencer, which generated reads 35 nt in length [6] . We excluded the last 11 nt of each read due to insufficient quality, so the effective read length was 24 nt. The reads were aligned to the mm9 build of the mouse genome using MAQ. We parsed the MAQ output using the ShortRead package [7] , and the alignments for three chromosomes (chr10, chr11 and chr12) were extracted for use as a demonstration data set in the chipseq package [8] .

> library(“chipseq”) > data(“cstest”) > ctcfReads <- cstest$ctcf The ctcfReads objectterm, listed in Table 3 , is a GRangesterm·mth holding the read alignmentsterm from the CTCF sample. The GRangesterm·mth stores the chromosome namesterm·var, rangesterm, and strandterm·var for each alignment, as well as a list of chromosome namesterm·var and lengths for the mm9 genome. Tracking the chromosome information guards against errors that could arise, for example, from mixing data across genome assemblies.

Rangesterm for the first three reads in the ctcfReads objectterm, storing the read alignmentsterm for the CTCF sample.

Each read represents only 24 nt from one endterm·var of a fragment of DNA. Since it was the fragment, but not necessarily the sequenced region, that was cross-linked to CTCF, we need to consider the entire fragment when predicting binding sites. We assume that the fragment length was approximately 120 nt and call resize to extend our read rangesterm to fragment-sized rangesterm: > ctcfFragments <- resize(ctcfReads, 120) Note that the strandterm·var of the alignment was automatically taken into account.

The coverageterm function calculates the coverageterm for a set of rangesterm. We calculate the coverageterm on the CTCF fragments from our ChIP-seq data set as follows: > ctcfCoverage <- coverageterm(ctcfFragments) > ctcfCoverage10 <- ctcfCoverage$chr10 The ctcfCoverage object is a list, with one coverage vectorterm per chromosome. For simplicity, we extract the element corresponding to “chr10”. The ctcfCoverage10 object is of class Rleterm·mth .

Vectors along the genome tend to have many repeated values. For the sake of compactness, we compress the data using a run-length encoding compression scheme. Through the R class system, we abstract this efficient implementation behind an API that supports the features of ordinary R vectors; the complexity is hidden from the user. The Rleterm·mth class represents a run-length encoded vectorterm and provides features beyond those of ordinary vectors. For example, one can use rangesterm to extract values from an Rleterm·mth . This integrates range-based datasets with data in chromosome-length vectors. To demonstrate, we find the position of the maximum coverageterm on chr10, and, in order to display the coverageterm in context, we extract a 5000 nt region centered around that position: > maxPos <- which.max(ctcfCoverage10) > roi <- resize(IRangesterm(maxPos, width = 1), 5000, “center”) > roiCoverage <- ctcfCoverage$chr10[roi] The roiCoverage vector is plotted in Figure 5 and shows how the coverageterm relates to the gene context.

Visualization of the coverageterm of bases by GFP- and CTCF-bound fragments (top) in the context of part of the gene model for Rrp1, Entrez gene 18114 (bottom).

Next, we wish to find peaks in the coverage vectorterm by slicing it at a fixed threshold. We call the slice function and pass it a cutoff of 8: > ctcfPeaks <- slice(ctcfCoverage10, lower = 8) The resulting coverageterm slices are rangesterm and and we rely on our data structures for managing and manipulating them. In particular, the ctcfPeaks object is a Views object, which combines the peak rangesterm with the original Rleterm·mth coverage vectorterm (another example of integrating rangesterm with vectors). Views has several utilities for summarizing the vector values within each rangeterm. We use these to summarize the coverageterm values within each peak: > ctcfMaxs <- viewMaxs(ctcfPeaks) > ctcfSums <- viewSums(ctcfPeaks) Accessing Read Alignmentsterm Read alignmentsterm may be loaded from a BAM fileterm as a GappedAlignmentsterm·mth or GappedAlignmentPairs object, depending on whether the reads should be treated as paired. Both data structures store short read alignmentterm results in terms of the position, chromosome, strandterm·var, CIGAR string (a compact representation of the gaps) and other information. Both also support some of the range operationsterm and can be coerced to GRangesterm·mth and GRangesListterm·mth . The GRangesterm·mth representation holds the ungapped extents of the read alignmentsterm, whereas the GRangesListterm·mth represents the alignments as rangesterm with gaps, including the inter-read gap of a pair and the skipped regions in the reference (e.g., introns). The choice of data structure depends on, for example, whether one wants to count the overlapping pairs, reads, or aligned segments separated by gaps.

A frequent goal of RNA-seq experiments is to estimate the levels of gene expression. For this demonstration, we will simply count the number of read alignmentsterm overlapping the exonic regions of each transcript. We begin by loading a BAM fileterm of read alignmentsterm from an RNA-seq experiment in yeast [9] . There are four samples, two of which are wildtype and two of which are RLP mutants. The leeBamViews package provides the alignments on chromosome XIII from positions 800000 to 900000.

> bams <- getBamsFromLeeBamViews() > ga <- readGappedAlignments(bams [ 1 ] ) > reads <- grglist(ga) The readGappedAlignments function loads the BAM fileterm as a GappedAlignmentsterm·mth object, which is then coerced to a GRangesListterm·mth , where each read consists of one or more rangesterm, separated by intronic gaps.

In the next step, we obtain the transcript annotationsterm for yeast. Bioconductor provides a TranscriptDbterm·mth object for the UCSC sacCer2 assembly, and we extract from it a GRangesListterm·mth representing the transcripts, by calling exonsBy . In order to detect overlaps, the alignments and gene annotationsterm need to have the same chromosome namesterm·var; we correct for that with the calls to keepSeqlevels and renameSeqlevels .

> library(“TxDb.Scerevisiae.UCSC.sacCer2.sgdGene”) > tx <- exonsBy(TxDb.Scerevisiae.UCSC.sacCer2.sgdGene) > tx <- renameSeqlevels(keepSeqlevels(tx, “chrXIII”), + c(chrXIII = “Scchr13”)) Now that the data and annotationsterm have been loaded, we count the number of read alignmentsterm in each genomic feature: > counts <- countOverlapsterm(tx, reads, ignore.strandterm·var = TRUE) More complex counting algorithms are available via the summarizeOverlapsterm function, which counts over multiple samples and returns the results as a SummarizedExperimentterm·mth object. Unlike the call to countOverlapsterm above, reads that map to multiple features are discarded.

Summarized Experiments The typical workflow in a genomic data analysisterm·dom is to reduce a complex raw data set, such as a set of RNA-seq read alignmentsterm, to a set of summaries, such as the number of reads aligned over each transcript or exon. Analogously, an exome-seq experiment yields variant calls at particular genomic rangesterm. It is often appropriate and convenient to store these summaries as a matrix, where the rows correspond to genes or some other genomic feature and the columns to samples. The SummarizedExperimentterm·mth class is designed to hold such values, along with annotationsterm on the genomic regions, the samples and the experiment as a whole. Its design follows the same pattern as the ExpressionSetterm·mth in the Bioconductor microarray infrastructure; the primary difference is that SummarizedExperimentterm·mth is based on IRangesterm data structures. A SummarizedExperimentterm·mth may be constructed directly or generated by a function that executes a high-level workflow. For example, for a list of BAM filesterm, the summarizeOverlapsterm function counts the overlaps between the read alignmentsterm and a database of transcripts.

Results In this section we describe range-based integrative computations related to the genetics of protein-DNA binding, and conclude with a topically organized list of Bioconductor packages that make essential use of the rangesterm infrastructure.

Exploring Genetics of CTCF Binding In the murine ChIP-seq example discussed previously, the rangesterm infrastructure was used to compute and display variation in read coverageterm over the mouse genome. In this example, we use tools based on the rangesterm infrastructure to examine both coverageterm and content of reads from a larger ChIP-seq experiment on human cell lines. Our basic intent is to show how the infrastructure can be used to evaluate the roles of genotype and genetic diversity in the genomic sequences where CTCF is reported to bind, with an understanding that ultimate inferences on protein binding locations and on base-call distributions over heterozygous loci will need to directly incorporate risks of base calling and read mapping errors, and will need to be followed up with wet-lab validation.

For an investigation of the prevalence of allele-specific protein-DNA bindingterm [10] , BAM filesterm on 22 ChIP-seq experiments addressing CTCF binding to DNA from immortalized B-cells were collected from the ENCODE project portal. Exclusion of files with aberrant quality score distributions left 16 BAM filesterm corresponding to 12 distinct individuals; two technical replicates were available for each of four individuals. Furthermore, different base-call qualityterm score scales were used for two batches of samples; by subtracting 31 from the reported mean quality scores for one set of samples, approximately identical medians and interquartile rangesterm were established for mean quality scores for all retained samples. Base calls for reads bound at all genomic locations with positive coverageterm were tabulated using the VariantTools tallyVariants function, and reduced to locations exhibiting statistical evidence of allele-dependent CTCF binding using the callVariants function. Each of these variant assessment tools makes use of infrastructure derived from GSNAP [11] with key results materialized as GRangesterm·mth instances. Table 4 depicts an excerpt from a callVariants output.

Partial output of countVariants applied to a BAM fileterm from an ENCODE CTCF ChIP-seq experiment.

The genome-wide searches for allele-dependent CTCF binding events employed default settings for variant calling by VariantTools callVariants , which include criteria on minimum coverageterm, minimum diversity of read cycles at which base is found, and limitation of risk of strandterm·var bias. This process yielded a total of 19655 locations with evidence of allele-dependent CTCF binding, corresponding to 50750 events over the 12 individuals. We obtained GRangesterm·mth representations of dbSNP build 137 with the scanVcf function of the VariantAnnotation . This facilitated distributed computation for partitioning allele-dependent CTCF binding events into 12691 coincident with known polymorphisms and 6964 at locations where no SNP has been reported in dbSNP.

Allelic imbalanceterm·def in CTCF binding corresponds to departure of the alternate nucleotide proportionterm (ANPterm) at a CTCF binding site from 50%. The upper panels of Figure 6 show the distributions of ANPterm stratified by coincidence of allele-dependent CTCF binding locations with locations of known SNP. The lower panels show identically stratified associations between ANPterm and mean base-call qualities for pileups over the allele-dependent binding locations. While the off-SNP locations show a proponderance of ANPterm below 20%, there is also an indication that base-call qualityterm for such binding events is relatively low, implying that these findings would be unlikely to replicate. For example, among on-SNP allele-dependent binding calls with ANPterm below 20%, 5% had mean quality less than 5; the corresponding frequency for low-quality on-SNP allele-dependent binding locations calls was 41%.

Top panels: distributions of alternate nucleotide proportionsterm for on- and off-SNP allele-dependent CTCF binding events. Bottom panels: relationships between average call quality values and alternate nucleotide proportionsterm are depicted using a 2D density estimate (darker regions correspond to higher density.).

Figure 6 was constructed using the packages and features described above with very little programming effort needed to specialize the computations to this example. Most computational biologists and other interested investigators could easily carry out these steps. We can conclude that allelic imbalanceterm·def in CTCF binding events is frequently detectable but, for the data considered, the apparent imbalances observed are likely a mix of real biology and technical artifactsterm originating from, e.g., sequencing and read mapping errors. Careful analysis of metadataterm collected in the variant-calling process may help to disentangle the key factors contributing to allele-dependent CTCF binding.

Software Based on the Infrastructure There is a growing ecosystem of packagesterm based on this infrastructure. By current count, more than 80 packages depend directly on the packages presented here. This includes packages for input and output of rangesterm ( rtracklayer , Rsamtools ) [12] , [13] , quality assessment ( ShortRead ) [7] , sequence analysis ( Biostrings ) [14] , variant calling ( VariantTools ) [15] , and other tasks. To summarize the different use cases addressed by dependent packages, Table 5 tabulates the descriptive labels chosen from a controlled vocabulary by the package authors.

Selected packages based on the Rangesterm infrastructure.

Availability and Future Directions All of the packages described, including IRangesterm , GenomicRangesterm and GenomicFeaturesterm , form the core infrastructure for sequence analysisterm in Bioconductor and are available from the project website: http://bioconductor.org (see also Software S1 - S3 ). We aim to continue to support scientists in their drive to further science by asking increasingly complex and integrative questions about increasingly complex and heterogeneous data. For example, we are working towards better support for detecting alternative and novel splicing, measuring isoform-specific expression, annotating sequence variants, mapping between genome, transcript and protein coordinate spaces, and integrating transcript annotationsterm with gene-level metadataterm. There is also an unmet need in the visualization of genomic rangesterm. In particular, we need better visualizations for relating RNA-seq coverageterm and junction counts to transcript structuresterm, and for diagnosing read alignmentsterm in the context of variant calling. Finally, as datasets continue to expand in size, we continue to seek more efficient algorithms and data structures, and we are vigilant for opportunities to leverage parallel computingterm.

Supporting Information The IRangesterm package. The IRangesterm package provides efficient low-level and reusable S4 classes for storing and manipulating rangesterm of integers and compressed, genome-length vectors.

(GZ) Click here for additional data file.

The GenomicRangesterm package. The GenomicRangesterm package defines general purpose containers for storing genomic rangesterm as well as more specialized containers for storing alignments against a reference genome.

(GZ) Click here for additional data file.

The GenomicFeaturesterm package. The GenomicFeaturesterm package is a set of tools and methodsterm·mth for making and manipulating transcript-centric annotationsterm.

(GZ) Click here for additional data file.

Funding Statement This work was funded by the National Institutes of Health, National Human Genome Research Group through grants P41 HG004059 and U41 HG004059 and (for VJC) by National Heart, Lung and Blood Institute grants R01 HL086601, R01 HL093076 and R01 HL094635. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

References Associated Data Supplementary Materials The IRangesterm package. The IRangesterm package provides efficient low-level and reusable S4 classes for storing and manipulating rangesterm of integers and compressed, genome-length vectors.

(GZ) Click here for additional data file.

The GenomicRangesterm package. The GenomicRangesterm package defines general purpose containers for storing genomic rangesterm as well as more specialized containers for storing alignments against a reference genome.

(GZ) Click here for additional data file.

The GenomicFeaturesterm package. The GenomicFeaturesterm package is a set of tools and methodsterm·mth for making and manipulating transcript-centric annotationsterm.

(GZ) Click here for additional data file.