--- title: STATegRa User's Guide author: The STATegra Consortium date: "2018-04-04" bibliography: bibliography.bib output: BiocStyle::html_document: toc: true fig_retina: false vignette: > %% \VignetteEngine{knitr::rmarkdown} %% \VignetteIndexEntry{STATegRa User's Guide} --- ```{r, echo=F, message=F, results="asis"} # this block is invisible BiocStyle::markdown() require(STATegRa) require(Biobase) require(gridExtra) require(ggplot2) g_legend<-function(a.gplot){ tmp <- ggplot_gtable(ggplot_build(a.gplot)) leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box") legend <- tmp$grobs[[leg]] return(legend)} ``` # Introduction Recent developments in high-throughput technologies for studying biological systems enables the researcher to simultaneously obtain several different types of data ("omics") over the course of an experiment. There exist many techniques for analysing the behaviour of these omics individually, but combining multiple classes of omics data can be used to give a better understanding of the biological system in question; the whole is greater than the sum of its parts. Integration of different types of omics data is an increasingly important technique for studying biological systems. The first step in this kind of integration analysis is to identify patterns in data shared by all the omics classes, and use these patterns to identify outliers. The most common techniques for this sort of analysis are clustering and principal components analysis. The `STATegRa` package provides several different techniques for the evaluation of reproducibility among samples and across experimental conditions by combining the information contained in multiple omics datasets. This is intended as a starting point for further integration analysis of any multi-omics dataset. The `STATegRa` package implements two main utilities for this purpose: component analysis and clustering. [Component Analysis](#omics-component-analysis) ~ Three different techniques for analysing the common and distinctive variability between two different, multi-omics datasets are provided, along with various utility functions and plotting tools to evaluate the results. [Clustering](#omics-clustering) ~ Methods are provided to cluster together features across different omics types, with a view to finding interesting similarities between features (rather than similarities between samples). Furthermore, the next important step is the identification of genes, which are differential expressed between the experimental conditions under study, according to the different data types considered as a whole. The `STATegRa` package provides the possibility to identify genes which are differential expressed in one or more data types. One main utility is implemented for this purpose: omicsNPC. [omicsNPC](#omicsNPC) ~ omicsNPC is specifically devised for identifying differentially expressed genes by ‘holistically’ combining different omics data. It applies the NonParametric Combination (NPC) methodology [@Pesarin2010] in order to combine all available information by taking into account possible between-datasets correlations. This guide provides an overview of the different techniques included in the package, some worked examples of using the tools and some guidance on interpretation of the results obtained. # Getting Started The `r Biocpkg("STATegRa")` package can be obtained from the [Bioconductor repository](http://www.bioconductor.org/). Load the `STATegRa` package into an `R` session by typing: ```{r, eval=F} library(STATegRa) # Load STATegRa package ``` General information about usage of the package and the algorithms used can be found in the package vignette. In addition, every public function in the package is documented and help can be found in the normal R fashion: ```{r, eval=F} help(package="STATegRa") ## Package help ?omicsCompAnalysis ## Specific function help ``` # Omics Component Analysis ## Overview The joint analysis of multiple omic datasets, both containing different classes of data and from different experimental conditions, could provide a "global" view on the biological system of interest. The major challenge in this type of analysis is to distinguish between the underlying mechanisms affecting all datasets, and the particular mechanisms which affect each omic dataset separately. Three different methods are provided to this end: *DISCO-SCA* [@van2012disco; @schouteden2013sca; @schouteden2014performing], *JIVE* [@lock2013joint] and *O2PLS* [@trygg2003o2]. Each method provides the user with a decomposition of the variability of the composite data into common and distinctive variability. All of them are based on singular value decomposition (SVD) of the data matrix, however they use different models to accomplish this. The *DISCO-SCA* [@van2012disco] approach consists of two steps. First a Simultaneous Components Analysis (SCA) is performed, then the scores obtained are rotated into a *DIS*tinctive and *CO*mmon structure (hence, *DISCO*). Therefore, by applying SCA approach, each block of data $X_k$ of size $I\times J_k$ becomes $$X_k=TP_k^T+E_k$$ with $T$ the $I\times R$ matrix of components scores that is shared between all blocks and $P_k$ the $J_k\times R$ matrix of components loadings for block $k$. Then, a rotation criterion is used where the target is the rotation which specifies distinctive components as components having zero scores in the positions that correspond to the data blocks the component does not underlie, and the remaining entries are arbitrary. The rotation matrix $B$ is found by minimizing $min(B)||W \circ (P_{target}-[P_1^TP_2^T]B)||^2 \mbox{ such that } B^TB=I=BB^T$, where $W$ is a binary matrix having ones in the positions of the entries in the target and zero elsewhere. The *JIVE* approach [@lock2013joint] model is as following: Let $X_1,X_2$ be two blocks of data and $X=[X_1,X_2]$ represent the joint data, then the JIVE decomposition is defined as: $$X_i=J_i+A_i+\epsilon_i \mbox{, }i=1,2$$ where $J=[J_1,J_2]$ is the $p\times n$ matrix of rank $r] data-set. The dataset used for this section was obtained from the dataset described in [@van2012disco] and is available from TCGA processed data. We extracted the [classification](https://tcga-data.nci.nih.gov/docs/publications/gbm_exp/TCGA_unified_CORE_ClaNC840.txt) and the [unified gene expression](https://tcga-data.nci.nih.gov/docs/publications/gbm_exp/unifiedScaledFiltered.txt) from TCGA. The miRNA was downloaded from TCGA directly. The full dataset can be loaded by typing: ```{r} data("STATegRa_S1") ls() ``` `Block1` includes mRNA data and `Block2` includes miRNA data. ### Computing the distance between genes by using mRNA data: the bioDistclass class Firstly, we generate an `ExpressionSet` object for both the miRNA and mRNA data. ```{r} # Block1 - Expression data mRNA.ds <- createOmicsExpressionSet(Data=Block1, pData=ed, pDataDescr=c("classname")) # Block2 - miRNA expression data miRNA.ds <- createOmicsExpressionSet(Data=Block2, pData=ed, pDataDescr=c("classname")) ``` Secondly, we compute the distance between all genes in `Block1` (mRNA data) using Spearman correlation. ```{r} # Create Gene-gene distance computed through mRNA data bioDistmRNA <- bioDistclass(name="mRNAbymRNA", distance=cor(t(exprs(mRNA.ds)), method="spearman"), map.name="id", map.metadata=list(), params=list()) ``` The `bioDistmRNA` object, generated with the `bioDistclass` function, is a `bioDistclass` object that contains both the original data and the computed distance between features. ### Loading the map between miRNA and genes: the bioMap class In this section we load and store the map between miRNA and mRNA. Data file (`STATegRa_S2`) contains, as a processed matrix, the information available from TargetScan [@targetscan], which provided a set of miRNA target predictions for humans. ```{r} data(STATegRa_S2) ls() ``` This data is stored in a `bioMap` class object generated through the `bioMap` function as follows: ```{r} MAP.SYMBOL <- bioMap(name = "Symbol-miRNA", metadata = list(type_v1="Gene", type_v2="miRNA", source_database="targetscan.Hs.eg.db", data_extraction="July2014"), map=mapdata) ``` ### miRNA-Surrogate gene Distances: the bioDist function The `bioDist` function returns a `bioDistclass` object. The input is a reference feature list (genes in this example), surrogate data (miRNA, in `Block2`) and the bioMap object between reference and surrogate features. ```{r} bioDistmiRNA <- bioDist(referenceFeatures=rownames(Block1), reference="Var1", mapping=MAP.SYMBOL, surrogateData=miRNA.ds, referenceData=mRNA.ds, maxitems=2, selectionRule="sd", expfac=NULL, aggregation="sum", distance="spearman", noMappingDist=0, filtering=NULL, name="mRNAbymiRNA") ``` ### Computing weighted distances: the bioDistW function Having `bioDistmiRNA` and `bioDistmRNA` `bioDistclass` objects containing distances between genes, we aim to use weighted combinations of them to compute an single distance matrix. First we make a list of `bioDistclass` objects: ```{r} bioDistList <- list(bioDistmRNA, bioDistmiRNA) ``` Secondly we make a matrix listing containing the weighted combinations to be generated. Each row is interpreted as a combination to generate, with the elements of the row interpreted as the weight for each of the input omics. ```{r} sample.weights <- matrix(0, 4, 2) sample.weights[, 1] <- c(0, 0.33, 0.67, 1) sample.weights[, 2] <- c(1, 0.67, 0.33, 0) sample.weights ``` This matrix corresponds to generating four combinations, with the first consisting of $0\times mRNA + 1\times miRNA$ and so on. Finally, the `bioDistW` function computes the weighted combinations in the weights matrix and stores it into a `bioDistWclass` list. ```{r} bioDistWList <- bioDistW(referenceFeatures=rownames(Block1), bioDistList=bioDistList, weights=sample.weights) length(bioDistWList) ``` ## Plots ### Plotting the feature distance of each weighted combination Each `bioDistWclass` object contains a distance matrix computed through a weighted combination of distances derived from different omics. By considering the distances between these distance matrices we can project in two dimensions using Multi-Dimensional Scaling. By this approach we can visualize the effect of the different weights on the feature-to-feature distance structure. To generate such a plot: ```{r, warning=F} bioDistWPlot(referenceFeatures=rownames(Block1), listDistW=bioDistWList, method.cor="spearman") ``` ### Plotting associated features The purpose of this analysis is to generate an overall distance measure between features, so it follows that given a feature of interest we will want to find other features that are near it. For this example we use the gene IDH1, which was shown to be of relevance in the original data analysis. In order to find all other genes for which at least one weighted combination has a correlation greater than $0.7$, we do: ```{r} IDH1.F <- bioDistFeature(Feature="IDH1", listDistW=bioDistWList, threshold.cor=0.7) ``` The `bioDistFeature` function generates a matrix of associated genes (columns) depending on weighted combinations (rows); rows are named by the `bioDistWclass`'s name slot. The `IDH1.F` matrix can be plotted with the `bioDistFeaturePlot` function, as shown below. (This function is a wrapper around `heatmap.2` from `r CRANpkg("gplots")` with appropriate options, but `IDH1.F` is a normal R matrix and can be used with most other matrix-plotting tools). ```{r, message=F} bioDistFeaturePlot(data=IDH1.F) ``` ### OmicsClustering Requirements The requirements for running OmicsClustering are minimal. The data considered has to be compatible with the distance measure selected. The relevant aspect is that the mapping between features needs to be `informative enough`. We do not consider the use of OmicsClustering when the mapping between features involves less than 15-25% of the reference set of features; this number was obtained from preliminary analysis over few data sets however further investigation is being conducted. # omicsNPC ```{r, echo=F, message=F} # clear the environment before starting the omicsNPC example rm(list=ls()) ``` ## Overview ### The Problem Recent advancements in omics technologies allow to measure several data modalities, e.g., RNA-sequencing, protein levels, etc., on the same biological samples. In such settings an interesting task is the identification of genes that, according to all modalities considered as a whole, are either deregulated or associated to an outcome of interest. ### The NPC Approach *omicsNPC* implements the Non Parametric Combination (NPC) methodology [@Pesarin2010] in a way that is specifically tailored for the idiosyncrasies of omics data. First, each data type is analyzed independently. Currently, *omicsNPC* uses the package `r Biocpkg("limma")` for computing deregulation / association statistics and p-values; count data, e.g., RNAseq, are first transformed using the voom function. The user can also specify custom functions for computing relevant statistics in each dataset. The resulting p-values are combined by employing appropriate combining functions. The Tippett combining function returns findings which are supported by at least one omics modality. The Liptak function returns findings which are supportd by most modalities. The Fisher function has an intermediate behavior between those of Tippett and Liptak. Several important features make the use of omicsNPC appealing. + *omicsNPC* makes **minimal assumptions**: as permutation is employed throughout the process, no parametric form is assumed for the null distribution of the statistical tests, and the main requirement is samples to be freely exchangeable under the null-hypothesis. + It provides **global p-values** for assessing the overall contribution of all data modalities. + It is characterized by **great flexibility**: it can address study designs with multiple factors, as well as partially overlapping measurements / samples. + It frees the researcher from the need of defining and **modeling between-datasets dependencies** [@Pesarin2010]. A worked example with code is given below. ## Usage ### Loading the data For this example we use an existing TCGA^[The Cancer Genome Atlas, ] data-set. We downloaded sixteen tumour samples and the sixteen matching normal, for Breast invasive carcinoma, BRCA, batch 93. We used three types of data modalities, “RNAseq”, “RNAseqV2” and “Expression-Gene”. `RNAseq` corresponds to RNA sequensing data, *IlluminaHiSeq-RNASeq* platform, `RNAseqV2` corresponds to RNA sequencing data, *IlluminaHiSeq-RNASeqV2* platform and `Expression-Gene` corresponds to array based expression data, *AgilentG4502A-07-3* platform. For `RNAseq` and `RNAseqV2` data we keep the "raw counts". For each data type, we pooled all data to one matrix, where rows corresponded to genes and columns to samples. We selected 100 genes to be used in this example. Finally, each matrix was converted to an `ExpressionSet` object and the resulting objects were saved in a `list`. The datasets can be loaded by typing: ```{r} data("TCGA_BRCA_Batch_93") ls() ``` `TCGA_BRCA_Data` is a list, which includes three objects of the `ExpressionSet` class. The matrices containing the expression values can be visualized by typing: ```{r, results='hide'} exprs(TCGA_BRCA_Data$RNAseq) # displays the RNAseq data exprs(TCGA_BRCA_Data$RNAseqV2) # displays the RNAseqV2 data exprs(TCGA_BRCA_Data$Microarray) # displays the Exp-Gene data ``` Each row corresponds to a gene and each column to a sample. Samples are subdivided between "tumour" and "normal". The class of each sample can be visualized by typing ```{r, results='hide'} pData(TCGA_BRCA_Data$RNAseq) # class of RNAseq samples pData(TCGA_BRCA_Data$RNAseqV2) # class of RNAseqV2 samples pData(TCGA_BRCA_Data$Microarray) # class of Exp-Gene samples ``` Please note that each dataset could have different design matrices; *omicsNPC* only assumes that the last column of each design matrix contains the factor whose association with each gene should be assessed. This factor can be either continuous (e.g., IG50 values), or categorical (tumour / normal) as in the present example. ### Setting the dataTypes variable The `dataTypes` vector specifies the type of data in input and consequently the method for computing the association between each measurement and the factor of interest. If given as a character vector, possible values for `dataTypes` are `'count'` and `'continuous'`. The first value requires the data to be transormed with the voom function (from `limma`) before being analyzed, while the latter directly apply limma moderated t-statistics. Alternatively, the user can specify a list of custom functions. ```{r} dataTypes <- c("count", "count", "continuous") ``` ### Setting the combMethods variable The `combMethods` variable is a character vector with possible values: `'Fisher', 'Liptak', 'Tippett'`. All combining functions can be used simultaneously. ```{r} combMethods <- c("Fisher", "Liptak", "Tippett") ``` ### Setting the numPerms, numCores and verbose variables `numPerms` is the number of permutations to perform. The number of permutations can be adjusted depending on the number of available samples. In general 1000 permutation is usually a good point to start. ```{r} numPerms <- 1000 ``` *omicsNPC* can perform the desired permutations in parallel, and `numCores` specifies the number of cores to be used. By default, one single core is used. ```{r} numCores <- 1 ``` `verbose` is a logical. If it is set to `TRUE`, *omicsNPC* prints out the steps that it performs. ```{r} verbose <- TRUE ``` ### Run omicsNPC analysis. One main function is available for running *omicsNPC* analysis: ```{r} results <- omicsNPC(dataInput=TCGA_BRCA_Data, dataTypes=dataTypes, combMethods=combMethods, numPerms=numPerms, numCores=numCores, verbose=verbose) ``` `results` is a list containing `pvalues0` and `pvaluesNPC`. + `pvalues0` is a matrix having one column for each data modality. Each column contains the p-values assessing the association of each measurement with the factor of interest corrected for the other factors present in the study design (pheno data). + `pvaluesNPC` is a matrix having one column for each combination function. Each column provides the global p-values obtained by combining single-datasets p-values with the corresponding combination function. # References