epivizr
: interactive visualization for genomic dataEpiviz
is an interactive visualization tool for functional genomics data. It supports genome navigation like other genome browsers, but allows multiple visualizations of data within genomic regions using scatterplots, heatmaps and other user-supplied visualizations. It also includes data from the Gene Expression Barcode project for transcriptome visualization. It has a flexible plugin framework so users can add d3 visualizations. You can find more information about Epiviz at http://epiviz.cbcb.umd.edu/help and see a video tour here.
The epivizr
package implements two-way communication between the R/Bioconductor
computational genomics environment and Epiviz
. Objects in an R
session can be displayed as tracks or plots on Epiviz. Epivizr uses Websockets for communication between the browser Javascript client and the R environment, the same technology underlying the popular Shiny system for authoring interactive web-based reports in R.
In this vignette we will look at colon cancer methylation data from the TCGA project and expression data from the gene expression barcode project. The epivizrData
package contains human chromosome 11 methylation data from the Illumina 450kHumanMethylation beadarray processed with the minfi
package. We use expression data from the antiProfilesData
bioconductor package.
library(epivizr)
library(antiProfilesData)
library(SummarizedExperiment)
data(tcga_colon_blocks)
data(tcga_colon_curves)
data(apColonData)
The tcga_colon_blocks
object is a GRanges
object containing chromosome 11 regions of hypo or hyper methylation in colon cancer identified using the blockFinder
function in the minfi
package.
show(tcga_colon_blocks)
## GRanges object with 129 ranges and 11 metadata columns:
## seqnames ranges strand | value
## <Rle> <IRanges> <Rle> | <numeric>
## [1] chr11 [ 4407026, 6435089] * | -0.142954896408814
## [2] chr11 [131239366, 133716186] * | -0.135137323912239
## [3] chr11 [ 55041873, 57022542] * | -0.17334711677526
## [4] chr11 [114645223, 116602403] * | -0.140934365473709
## [5] chr11 [ 78357700, 80184550] * | -0.15603691753532
## ... ... ... ... . ...
## [125] chr11 [29644815, 29650449] * | -0.0940450729648221
## [126] chr11 [41943963, 41956273] * | -0.0760193132782017
## [127] chr11 [16298618, 16314417] * | -0.0748443876820716
## [128] chr11 [38740154, 38740154] * | -0.117248761155248
## [129] chr11 [81757203, 81757203] * | -0.0710557720226853
## area cluster indexStart indexEnd L
## <numeric> <numeric> <integer> <integer> <numeric>
## [1] 30.3064380386685 495 130755 131000 212
## [2] 23.9193063324663 520 141959 142173 177
## [3] 19.5882241956043 507 134251 134374 113
## [4] 14.0934365473709 520 140000 140138 100
## [5] 14.0433225781788 507 138132 138249 90
## ... ... ... ... ... ...
## [125] 0.188090145929644 497 132733 132734 2
## [126] 0.152038626556403 502 133379 133380 2
## [127] 0.149688775364143 495 131986 131987 2
## [128] 0.117248761155248 499 133331 133331 1
## [129] 0.0710557720226853 508 138263 138263 1
## clusterL p.value fwer p.valueArea fwerArea
## <integer> <numeric> <numeric> <numeric> <numeric>
## [1] 1629 0 0 0 0
## [2] 1759 0 0 0 0
## [3] 1816 0 0 0 0
## [4] 1759 0 0 0 0
## [5] 1816 0 0 0 0
## ... ... ... ... ... ...
## [125] 367 0.0949029198604193 1 0.416952488890214 1
## [126] 45 0.316487220018491 1 0.485758597035402 1
## [127] 1629 0.342091920427093 1 0.494228876494974 1
## [128] 7 0.0564138506964121 1 0.592024814339825 1
## [129] 22 0.729905454979272 1 0.860076351814847 1
## -------
## seqinfo: 23 sequences from an unspecified genome; no seqlengths
The columns value
and p.value
can be used to determine which of these regions, or blocks, are interesting by looking at a volcano plot for instance.
plot(tcga_colon_blocks$value, -log10(tcga_colon_blocks$p.value), main="Volcano plot", xlab="Avg. methylation difference", ylab="-log10 p-value",xlim=c(-.5,.5))
The tcga_colon_curves
object is another GRanges
object which contains the basepair resolution methylation data used to define these regions.
show(tcga_colon_curves)
## GRanges object with 7135 ranges and 7 metadata columns:
## seqnames ranges strand | id type
## <Rle> <IRanges> <Rle> | <numeric> <array>
## [1] chr11 [131996, 132411] * | 129466 OpenSea
## [2] chr11 [189654, 189654] * | 129467 OpenSea
## [3] chr11 [190242, 190242] * | 129468 OpenSea
## [4] chr11 [192096, 192141] * | 129469 OpenSea
## [5] chr11 [192763, 193112] * | 129470 OpenSea
## ... ... ... ... . ... ...
## [7131] chr11 [134892703, 134892703] * | 142360 OpenSea
## [7132] chr11 [134903175, 134903175] * | 142361 OpenSea
## [7133] chr11 [134910774, 134910774] * | 142362 OpenSea
## [7134] chr11 [134911302, 134911302] * | 142363 OpenSea
## [7135] chr11 [134945848, 134945848] * | 142364 OpenSea
## blockgroup diff smooth normalMean cancerMean
## <numeric> <matrix> <numeric> <numeric> <numeric>
## [1] 495 -0.088871013 -0.12366708 0.75449261 0.66562159
## [2] 495 -0.001244860 -0.08656241 0.95952884 0.95828398
## [3] 495 -0.164759184 -0.08625123 0.68006830 0.51530911
## [4] 495 0.003094723 -0.08526643 0.05134196 0.05443669
## [5] 495 -0.102463563 -0.08484055 0.55235839 0.44989482
## ... ... ... ... ... ...
## [7131] 520 -0.11394759 -0.1421404 0.4013819 0.2874343
## [7132] 520 -0.01148901 -0.1530226 0.8482137 0.8367247
## [7133] 520 -0.20066455 -0.1627672 0.6365234 0.4358588
## [7134] 520 -0.14013791 -0.1635050 0.5831618 0.4430239
## [7135] 520 -0.24670370 -0.2308602 0.6757470 0.4290433
## -------
## seqinfo: 24 sequences from an unspecified genome; no seqlengths
This basepair resolution data includes mean methylation levels for normal and cancer and a smoothed estimate of methylation difference. This smoothed difference estimate is used to define regions in the tcga_colon_blocks
object.
Finally, the apColonData
object is an ExpressionSet
containing gene expression data for colon normal and tumor samples for genes within regions of methylation loss identified this paper. Our goal in this vignette is to visualize this data as we browse the genome.
The connection to Epiviz
is managed through a session manager object of class EpivizApp
. We can create this object and open Epiviz
using the startEpiviz
function.
app <- startEpiviz(workspace="qyOTB6vVnff", gists="2caf8e891201130c7daa")
This opens a websocket connection between the interactive R
session and the browser client. This will allow us to visualize data stored in the Epiviz
server along with data in the interactive R
session.
Windows users: In Windows platforms we need to use the service
function to let the interactive R
session connect to the epiviz
web app and serve data requests. We then escape (using ctl-c
or esc
depending on your environment) to continue with the interactive R
session. This is required anytime you want epivizr
to serve data to the web app, for example, when interacting with the UI. (We are actively developing support for non-blocking sessions in Windows platforms).
app$server$service()
Once the browser is open and is connected to an active epivizr session, the epivizr session registers available chart types supported in the epiviz JS app. To list available chart types we use method list_chart_types
. In this vignette, we only show three chart types available for use in the epiviz JS app. The production epiviz JS app has a larger number of chart types available. Chart types can also be added dynamically as described in the epiviz JS documentation: http://epiviz.cbcb.umd.edu/help. Chart types that are added dynamically are listed and usable within an epivizr
session.
app$chart_mgr$list_chart_types()
## type js_class num_settings
## 1 BlocksTrack epiviz.plugins.charts.BlocksTrack 6
## 2 LineTrack epiviz.plugins.charts.LineTrack 15
## 3 ScatterPlot epiviz.plugins.charts.ScatterPlot 11
## settings
## 1 title,marginTop,marginBottom,marginLeft,marginRight,minBlockDistance
## 2 title,marginTop,marginBottom,marginLeft,marginRight,measurementGroupsAggregator,...
## 3 title,marginTop,marginBottom,marginLeft,marginRight,measurementGroupsAggregator,...
## num_colors
## 1 10
## 2 8
## 3 10
Chart type settings, e.g., line widths, colors, margins etc. can be modified dynamically. Settings available for each chart type are briefly listed in this call. See below for further detail on customizing charts through settings and colors.
Once the browser is open we can visualize the tcga_colon_blocks
object containing blocks of methylation modifications in colon cancer. We use the plot
method to do so.
blocks_chart <- app$plot(tcga_colon_blocks, datasource_name="450k colon_blocks")
Windows users: We need the service
function to let the interactive R
session serve data requests from the browser client as you interact with the UI. Escape (using ctl-c
or esc
depending on your environment) to continue with the interactive R
session.
app$server$service()
You should now see that a new track is added to the Epiviz
web app. You can think of this track as an interactive display device in R
. As you navigate on the browser, data is requested from the R
session through the websocket connection. Remember to escape to continue with your interactive R
session if you are not running “non-blocking” mode. The blocks_chart
object inherits from class EpivizChart
, which we can use to get information about the data being displayed and the chart used to display it. Note that the “brushing” effect we implement in Epiviz
works for epivizr
tracks as well.
Now that we have interactive data connections to Epiviz
we may want to iteratively compute and visualize our data. For example, we may want to only display methylation blocks inferred at a certain statistical significance level. In this case, we will filter by block size.
# subset to those with length > 250Kbp
keep <- width(tcga_colon_blocks) > 250000
# get the data object for chart
ms_obj <- app$get_ms_object(blocks_chart)
app$update_measurements(ms_obj, tcga_colon_blocks[keep,])
Now, only this subset of blocks will be displayed in the already existing track.
To modify default chart colors or settings, we need to know what settings can be applied to the chart. Again, as these are defined dynamically in the epiviz JS app we use a method, print_chart_type_info
on the EpivizChartMgr
class to list settings that can be applied to specific chart type. For example, to list chart settings available for a BlocksTrack
chart we use the following:
app$chart_mgr$print_chart_type_info("BlocksTrack")
## Settings for chart type BlocksTrack
## id label default_value possible_values
## 1 title Title
## 2 marginTop Top margin 25
## 3 marginBottom Bottom margin 23
## 4 marginLeft Left margin 20
## 5 marginRight Right margin 10
## 6 minBlockDistance Minimum block distance 5
## type
## 1 string
## 2 number
## 3 number
## 4 number
## 5 number
## 6 number
## Colors: #1f77b4, #ff7f0e, #2ca02c, #d62728, #9467bd, #8c564b, #e377c2, #7f7f7f, #bcbd22, #17becf
There are two ways to set settings and colors to a chart: when the chart is initially created using the plot
method, or after the chart is created using the set
method in class EpivizChart
. We will illustrate both ways of doing this:
If using the plot
method, use the list_chart_settings
method as above to list settings that can be applied to a chart type. For example, for a BlocksTrack
chart, we can set minBlockDistance
as a setting for the plot function which controls how close genomic regions can be before they are merged into a single rectangle in the chart. We can also change colors used in the chart this way.
settings <- list(minBlockDistance=50)
colors <- c("#d15014", "#5e97eb", "#e81ccd")
blocks_chart <- app$plot(tcga_colon_blocks, datasource_name="450k colon_blocks", settings=settings, colors=colors)
This will create a second blocks track using a different color map.
On the other hand, if the BlockChart
is already added to the epiviz JS application session, use the set_chart_settings
method to update settings and colors.
# create a list of settings and colors to update the plot
settings <- list(minBlockDistance=100)
colors <- c("#5e97eb", "#d15014", "#e81ccd")
app$chart_mgr$set_chart_settings(blocks_chart, settings=settings, colors=colors)
This changes the color map for the second blocks track added. Use the print_chart_info
method to list settings and colors currently used in the chart.
# to list applied chart settings
app$chart_mgr$print_chart_info(blocks_chart)
## Chart settings for chart id epivizChart_2 :
## id label default_value possible_values
## 1 title Title
## 2 marginTop Top margin 25
## 3 marginBottom Bottom margin 23
## 4 marginLeft Left margin 20
## 5 marginRight Right margin 10
## 6 minBlockDistance Minimum block distance 100
## type
## 1 string
## 2 number
## 3 number
## 4 number
## 5 number
## 6 number
## Colors: #5e97eb, #d15014, #e81ccd
Methods are also available directly on the EpivizChart
objects. E.g., calls blocks_chart$set(settings=settings, colors=colors)
and blocks_chart$print_info()
yield the same result.
Charts can be printed as a pdf or png file through the epivizr session. To do so, use the print_chart
method:
app$chart_mgr$print_chart(blocks_chart, file_name="blocks_chart", file_type="pdf")
This will create a file named blocks_chart.pdf
which will be downloaded through your web browser. It will be found in the default file download location for your web browser.
There are a number of different data types available to use through epivizr
. You can see a list of R/BioC classes that are supported using ?register
. In the previous section, we used the block
data type to register a GenomicRanges
object. To visualize methylation data at base-pair resolution from the tcga_colon_curves
GenomicRanges
object, we will use the bp
type.
# add low-filter smoothed methylation estimates
means_track <- app$plot(tcga_colon_curves, datasource_name="450kMeth",type="bp",columns=c("cancerMean","normalMean"))
NOTE: You can adjust track settings to change how this new track looks like. For instance, to show all points in the window set the step
parameter to 1, and to see a smooth interpolation of the data set the interpolation
parameter to basis
:
means_track$set(settings=list(step=1, interpolation="basis"))
Notice that we added two lines in this plot, one for mean methylation in cancer and another for mean methylation in normal. The columns
argument specifies which columns in mcols(colon_curves)
will be displayed.
We can also add a track containing the smooth methylation difference estimate used to define methylation blocks.
diff_chart <- app$plot(tcga_colon_curves, datasource_name="450kMethDiff",type="bp",columns=c("smooth"),ylim=matrix(c(-.5,.5),nc=1))
We pass limits for the y axis in this case. To see other arguments supported, you can use the help framework in R ?"EpivizApp"
. As before, we can specify settings for this new track.
We can use the app connection object to list charts we have added so far, or to remove charts.
app$chart_mgr$list_charts()
## id type
## 1 epivizChart_1 epiviz.plugins.charts.BlocksTrack
## 2 epivizChart_2 epiviz.plugins.charts.BlocksTrack
## 3 epivizChart_3 epiviz.plugins.charts.LineTrack
## 4 epivizChart_4 epiviz.plugins.charts.LineTrack
## measurements connected
## 1 450k colon_blocks_1:450k colon_blocks
## 2 450k colon_blocks_2:450k colon_blocks
## 3 450kMeth_3:cancerMean,450kMeth_3:normalMean
## 4 450kMethDiff_4:smooth
app$chart_mgr$rm_chart(means_track)
app$chart_mgr$list_charts()
## id type
## 1 epivizChart_1 epiviz.plugins.charts.BlocksTrack
## 2 epivizChart_2 epiviz.plugins.charts.BlocksTrack
## 3 epivizChart_4 epiviz.plugins.charts.LineTrack
## measurements connected
## 1 450k colon_blocks_1:450k colon_blocks
## 2 450k colon_blocks_2:450k colon_blocks
## 3 450kMethDiff_4:smooth
Now we want to visualize the colon expression data in apColonData
object as an MA plot in Epiviz
. First, we add an "MA"
assay to the ExpressionSet
:
keep <- pData(apColonData)$SubType!="adenoma"
apColonData <- apColonData[,keep]
status <- pData(apColonData)$Status
Indexes <- split(seq(along=status),status)
exprMat <- exprs(apColonData)
mns <- sapply(Indexes, function(ind) rowMeans(exprMat[,ind]))
mat <- cbind(colonM=mns[,"1"]-mns[,"0"], colonA=0.5*(mns[,"1"]+mns[,"0"]))
assayDataElement(apColonData, "MA") <- mat
show(apColonData)
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 5339 features, 2 samples
## element names: MA, exprs
## protocolData: none
## phenoData
## sampleNames: GSM95473 GSM95474 ... GSM215114 (53 total)
## varLabels: filename DB_ID ... Status (7 total)
## varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## Annotation: hgu133plus2
epivizr
will use the annotation(apColonData)
annotation to determine genomic locations using the AnnotationDbi
package so that only probesets inside the current browser window are displayed.
eset_chart <- app$plot(apColonData, datasource_name="MAPlot", columns=c("colonA","colonM"), assay="MA")
## Loading required package: hgu133plus2.db
## Loading required package: AnnotationDbi
## Loading required package: org.Hs.eg.db
##
##
## 'select()' returned 1:many mapping between keys and columns
In this case, we specify which data is displayed in each axis of the scatter plot using the columns
argument. The assay
arguments indicates where data is obtained.
Epiviz
is also able to display plots of data in the form of a RangedSummarizedExperiment
object. After loading the tcga_colon_expression
dataset in the epivizrData
package, we can see that this object contains information on 239322 exons in 40 samples.
data(tcga_colon_expression)
show(tcga_colon_expression)
## class: RangedSummarizedExperiment
## dim: 12800 40
## metadata(0):
## assays(1): ''
## rownames(12800): 143686 149186 ... 149184 149185
## rowData names(2): exon_id gene_id
## colnames(40): TCGA-A6-5659-01A-01R-1653-07
## TCGA-A6-5659-11A-01R-1653-07 ... TCGA-D5-5540-01A-01R-1653-07
## TCGA-D5-5541-01A-01R-1653-07
## colData names(110): bcr_patient_barcode bcr_sample_uuid ...
## Basename fullID
The assay
slot holds a matrix of raw sequencing counts, so before we can plot a scatterplot showing expression, we must first normalize the count data. We use the geometric mean of each row as a reference sample to divide each column (sample) by, then use the median of each column as a scaling factor to divide each row (exon) by.
ref_sample <- 2 ^ rowMeans(log2(assay(tcga_colon_expression) + 1))
scaled <- (assay(tcga_colon_expression) + 1) / ref_sample
scaleFactor <- Biobase::rowMedians(t(scaled))
assay_normalized <- sweep(assay(tcga_colon_expression), 2, scaleFactor, "/")
assay(tcga_colon_expression) <- assay_normalized
Now, using the expression data in the assay
slot and the sample data in the colData
slot, we can compute mean exon expression by sample type.
status <- colData(tcga_colon_expression)$sample_type
index <- split(seq(along = status), status)
logCounts <- log2(assay(tcga_colon_expression) + 1)
means <- sapply(index, function(ind) rowMeans(logCounts[, ind]))
mat <- cbind(cancer = means[, "Primary Tumor"], normal = means[, "Solid Tissue Normal"])
Now, create a new RangedSummarizedExperiment
object with the two column matrix, and all the information about the features of interest, in this case exons, are stored in the rowRanges
slot to be queried by Epiviz
.
sumexp <- SummarizedExperiment(mat, rowRanges=rowRanges(tcga_colon_expression))
se_chart <- app$plot(sumexp, datasource_name="Mean by Sample Type", columns=c("normal", "cancer"))
Again, the columns
argument specifies what data will be displayed along which axis.
We can navigate to a location on the genome using the navigate
method of the app object:
app$navigate("chr11", 110000000, 120000000)
There is a convenience function to quickly navigate to a series of locations in succession. We can use that to browse the genome along a ranked list of regions. Let’s navigate to the 5 most up-regulated exons in the colon exon expression data.
foldChange <- mat[,"cancer"]-mat[,"normal"]
ind <- order(foldChange,decreasing=TRUE)
# bounding 1Mb around each exon
slideshowRegions <- trim(rowRanges(sumexp)[ind] + 1000000L)
app$slideshow(slideshowRegions, n=5)
To print the current epiviz workspace as a pdf or png, we use the print_workspace
method:
app$print_workspace(file_name="workspace", file_type="pdf")
This will create a file named workspace.pdf
and will be downloaded through your web browser. It will be saved to the default file download location for your web browser.
To close the connection to Epiviz
and remove all tracks added during the interactive session, we use the stop_app
function.
app$stop_app()
The epivizrStandalone
all files required to run the web app UI locally. This feature can be used to browse any genome of interest. See that packages vignette for more information.
sessionInfo()
## R version 3.3.1 (2016-06-21)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 14.04.4 LTS
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] hgu133plus2.db_3.2.3 org.Hs.eg.db_3.3.0
## [3] AnnotationDbi_1.34.4 SummarizedExperiment_1.2.3
## [5] GenomicRanges_1.24.2 GenomeInfoDb_1.8.3
## [7] IRanges_2.6.1 S4Vectors_0.10.2
## [9] antiProfilesData_1.8.0 Biobase_2.32.0
## [11] BiocGenerics_0.18.0 epivizr_2.2.4
##
## loaded via a namespace (and not attached):
## [1] Rcpp_0.12.6 BiocInstaller_1.22.3
## [3] formatR_1.4 XVector_0.12.1
## [5] GenomicFeatures_1.24.5 bitops_1.0-6
## [7] tools_3.3.1 zlibbioc_1.18.0
## [9] biomaRt_2.28.0 digest_0.6.9
## [11] RSQLite_1.0.0 evaluate_0.9
## [13] graph_1.50.0 DBI_0.4-1
## [15] OrganismDbi_1.14.1 yaml_2.1.13
## [17] epivizrData_1.0.3 rtracklayer_1.32.2
## [19] stringr_1.0.0 knitr_1.13
## [21] Biostrings_2.40.2 R6_2.1.2
## [23] XML_3.98-1.4 RBGL_1.48.1
## [25] BiocParallel_1.6.3 rmarkdown_1.0
## [27] magrittr_1.5 epivizrServer_1.0.3
## [29] Rsamtools_1.24.0 htmltools_0.3.5
## [31] GenomicAlignments_1.8.4 mime_0.5
## [33] httpuv_1.3.3 stringi_1.1.1
## [35] RCurl_1.95-4.8 rjson_0.2.15