restfulSE: A semantically rich interface for cloud-scale genomics with Bioconductor

The SummarizedExperiment class and related methods

Let Q denote a matrix of quantifications arising from a genome scale assay with G assay features measured on N experimental samples. The elements of Q are the numbers q_ij, i = 1, … , G, j = 1, …, N. Bioconductor’s SummarizedExperiment structure manages feature quantifications with associated metadata about assay features and samples.

In the 10x mouse neuron dataset, G = 27998 and N = 1.3 million. Each of the G features is a gene, and it is useful to have handy a number of feature annotations like gene name, location, functional role; suppose each gene has F such features recorded. When these quantifications and associated annotations are managed in a Bioconductor SummarizedExperiment X, the matrix Q is programmatically bound to a G × F table of feature-level metadata accessible by the rowData method, and to an N × R table of sample-level metadata accessible by colData, where R denotes the number of sample-level metadata features recorded (Huber et al.¹). See Figure 1.

Figure 1. Schematic of SummarizedExperiment class structure.

Colored regions of panels within the schematic are linked with command examples in colored text beneath the panels. For example, the purple command subsetByOverlaps(se, roi) would produce a restricted RangedSummarizedExperiment instance with features limited to those colored purple. The sizeFactors component is specific to a subclass for single cell data.

In the context of R programming, let K denote a vector of feature identifiers, S denote a vector of sample identifiers. The standard subsetting idiom X[K,S] expresses filtering of the all the information in Q and the associated metadata to features K and samples S. A GRanges instance (Lawrence et al.²) defining genomic coordinates for features may be bound to X, facilitating queries defined by genomic location (using, for example, subsetByOverlaps) to isolate features coincident with or near the elements of a set of query genomic ranges (eg., binding peaks). This outline of genomic data representation and analysis is characteristic of Bioconductor.

Examples of remote back ends

Google BigQuery. The Institute for Systems Biology Cancer Genomics Cloud project (ISB-CGC) (ISB³) uses Google BigQuery to provide access to various public cancer genomics resources including TCGA and the PanCancer Atlas (Hoadley et al.⁴). The pancan_SE function of restfulSE constructs queries that derive SummarizedExperiment instances using quantifications and annotations for PanCancer atlas experiments managed in BigQuery tables.

HDF Scalable Data Service (HSDS). An AWS S3-based distributed data object model for HDF5 datasets, including a RESTful API to structure, populate, and query HDF5 archives, has been implemented by the HDF Group. A number of datasets of interest in bioinformatics are served through HDF Kita Lab in the /shared/bioconductor folder.

Lazy data retrieval via DelayedArray

The restfulSE package provides interfaces to BigQuery and HSDS so that the numerical content housed in these services satisfies the API of the Bioconductor DelayedArray (Pagès and Hickey⁵). Any DelayedArray instance can serve as the assay component of a SummarizedExperiment instance. Thus the capacities of SummarizedExperiment to bind semantically rich metadata to genome-scale assays are extended implicitly to data resources for which no standards exist for associating substantive metadata.

In conjunction with the rhdf5client and bigrquery packages, restfulSE functions translate filtering and selection operations which are readily defined using rowData, rowRanges, colData into formal queries resolvable by the HDF5 and BigQuery services. Numerical results are transmitted from server to client only when needed.

Hybrid data/annotation strategy for integrative analysis

The following code chunk, which generates Figure 2, illustrates the use of the restfulSE protocol with the ISB-CGC BigQuery back end.

library(SummarizedExperiment)
library(BiocOncoTK)       # uses restfulSE for cancer bioinformatics
bq = pancan_BQ()          # need CGC_BILLING to authenticate
seCOAD = buildPancanSE(bq, acronym="COAD", assay="RNASeqv2")
seCOAD = bindMSI(seCOAD)  # update to include MSIsensor scores
par(mfrow=c(1,2))         # figure layout
amap = c("29126"="PD-L1", "925"="CD8A") # entrez:symbol mapping
bxs <- lapply( c("29126", "925"),       # for genes of interest
  function(x) boxplot(split(log2(as.numeric(assay( seCOAD[x,]))+1),
      seCOAD$msiTest >= 4), names = c("<4", ">=4"), ylab=amap[x],
      xlab="MSIsensor score")
  )

Our interest is in replicating part of Figure 5C of Bailey et al.⁶. In that paper, it is shown that microsatellite instability (MSI) is associated with different expression signatures of immune cell infiltration for adenocarcinomas of colon (COAD) and stomach (STAD), and uterine corpus endometrial carcinoma (UCEC). The MSI scores developed using MSIsensor are found in Table S5 of Ding et al.⁷. These scores are not available in BigQuery, but can be combined with the assay data using standard R programming, leading to a hybrid data/annotation strategy.

Figure 2. Association of MSI sensor scores with distributions of PDL-1 and CD8A in TCGA colorectal adenocarcinoma samples (COAD).

Functions in the BiocOncoTK package (Carey⁸) build on restfulSE functionality to a) authenticate the user to the BigQuery platform, b) select a tumor type (COAD) and assay for SummarizedExperiment construction, c) bind Ding et al.’s MSI values as sample-level data variable msiTest, d) acquire and transform the PD-L1 and CD8A (Entrez IDs 29126 and 925) expression values, and e) form the stratified boxplot. The basic findings of Bailey et al. are replicated. Enhancement of the code to produce a display covering more genes and tumor types is demonstrated in the BiocOncoTK package vignette. Note that in this example, expression values are only downloaded for the genes requested, without altering the end user programming paradigm of working with a SummarizedExperiment instance.

HDF Scalable Data Service

Figure 3 demonstrates use of a RESTful SummarizedExperiment, with assay data provided in the object /shared/bioconductor/darmgcls.h5 at hsdshdflab.hdfgroup.org. Briefly, as a prelude to single-cell RNA-sequencing of glioblastoma (GBM) tumors from four patients, Darmanis et al.⁹ used immunopanning to increase the proportion of non-neoplastic cells that constitute the “migrating front” of progression of glioblastoma. Antibody to CD45 was used to capture microglial cells. Figure 3 provides code to compare the distribution of CD45 expression among the classes of cells as labeled in the metadata of GSE84465, the NCBI GEO archive from which the quantifications were derived. In this example, data on one gene from all cells is retrieved when the statement defining vector vals is executed. The display can be recapitulated for other genes by substituting different symbols in the statement computing ind. The DelayedArray framework leveraged here enables basic computations of this kind without loading the entire matrix into memory.

library(rhdf5client)
library(SummarizedExperiment)
library(ggplot2)
cdar = BiocOncoTK::darmGBMcls
ind = match("PTPRC", rowData(cdar)$symbol)
var = gsub("selection: ", "",
       cdar$characteristics_ch1.8)
vals = log10(assay(cdar[ind,])+1)
ddd = data.frame(log10norm=vals, pan=var)
ggplot(ddd, aes(x=log10norm, colour=pan)) +
  geom_density() + ylim(0,1) +
  xlab("log10 CD45+1")

Figure 3. Density estimates for log10 CD45 expression in single-cell RNA-seq studies of glioblastoma.