\documentclass{article} %\VignetteIndexEntry{rtracklayer} \title{The \textbf{rtracklayer} package} \author{Michael Lawrence} \begin{document} \maketitle <>= options(width=70) @ \section{Introduction} The \textbf{rtracklayer} package is an interface (or \emph{layer}) between \textbf{R} and genome browsers. Its main purpose is the visualization of genomic annotation \emph{tracks}, whether generated through experimental data analysis performed in R or loaded from an external data source. The features of \textbf{rtracklayer} may be divided into two categories: 1) the representation and import/export of track data and 2) the control and querying of external genome browser sessions and views. For working with track data, the package defines a data structure named \textit{trackSet}, which extends the Bioconductor \textit{eSet} class. A \textit{trackSet} contains information on the features of a track, including their chromosome, start and end positions in the genome, as well as any associated experimental measurements. Import and export from and to files is supported in various formats, see Section~\ref{sec:import-export}. The \textbf{rtracklayer} package currently interfaces with the \textbf{UCSC} web-based genome browser as well as the locally installed Java-based \textbf{Argo} browser. Other packages may provide drivers for other genome browsers through a plugin system. With \textbf{rtracklayer}, the user may start a genome browser session, create and manipulate genomic views, and import/export tracks and sequences to and from a browser. Please note that not all features are necessarily supported by every browser interface. The rest of this vignette will consist of a number of case studies. First, we consider an experiment investigating microRNA regulation of gene expression, where the microRNA target sites are the primary genomic features of interest. \section{Gene expression and microRNA target sites} This section will demonstrate the features of \textbf{rtracklayer} on a microarray dataset from a larger experiment investigating the regulation of human stem cell differentiation by microRNAs. The transcriptome of the cells was measured before and after differentiation by HG-U133plus2 Affymetrix GeneChip arrays. We begin our demonstration by constructing an annotation dataset from the experimental data, and then illustrate the use of the genome browser interface to display interesting genomic regions in the UCSC browser. \subsection{Creating a target site track} For the analysis of the stem cell microarray data, we are interested in the genomic regions corresponding to differentially expressed genes that are known to be targeted by a microRNA. We will represent this information as an annotation track, so that we may view it in the UCSC genome browser. \subsubsection{Constructing the \textit{trackSet}} In preparation for creating the microRNA target track, we first used \textbf{limma} to detect the differentially expressed genes in the microarray experiment. The locations of the microRNA target sites were obtained from MiRBase. Information about the target sites on differentially expressed genes was stored in the \textit{data.frame} called \texttt{targets}. % <>= library("humanStemCell") data(fhesc) library("genefilter") filtFhesc <- nsFilter(fhesc)[[1]] library("limma") design <- model.matrix(~filtFhesc$Diff) hesclim <- lmFit(filtFhesc, design) hesceb <- eBayes(hesclim) tab <- topTable(hesceb, coef = 2, adjust.method = "BH", n = 7676) tab2 <- tab[(tab$logFC > 1) & (tab$adj.P.Val < 0.01),] affyIDs <- tab2$ID library("microRNA") data(hsTargets) library("hgu133plus2.db") entrezIDs <- mappedRkeys(hgu133plus2ENTREZID[affyIDs]) library("org.Hs.eg.db") mappedEntrezIDs <- entrezIDs[entrezIDs %in% mappedkeys(org.Hs.egENSEMBLTRANS)] ensemblIDs <- mappedRkeys(org.Hs.egENSEMBLTRANS[mappedEntrezIDs]) targetMatches <- match(ensemblIDs, hsTargets$target, 0) targets <- hsTargets[targetMatches,] @ % The following code creates the track from the \texttt{targets} dataset by passing it to the \texttt{trackSet} constructor that is contained in the \textbf{rtracklayer} package. The columns of \texttt{targets} are named according to the convention expected by the constructor. % <>= head(targets) library(rtracklayer) targetTrack <- trackSet(targets) @ % \subsubsection{Accessing track information} The track information is now stored in the R session as a \textit{trackSet} object. As the \textit{trackSet} class inherits from \textit{eSet} defined in the \textbf{Biobase} package, it has slots for storing feature information (\textit{featureData}), experimental design information (\textit{phenoData}) and experimental measurements (\textit{assayData}). Most of the information for each feature in a track is stored in the \textit{featureData}. This includes the chromosome name, numeric start and end positions, the DNA strand (+/-/NA), and any other information available. There are methods for accessing the most commonly used fields: \texttt{chrom}, \texttt{start}, \texttt{end}, and \texttt{strand}. For example, the following code retrieves the chromosome names and then start positions for each feature in the track. <>= head(chrom(targetTrack)) head(start(targetTrack)) @ \paragraph{Exercises} \begin{enumerate} \item Get the strand of each feature in the track \item Calculate the length of each feature \item Reconstruct (partially) the \texttt{targets} \textit{data.frame} \end{enumerate} <>= head(strand(targetTrack)) head(end(targetTrack) - start(targetTrack) + 1) data.frame(chrom = chrom(targetTrack), start = start(targetTrack), end = end(targetTrack), strand = strand(targetTrack)) @ \subsubsection{Subsetting a \textit{trackSet}} It is often helpful to extract subsets from \textit{trackSet} instances, especially when uploading to a genome browser. Subsets are specified by a set of features, as well as a set of columns in the \textit{dataVals} matrix, if any. The conventional \texttt{[} method is employed for subsetting, where the first parameter, \textit{i}, indexes the features and \textit{j} indexes the data columns. Both \textit{i} and \textit{j} may contain numeric, logical and character indices, which behave as expected. % <>= ## get the first 10 targets first10 <- targetTrack[1:10] ## get pos strand targets posTargets <- targetTrack[strand(targetTrack) == "+"] @ % The features may also be indexed by a \textit{genomeSegment} instance (explained later) or by chromosome. We give an example of the latter below for obtaining all target sites on chromosome 1. In order to distinguish between a chromosome name and a feature name, we cast the chromosome name to a \textit{chrid} instance with the \texttt{chrid} function. % <>= chr1Targets <- targetTrack[chrid("chr1")] @ % \paragraph{Exercises} \begin{enumerate} \item Subset the track for all features on the negative strand of chromosome 2. \end{enumerate} <>= negTargetTrack <- targetTrack[strand(targetTrack) == "-"] @ \subsubsection{Exporting and importing tracks} \label{sec:import-export} Import and export of \textit{trackSet} instances is supported in the following formats: Browser Extended Display (BED), versions 1, 2 and 3 of the General Feature Format (GFF), and Wiggle (WIG). Support for additional formats may be provided by other packages through a plugin system. To save the microRNA target track created above in a format understood by other tools, we could export it as BED. This is done with the \texttt{export} function, which accepts a filename or any R connection object as its target. If a target is not given, the serialized string is returned. The desired format is derived, by default, from the extension of the filename. Use the \texttt{format} parameter to explicitly specify a format. <>= export(targetTrack, "targets.bed") @ To read the data back in a future session, we could use the \texttt{import} function. The source of the data may be given as a connection, a filename or a character vector containing the data. Like the \texttt{export} function, the format is determined from the filename, by default. <>= restoredTrack <- import("targets.bed") @ \paragraph{Exercises} \begin{enumerate} \item Output the track to a file in the ``gff'' format. \item Read the track back into R. \item Export the track as a character vector. \end{enumerate} <>= export(targetTrack, "targets.gff") targetGff <- import("targets.gff") targetChar <- export(targetTrack, format = "gff1") @ \subsection{Viewing the targets in a genome browser} For the next step in our example, we will load the track into a genome browser for visualization with other genomic annotations. The \textbf{rtracklayer} package is capable of interfacing with any genome browser for which a driver exists. In this case, we will interact with the web-based \textbf{UCSC} browser, but the same code should work for any browser. \subsubsection{Starting a session} The first step towards interfacing with a browser is to start a browser session, represented in R as a \textit{browserSession} object. A \textit{browserSession} is primarily a container of tracks and genomic views. The following code creates a \textit{browserSession} for the \textbf{UCSC} browser: <>= session <- browserSession("ucsc") @ Note that the name of any other supported browser could have been given here instead of ``ucsc''. To see the names of supported browsers, enter: <>= genomeBrowsers() @ \subsubsection{Laying the track} Before a track can be viewed on the genome, it must be loaded into the session using the \texttt{layTrack} function, as demonstrated below: <>= session <- layTrack(session, targetTrack, name = "targets") @ The \textit{name} argument should be a character vector that will help identify the track within \texttt{session}. Note that the invocation of \texttt{layTrack} above does not specify an upload format. Thus, the default, ``auto'', is used. Since the track does not contain any data values, the track is uploaded as BED. To make this explicit, we could pass ``bed'' as the \textit{format} parameter. \paragraph{Exercises} \begin{enumerate} \item Lay a track with the first 100 features of \texttt{targetTrack} \end{enumerate} <>= session <- layTrack(session, targetTrack[1:100], name = "target100") @ \subsubsection{Viewing the track} By default, the \texttt{layTrack} function creates a \textit{browserView}, an object that represents a view of a particular set of tracks along a particular region of the genome. For \textbf{UCSC}, this roughly corresponds to one tab or window in the web browser. To override the automatic creation of a view, pass \texttt{view = FALSE} to the \texttt{layTrack} function. The default view of the track attempts to show the entire track. The view region is determined by a call to \texttt{genomeSegment} on the track object. <>= genomeSegment(chr1Targets) @ The returned value from \texttt{genomeSegment} is an instance of the \textit{genomeSegment} class, which specifies a segment of a genome by its genome name, chromosome name and start and end positions. If any of the fields are empty, they stand for wildcards. \paragraph{Exercises} \begin{enumerate} \item Get the \textit{genomeSegment} describing the first feature in \texttt{targetTrack}. \item Get the chromosome ID of that segment. \end{enumerate} <>= segment <- genomeSegment(targetTrack[1]) chrom(segment) @ The target sites are distributed throughout the genome, so we will only be able to view a few features at a time. In this case, we will view only the first feature in the track. A convenient way to focus a view on a particular set of features is to subset the track and pass the range of the subtrack to the constructor of the view. Below we take a track subset that contains only the first feature. <>= subTargetTrack <- targetTrack[1] # get first feature @ Now we call the \texttt{browserView} function to construct the view and pass the \textit{genomeSegment} of the subtrack, zoomed out by a factor of 10, as the segment to view. By passing the name of the targets track in the \textit{pack} parameter, we instruct the browser to use the ``pack'' mode for viewing the track. This results in the name of the microRNA appearing next to the target site glyph. <>= view <- browserView(session, genomeSegment(subTargetTrack) / 10, pack = "targets") @ \paragraph{Exercises} \begin{enumerate} \item Create a new view with the same region as \texttt{view}, except zoomed out 2X. \item Create a view with the ``targets'' track displayed in ``full'' mode, instead of ``packed''. \end{enumerate} <>= viewOut <- browserView(session, genomeSegment(view) / 2) viewFull <- browserView(session, full = "targets") @ \subsubsection{A shortcut} There is also a shortcut to the above steps. The \texttt{browseGenome} function creates a session for a specified browser, loads one or more tracks into the session and creates a view of a given genome segment. In the following code, we create a new \textbf{UCSC} session, load the track and view the first two features, all in one call: <>= browseGenome(tracks = targetTrack, segment = genomeSegment(subTargetTrack) / 10) @ It is even simpler to view the subtrack in \textbf{UCSC} by relying on parameter defaults: <>= browseGenome(subTargetTrack) @ \subsubsection{Downloading tracks} It is possible to query the browser to obtain the names of the loaded tracks and to download the tracks into R. To list the tracks loaded in the browser, enter the following: % <>= loaded_tracks <- tracks(session) @ % One may download any of the tracks, such as the ``targets'' track that was loaded previously in this example. % <>= subTargetTrack <- trackSet(session, "targets") @ % By default, the segment of the track downloaded is the current default genome segment associated with the session. One may download track data for any genome segment, such as those on a particular chromosome. % <>= posTargets <- trackSet(session, "targets", genomeSegment(chr1Targets)) @ \paragraph{Exercises} \begin{enumerate} \item Get the SNP under the first target, displayed in \texttt{view}. \item Get the UCSC gene for the same target. \end{enumerate} <>= targetSNP <- trackSet(session, "snp128", genomeSegment(subTargetTrack)) featureNames(targetSNP) targetGene <- trackSet(session, "knownGene", genomeSegment(subTargetTrack)) featureNames(targetGene) @ \subsubsection{Accessing view state} The \texttt{view} variable is an instance of \textit{browserView}, which provides an interface for getting and setting view attributes. Note that for the UCSC browser, changing the view state opens a new view, as a new page must be opened in the web browser. To programmatically query the segment displayed by a view, use the \texttt{genomeSegment} method for a \textit{browserView}. % <>= segment <- genomeSegment(view) @ % Similarly, one may get and set the names of the visible tracks in the view. <>= visible_tracks <- tracks(view) tracks(view) <- visible_tracks @ The visibility mode (hide, dense, pack, squish, full) of the tracks may be retrieved with the \texttt{ucscTrackModes} method. % <>= modes <- ucscTrackModes(view) @ % The returned value, \texttt{modes}, is of class \textit{ucscTrackModes}. The modes may be accessed using the \texttt{[} function. Here, we set the mode of our ``targets'' track to ``full'' visibility. <>= modes["targets"] modes["targets"] <- "full" ucscTrackModes(view) <- modes @ Existing browser views for a session may be retrieved by calling the \texttt{browserViews} method on the \textit{browserSession} instance. % <>= views <- browserViews(session) length(views) @ % \paragraph{Exercises} \begin{enumerate} \item Retrieve target currently visible in the view. \item Limit the view to display only the SNP, UCSC gene and target track. \item Hide the UCSC gene track. \end{enumerate} <>= viewTarget <- trackSet(session, "targets", genomeSegment(view)) tracks(view) <- c("snp128", "knownGene", "targets") ucscTrackModes(view)["knownGene"] <- "hide" @ \section{CPNE1 expression and HapMap SNPs} Included with the \textbf{rtracklayer} package is a track object (created by the \textbf{GGtools} package) with features from a subset of the SNPs on chromosome 20 from 60 HapMap founders in the CEU cohort. Each SNP has an associated data value indicating its association with the expression of the CPNE1 gene according to a Cochran-Armitage 1df test. The top 5000 scoring SNPs were selected for the track. We load the track presently. <>= library(rtracklayer) data(cpneTrack) @ \subsection{Loading and manipulating the track} The data values for a track are stored as a numeric matrix under the name \textit{dataVals} within the \textit{assayData} slot. The values may retrieved with the \texttt{dataVals} method. <>= head(dataVals(cpneTrack)) @ % Sometimes, it may be convenient to extract the track information as a % \textit{data.frame}. The \textit{trackData} function does this by % combining the \textit{featureData} matrix with the \textit{dataVals}. % It also adds a column named \textit{featMid}, which gives the % mid-points (the mean of the start and end positions) of each feature % in the track. Here is an example of using \textit{trackData} to plot % the test value for each SNP vs. its position. One use of extracting the data values is to plot the data. <>= plot(start(cpneTrack), dataVals(cpneTrack)) @ \subsection{Browsing the SNPs} We now aim to view some of the SNPs in the UCSC browser. Unlike the microRNA target site example above, this track has quantitative information, which requires special consideration for visualization. \subsubsection{Laying a WIG track} To view the SNP locations as a track in a genome browser, we first need to upload the track to a fresh session. <>= session <- browserSession() session <- layTrack(session, cpneTrack, name = "cpne") @ % Note that because \texttt{cpneTrack} contains data values and its features do not overlap, it is uploaded to the browser in the WIG format. One limitation of the WIG format is that it is not possible to encode strand information. Thus, each strand needs to have its own track, and \textbf{rtracklayer} does this automatically, unless only one strand is represented in the track (as in this case). One could pass ``bed'' to the \textit{format} parameter of \texttt{layTrack} to prevent the split, but tracks uploaded as BED are much more limited compared to WIG tracks in terms of visualization options. To form the labels for the WIG subtracks, `` p'' is concatenated onto the plus track and `` m'' onto the minus track. Features with missing track information are placed in a track named with the `` na'' postfix. It is important to note that the subtracks must be identified individually when, for example, downloading the track or changing track visibility. \subsubsection{Plotting the SNP track} To plot the data values for the SNP's in a track, we need to create a \textit{browserView}. We will view the first 5 SNPs in the track, which will be displayed in the ``full'' mode. % <>= view <- browserView(session, genomeSegment(cpneTrack[1:5]), full = "cpne") @ % The UCSC browser will plot the data values as bars. There are several options available for tweaking the plot, as described in the help for the \textit{wigTrackLine} class. These need to be specified laying the track, so we will lay a new track named ``cpne2''. First, we will turn the \textit{autoScale} option off, so that the bars will be scaled globally, rather than locally to the current view. Then we could turn on the \textit{yLineOnOff} option to add horizontal line that could represent some sort of cut-off. The position of the line is specified by \textit{yLineMark}. We set it arbitrarily to the 25\% quantile. % <>= session <- layTrack(session, cpneTrack, name = "cpne2", autoScale = FALSE, yLineOnOff = TRUE, yLineMark = quantile(dataVals(cpneTrack), .25)) view <- browserView(session, genomeSegment(cpneTrack[1:5]), full = "cpne2") @ % \section{Binding sites for NRSF} \label{sec:binding} Another common type of genomic feature is transcription factor binding sites. Here we will use the \textbf{Biostrings} package to search for matches to the binding motif for NRSF, convert the result to a track, and display a portion of it in the \textbf{UCSC} browser. \subsection{Creating the binding site track} We will use the \textbf{Biostrings} package to search human chromosome 1 for NRSF binding sites. The binding sequence motif is assumed to be \textit{TCAGCACCATGGACAG}, though in reality it is more variable. To perform the search, we construct a \textit{PDict} dictionary and run \texttt{matchPDict} on the positive strand of chromosome 1. % <>= library(BSgenome.Hsapiens.UCSC.hg18) nrsfPD <- PDict("TCAGCACCATGGACAG") # TCAGCACC trusted band nrsfHits <- matchPDict(nrsfPD, Hsapiens[[1]]) countIndex(nrsfHits) @ % We then convert the hits, stored as an \textit{MIndex} object, to a \textit{trackSet} instance, via an \textit{IRanges} object. % <>= nrsfTrack <- trackSet(unlist(nrsfHits), "chr1", "+") @ % \subsection{Browsing the binding sites} Now that the NRSF binding sites are stored as a track, we can upload them to the UCSC browser and view them. Below, load the track and we view the region around the first hit in a single call to \texttt{browseGenome}. % <>= session <- browseGenome(nrsfTrack, segment = genomeSegment(nrsfTrack[1]) / 10) @ % We observe significant conservation across mammal species in the region of the motif. \section{Conclusion} These case studies have demonstrated a few of the most important features of \textbf{rtracklayer}. Please see the package documentation for more details. The following is the session info that generated this vignette: <>= sessionInfo() @ \end{document}