---
title: "GeneExpressionSignature: Computing pairwise distances between different biological states"
author:
- Yang Cao
- Fei Li
package: GeneExpressionSignature
output: BiocStyle::html_document
vignette: >
%\VignetteIndexEntry{GeneExpressionSignature}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
bibliography: ref.bib
---
```{r, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(BiocStyle)
```
# Introduction
The `r Biocpkg("GeneExpressionSignature")` package utilizes gene expression
profiles to measure the similarity between different biological states. It
provides two algorithms for similarity measurement: the GSEA algorithm which is
mentioned in [@iorio2010discovery] and the PGSEA algorithm in
`r Rpackage("PGSEA")` package. A further description of the measurement methods
based on gene expression signature can be found in Lamb [@lamb2006connectivity],
Hu [@hu2009human] and Iorio [@iorio2010discovery].
This manual is a brief introduction to structure, functions and usage of
`r Biocpkg("GeneExpressionSignature")` package. It shows how the biological
similarity is determined through a series of calculation steps and how that
information can be used for further cluster analysis.
The current version of `r Biocpkg("GeneExpressionSignature")` can be used only
with data coming from the same platform, examples are on the HG-U133A platform.
# Getting Started
A complete analysis procedure accepts a set of gene expression profiles
representing different biological states as input, and generates a similarity
matrix as output. It can be divided into three steps: 1)data ranking, 2)rank
merging, and 3)similarity measuring.
First, we load the package by entering the following command in your R session:
```{r load-ges}
library(GeneExpressionSignature)
```
## Data Ranking
Gene expression profiles should be properly preprocessed before analysis as
prerequisite, including background correction, normalization and summarization.
Instead of the exact values, ranks of gene expression levels are used in the
following procedure. A ranked list of genes was obtained first by sorting the
micro-array probe-set identifiers according to the different expression
values (count or ratio). It should be noticed that there is no standard methods
for data preprocessing, and there is a function `getRLs()` which takes the
method in C-MAP for data preprocessing just for reference. We can obtain ranked
lists matrix by calling `getRLs()`.
Your experimental data could be used for analyzing, or users can download
gene-expression profiles from the GEO database with R package
`r Biocpkg("GEOquery")`. Users can see the doc in the `r Biocpkg("GEOquery")`
for more details.
As an example, we download data from GEO database with package
`r Biocpkg("GEOquery")`. Then combined the treatment expression values to form
a treatment matrix as well as the control expression values.
```{r get-ranklist,message=FALSE}
# If you have network access
#GSM118720 <- getGEO('GSM118720')
# GSM118721 <- getGEO('GSM118721')
if (require(GEOquery)){
#treatment gene-expression profiles
GSM118720 <- getGEO(
filename = system.file(
"extdata/GSM118720.soft",
package = "GeneExpressionSignature")
)
#control gene-expression profiles
GSM118721 <- getGEO(
filename=system.file(
"extdata/GSM118721.soft",
package = "GeneExpressionSignature")
)
#data ranking according to the different expression values
control <- as.matrix(as.numeric(Table(GSM118721)[, 2]))
treatment <- as.matrix(as.numeric(Table(GSM118720)[, 2]))
ranked_list <- getRLs(control, treatment)
}
```
## Rank Merging
By rank merging, multiple ranked lists are merging into a single ranked list,
referred as prototype ranked list (PRL), representing certain kind of biological
state. This procedure is mainly performed before similarity measuring, and
applied to specific situations that occur when multiple ranked list are
assigned to one single biological state with different cell types or
experimental condition.
However, two different cases should be considered: 1) all ranked list with the
same biological state are treated equally important; 2) each individual ranked
lists has its own ranked weights. This package provides two commonly employed
algorithms: one utilizes the Kruskal algorithm proposed by Iorio
[@iorio2010discovery] for the former case and another takes the average ranking
technique a simple but rather useful method. Function `RankMering()` is provided
for aggregating the ranked lists into one or many PRLs according their
phenotypic data. All the things that we need to do is construct a `ExpressionSet`
object as input, with ranked lists as assay data and corresponding biological
states as phenotypic data.
For convenience, ranking data stored as `ExpressionSet` class in `eset` object
as input data, with ranked lists (obtained by calling `getRLs()`) as assay data
and corresponding biological states as phenotypic data. As an example, we start
from loading cultured the exampleSet data, a subset of C-MAP
[@lamb2006connectivity] as sample data, which is a large reference catalogue of
gene expression data from cultured human cells perturbed with many chemicals
and genetic reagents. The sub dataset is composed of 50 paired gene expression
profiles involving 22283 genes. This profiles are obtained from cells treated
15 compounds respectively, the values of which already converted to rank orders.
```{r sample-data}
data(exampleSet)
show(exampleSet)
exprs(exampleSet)[c(1:10), c(1:3)]
levels(as(phenoData(exampleSet), "data.frame")[, 1])
```
Rank merging process will generate a mergingSet of 15 PRLs from 50 paired
expression profiles with each PRL corresponding one of 15 compounds
respectively.
```{r rank-merge}
MergingSet <- RankMerging(exampleSet, "Spearman", weighted = TRUE)
show(MergingSet)
```
## Similarity Measuring
One single combined PRL for a state was obtained after rank merging procedure.
These PRLs are used to measure the similarity of the gene signature across
different biological states by scoring functions `ScoreGSEA()` and
`ScorePGSEA()`. Not all the genes are involved in similarity measuring, but
only a subset of genes called gene signature whose combined expression pattern
is uniquely characteristic of the biological state. Generally the genes used as
gene signatures in the similarity scoring procedure are predefined by priori
knowledge. Iorio [@iorio2010discovery] proposed an "optimal signature"
approach by taking the most up-regulated genes and the most down-regulated
genes as gene signature.The size of gene signatures need to be considered, which
is taken as another parameter besides the PRLs in similarity measuring. In most
cases, the default size of gene signature is 250 for genome-wide expression
profile.
Suppose `N` is the number of PRLs (also same as the number of biological states),
an `N x N` distance matrix is generated by similarity measurement. For mergingSet,
we will get a `15 x 15` matrix corresponding to the similarity distances
between these compounds.
```{r gsea}
ds <- ScoreGSEA(MergingSet, 250, "avg")
ds[1:5, 1:5]
```
## Signature Distance
As we mentioned above, four algorithms implemented as functions `getRLs()`,
`RankMerging()`, `ScoreGSEA()`, and `ScorePGSEA()`, one is for data
preprocessing, one called Iorio algorithm is for rank merging, the other two
algorithms called GSEA and PGSEA are for similarity measuring. Moreover,
function `SignatureDistance()` is provided to serve as a single entry and easy
access point to rank merging and similarity measuring, which runs through the
including rank merging and scoring, and is recommended to use in most cases.
Data ranking is not integration into this function for no standard methods for
data preprocessing and gene-expression data types is uncertain. Furthermore,
there is no effective method to integrate data from different platforms.
Function `getRLs()` which takes the method in C-MAP for data preprocessing just
for reference.
```{r sig-dis}
SignatureDistance(
exampleSet,
SignatureLength = 250,
MergingDistance = "Spearman",
ScoringMethod = "GSEA",
ScoringDistance = "avg",
weighted = TRUE
)
```
# Implementation Details
## Adaptively Weighted Rank Merging
The Iorio's rank merging algorithm utilizes Kruskal algorithm
[@cormen1990introduction] to merge the ranked lists which corresponding to a
same biological state. The distance of these ranked lists must be calculated
first, a measure of the distance between two ranked lists is computed using
"Spearman" algorithm or "Kendall tau" algorithm. It should be noticed that
rank merging with Kendall tau distance is time consuming, so we recommend
selecting the "Spearman" distance. Next, merge the two or more ranked lists
with the same biological state using "Borda merging" algorithm.
According to the "Kruskal" algorithm method [@cormen1990introduction], this
rank merging algorithm searches for the two ranked lists with the smallest
Spearman's Footrule distance first, and then merges them using the Borda
Merging method, obtaining a new ranked list. Finally, the new list replaces the
two unmerged lists. This process won't terminate until only one list remains.
For convenience, users can directly obtain a PRL for each state by the function
`RankMerging()`, which uses "Sprearman", "BordaMerging", and "Kruskal"
algorithms to aggregate the ranked lists obtained with the same biological
state. For instance, we will merge the sample data which with 50 samples into
15 samples.
```{r merge-detail}
MergingSet <- RankMerging(exampleSet, "Spearman", weighted = TRUE)
show(MergingSet)
```
### Eqully Weighted Rank Merging
A simple but rather useful method for this problem is the average ranking
technique. The technique is a two step process when we are under the assumption
that importance is equally weighted for each ranked list. First step is
to calculate average rank for each ranked list and then the second step is to
construct their final rankings.
### Similarity Measuring
Once ranked lists with same biological states are merged to one single PRL,
Gene Set Enrichment Analysis (GSEA) and Parametric Gene Set Enrichment Analysis
(PGSEA) are adopted to measure the similarity among these PRLs.
GSEA algorithm [@subramanian2005gene] is a nonparametric, rank-based method for
similarity measuring to determine whether a priori defined set of genes shows
statistically significant, concordant differences between two biological states,
whereas PGSEA algorithm is a modified gene set enrichment analysis method based
on a parametric statistical analysis model. Both of these two functions gives
the corresponding p value, function `ScoreGSEA()` calcutes the empirical p
values from Monte Carlo Procedures [@north2002note].
```{r gsea-detail}
ds <- ScoreGSEA(MergingSet,250,"avg")
ds[1:5,1:5]
```
```{r pgsea-detail}
ds <- ScorePGSEA(MergingSet,250,"avg")
ds[1:5,1:5]
```
# Futher Analysis
To illustrate how to use `r Biocpkg("GeneExpressionSignature")` in analysis of
gene expression signatures, affinity propagation clustering can be used to group
these biological states by the similarity of gene signature. Affinity
propagation cluster algorithm iteratively searches for optimal clustering by
maximizing an objective function called net similarity. Here, we use function
in `r CRANpkg("apcluster")` package to classify the 15 biological states into
3 groups. In this step, R package `r CRANpkg("apcluster")` should also be
installed on your computer.
```{r cluster}
if (require(apcluster)){
library(apcluster)
clusterResult <- apcluster(1 - ds)
show(clusterResult)
}
```
Cytoscape is used to visualize the result of clustering. In the network, nodes
denotes different compounds (cell states treated with different compounds), and
the edge means the similarity distance between these two compounds is lower
than a threshold, which is 0.68 here. Different colors denote different groups,
as the classification of compounds. We note that the largest group is numbered
9 nodes, and the other two consist of 3 nodes for each group.
```{r cluster-graph,echo=FALSE}
knitr::include_graphics("cluster.png")
```
# Session Information {-}
```{r session}
sessionInfo()
```
# References