Contents

1 Setup

1.1 Install and load package

if (!require("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("GenomicSuperSignature")
BiocManager::install("bcellViper")
library(GenomicSuperSignature)
library(bcellViper)

1.2 Download RAVmodel

You can download GenomicSuperSignature from Google Cloud bucket using GenomicSuperSignature::getModel function. Currently available models are built from top 20 PCs of 536 studies (containing 44,890 samples) containing 13,934 common genes from each of 536 study’s top 90% varying genes based on their study-level standard deviation. There are two versions of this RAVmodel annotated with different gene sets for GSEA - MSigDB C2 (C2) and three priors from PLIER package (PLIERpriors).

The demo in this vignette is based on human B-cell expression data, so we are using the PLIERpriors model annotated with blood-associated gene sets.

Note that the first interactive run of this code, you will be asked to allow R to create a cache directory. The model file will be stored there and subsequent calls to getModel will read from the cache.

RAVmodel <- getModel("PLIERpriors", load=TRUE)
#> [1] "downloading"
RAVmodel
#> class: PCAGenomicSignatures 
#> dim: 13934 4764 
#> metadata(8): cluster size ... version geneSets
#> assays(1): RAVindex
#> rownames(13934): CASKIN1 DDX3Y ... CTC-457E21.9 AC007966.1
#> rowData names(0):
#> colnames(4764): RAV1 RAV2 ... RAV4763 RAV4764
#> colData names(4): RAV studies silhouetteWidth gsea
#> trainingData(2): PCAsummary MeSH
#> trainingData names(536): DRP000987 SRP059172 ... SRP164913 SRP188526

version(RAVmodel)
#> [1] "1.1.1"

1.3 Example dataset

The human B-cell dataset (Gene Expression Omnibus series GSE2350) consists of 211 normal and tumor human B-cell phenotypes. This dataset was generated on Affymatrix HG-U95Av2 arrays and stored in an ExpressionSet object with 6,249 features x 211 samples.

data(bcellViper)
dset
#> ExpressionSet (storageMode: lockedEnvironment)
#> assayData: 6249 features, 211 samples 
#>   element names: exprs 
#> protocolData: none
#> phenoData
#>   sampleNames: GSM44075 GSM44078 ... GSM44302 (211 total)
#>   varLabels: sampleID description detailed_description
#>   varMetadata: labelDescription
#> featureData: none
#> experimentData: use 'experimentData(object)'
#> Annotation:

You can provide your own expression dataset in any of these formats: simple matrix, ExpressionSet, or SummarizedExperiment. Just make sure that genes are in a ‘symbol’ format.

2 Which RAV best represents the dataset?

validate function calculates validation score, which provides a quantitative representation of the relevance between a new dataset and RAV. RAVs that give the validation score is called validated RAV. The validation results can be displayed in different ways for more intuitive interpretation.

val_all <- validate(dset, RAVmodel)
head(val_all)
#>          score PC          sw cl_size cl_num
#> RAV1 0.2249382  2 -0.05470163       6      1
#> RAV2 0.3899341  2  0.06426256      21      2
#> RAV3 0.1887409  2 -0.01800335       4      3
#> RAV4 0.2309703  6 -0.04005584       7      4
#> RAV5 0.2017431  2  0.05786189       3      5
#> RAV6 0.1903602  8 -0.02520973       3      6

2.1 HeatmapTable

heatmapTable takes validation results as its input and displays them into a two panel table: the top panel shows the average silhouette width (avg.sw) and the bottom panel displays the validation score.

heatmapTable can display different subsets of the validation output. For example, if you specify scoreCutoff, any validation result above that score will be shown. If you specify the number (n) of top validation results through num.out, the output will be a n-columned heatmap table. You can also use the average silhouette width (swCutoff), the size of cluster (clsizecutoff), one of the top 8 PCs from the dataset (whichPC).

Here, we print out top 5 validated RAVs with average silhouette width above 0.

heatmapTable(val_all, RAVmodel, num.out = 5, swCutoff = 0)

2.2 Interactive Graph

Under the default condition, plotValidate plots validation results of all non single-element RAVs in one graph, where x-axis represents average silhouette width of the RAVs (a quality control measure of RAVs) and y-axis represents validation score. We recommend users to focus on RAVs with higher validation score and use average silhouette width as a secondary criteria.

plotValidate(val_all, interactive = FALSE)