valr: Reproducible genome interval analysis in R

New tools for reproducible exploratory data analysis of large datasets are important to address the rising size and complexity of genomic data. We developed the valr R package to enable flexible and efficient genomic interval analysis. valr leverages new tools available in the ”tidyverse”, including dplyr. Benchmarks of valr show it performs similar to BEDtools and can be used for interactive analyses and incorporated into existing analysis pipelines.

This article is included in the gateway. RPackage v1 Introduction A routine bioinformatic task is the analysis of the relationships between sets of genomic intervals, including the identification of DNA variants within protein coding regions, annotation of regions enriched for nucleic acid binding proteins, and computation of read density within a set of exons. Command-line tools for interval analysis such as BEDtools 1 and BEDOPS 2 enable analyses of genome-wide datasets and are key components of analysis pipelines. Analyses with these tools commonly combine processing intervals on the command-line with visualization and statistical analysis in R. However, the need to master both the command-line and R hinders exploratory data analysis, and the development of reproducible research workflows built in the RMarkdown framework.
Existing R packages developed for interval analysis include IRanges 3 , bedr 4 , and GenometriCorr 5 . IRanges is a Bioconductor package that provides interval classes and methods to perform interval arithmetic, and is used by many Bioconductor packages. bedr is a CRAN-distributed package that provides wrapper R functions to call the BEDtools, BEDOPS, and tabix command-line utilities, providing out-of-memory support for interval analysis.
Finally, GenometriCorr provides a set of statistical tests to determine the relationships between interval sets using IRanges data structures. These packages provide functionality for processing and statistical inference of interval data, however they require a detailed understanding of S4 classes (IRanges) or the installation of external command-line dependencies (bedr). Additionally, these packages do not easily integrate with the recent advances provided by the popular tidyverse suite of data processing and visualization tools (e.g. dplyr, purrr, broom and ggplot2) 6 . We therefore sought to develop a flexible R package for genomic interval arithmetic built to incorporate new R programming, visualization, and interactivity features.

Methods
Implementation valr is an R package that makes extensive use of dplyr, a flexible and high-performance framework for data manipulation in R 7 . Additionally, compute intensive functions in valr are written in C++ using Rcpp to enable fluid interactive analysis of large datasets 8 . Interval intersections and related operations use an interval tree algorithm to efficiently search for overlapping intervals 9 . BED files are imported and handled in R as data_frame objects, requiring minimal pre or post-processing to integrate with additional R packages or command-line tools.
Operation valr is distributed as part of the CRAN R package repository and is compatible with Mac OS X, Windows, and major Linux operating systems. Package dependencies and system requirements are documented in the valr CRAN repository.

Use cases
To demonstrate the functionality and utility of valr, we present a basic tutorial for using valr and additional common use cases for genomic interval analysis.
Basic usage Input data. valr provides a set of functions to read BED, BEDgraph, and VCF formats into R as convenient tibble (tbl) data_frame objects. All tbls have chrom, start, and end columns, and tbls from multi-column formats have additional pre-determined column names. Standards methods for importing data (e.g. read. Example of combining valr tools. The functions in valr have similar names to their BEDtools counterparts, and so will be familiar to users of the BEDtools suite. Also, similar to pybedtools 10 , a python wrapper for BEDtools, valr has a terse syntax. For example, shown below is a demonstration of how to find all intergenic SNPs within 1 kilobase of genes using valr. The BED files used in the following examples are described in the Data Availability section.

library(dplyr)
snps <-read_bed(valr_example("hg19.snps147.chr22.bed.gz"), n_fields = 6) genes <-read_bed(valr_example("genes.hg19.chr22.bed.gz"), n_fields = 6) # find snps in intergenic regions intergenic <-bed_subtract(snps, genes) # distance from intergenic snps to nearest gene nearby <-bed_closest(intergenic, genes) Visual documentation. By conducting interval arithmetic entirely in R, valr is also an effective teaching tool for introducing interval analysis to early-stage analysts without requiring familiarity with both command-line tools and R. To aid in demonstrating the interval operations available in valr, we developed the bed_glyph() tool which produces plots demonstrating the input and output of operations in valr in a manner similar to those found in the BEDtools documentation. Shown below is the code required to produce glyphs displaying the results of intersecting x and y intervals with bed_intersect(), and the result of merging x intervals with bed_merge() ( Figure 1).
Column specification. Columns in BEDtools are referred to by position: # calculate the mean of column 6 for intervals in 'b' that overlap with 'a' bedtools map -a a.bed -b b.bed -c 6 -o mean In valr, columns are referred to by name and can be used in multiple name/value expressions for summaries.

Summarizing interval coverage across genomic features
This demonstration illustrates how to use valr tools to perform a "meta-analysis" of signals relative to genomic features. Here we analyze the distribution of histone marks surrounding transcription start sites, using H3K4Me3 Chip-Seq data from the ENCODE project.
First we load packages and relevant data.

Interval statistics
Estimates of significance for interval overlaps can be obtained by combining bed_shuffle(), bed_random() and the sample_ functions from dplyr with interval statistics in valr.
Here, we examine the extent of overlap of repeat classes (repeatmasker track obtained from the UCSC genome browser) with exons in the human genome (hg19 build, on chr22 only, for simplicity) using the jaccard similarity index. bed_jaccard() implements the jaccard test to examine the similarity between two sets of genomic intervals. Using bed_shuffle() and replicate() we generate a data_frame containing 100 sets of randomly selected intervals then calculate the jaccard index for each set against the repeat intervals to generate a null-distribution of jaccard scores. Finally, an empirical p-value is calculated from the null-distribution.

Benchmarking against bedtools
In order to ensure that valr performs fast enough to enable interactive analysis, key functionality is implemented in C++. To test the speed of major valr functions we generated two data_frames containing 1 million randomly selected 1 kilobase intervals derived from the human genome (hg19). Most of the major valr functions complete execution in less than 1 second, demonstrating that valr can process large interval datasets efficiently ( Figure 3A).
We also benchmarked major valr functions against corresponding commands in BEDtools. valr operates on data_frames already loaded into RAM, whereas BEDtools performs file-reading, processing, and writing. To compare valr against BEDtools we generated two BED files containing 1 million randomly selected 1 kilobase intervals derived from the human genome (hg19). For valr functions, we timed reading the table into R (e.g. with read_bed()) and performing the respective function. For BEDtools commands we timed executing the command with the output written to /dev/null. valr functions performed similarly or faster than BEDtools commands, with the exception of bed_map and bed_fisher ( Figure 3B).

Reproducible reports and interactive visualizations
Command-line tools like BEDtools and bedops can be incorporated into reproducible workflows (e.g., with snakemake 11 ), but it is cumbersome to transition from command-line tools to exploratory analysis and plotting software. RMarkdown documents are plain text files, amenable to version control, which provide an interface to generate feature rich PDF and HTML reports that combine text, executable code, and figures in a single document. valr can be used in RMarkdown documents to provide rapid documentation of exploratory data analyses and generate reproducible work-flows for data processing. Moreover, new features in RStudio, such as notebook viewing, and multiple language support enable similar functionality to another popular notebook platform jupyter notebooks.
Additionally, valr seamlessly integrates into R shiny 12 applications allowing for complex interactive visualizations relating to genomic interval analyses. We have developed a shiny application (available on Gitub) that explores ChiP-Seq signal density surrounding transcription start sites and demonstrates the ease of implementing valr to power dynamic visualizations.

Summary
valr provides a flexible framework for interval arithmetic in R/Rstudio. valr functions are written with a simple and terse syntax that promotes flexible interactive analysis. Additionally by providing an easy-to-use interface for interval arithmetic in R, valr is also a useful teaching tool to introduce the analyses necessary to investigate correlations between genomic intervals, without requiring familiarity with the command-line. We envision that valr will help researchers quickly and reproducibly analyze genome interval datasets.

Data and software availability
The valr package includes external datasets stored in the inst/extdata/ directory that were used in this manuscript. These datasets were obtained from the ENCODE Project 13 or the UCSC genome browser 14 . BED files were generated by converting the UCSC tables into BED format. BED and BEDgraph data was only kept from chromosome 22, and was subsampled to produce file sizes suitable for submission to the CRAN repository. The original raw data is available from the following sources:  valr can be installed via CRAN using install.packages("valr").

Are sufficient details of the code, methods and analysis (if applicable) provided to allow replication of the software development and its use by others? Yes
Is sufficient information provided to allow interpretation of the expected output datasets and any results generated using the tool? Yes Are the conclusions about the tool and its performance adequately supported by the findings presented in the article? Yes No competing interests were disclosed.

Competing Interests:
I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard. Riemondy and colleagues have made a valuable tool, `valr`, available to the general public, one which deserves much merit for bringing modern R idioms from the `tidyverse` into the world of genomic research. As a fellow bioinformatician who has struggled with the idiosyncrasies of the aforementioned tools for interval manipulation in R, `valr` addresses many of the usability issues associated with these legacy methods by fundamentally altering the user experience. Consequently, there are two main advantages to adopting `valr` for interval manipulation in R: ease of writing, and ease of reading, code. Thus, this referee wholeheartedly endorses `valr`, and hopes to see more work that brings many of thè tidyverse` philosophies over to working with genomics in R.
Nonetheless, there is some room for improvement of the associated manuscript to better help explain the philosophy and usage of `valr`, and its place amongst the many tools for the manipulation of genomic intervals.
Firstly, in the introduction, it is mentioned that there exist `IRanges` methods that utilize the S4 convention, whereas `valr` utilizes a less formal schema where 3 columns, `chrom`, `start`, `end`, are present in the `data_frame` object. Indeed, it may be of use to expand upon such design choices that were made, and what advantages/disadvantages are made in using this less formal schema, and any other highly pertinent choices that affect user experience. In addition, one line mentions integration with other `tidyverse` tools, and should expand upon this with either one to a few specific examples or explain this point in more detail. Additionally, it should also be pointed out how `valr` builds upon these existing toolkits, and either expands upon/adopts their conventions. One way might be to create a table comparing functions between `valr`/bedtools/GenomicRanges might be helpful for a reader to see that the toolkit will