```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r message=FALSE, warning=FALSE, include=FALSE, paged.print=FALSE}
library(kableExtra)
library(deconica)
```
<div style= "float:left;position: relative; top: -3px; left: -2px;">
```{r echo=FALSE, out.height='192.1875', out.width='210.4688'}
knitr::include_graphics("./figures-ext/deconicalogo.png")
```
</div>
This is an introduction to the [`deconICA`](https://github.com/UrszulaCzerwinska/DeconICA) R package. 

**DeconICA** stands for **Deconvolution of transcriptome through Immune Component Analysis**. 

The aim of the project is to adapt blind source separation techniques to extract immune-related signals from mixed biological samples. A great example of mixed biological sample is transcriptome measured in heterogenous tissue such as blood or tumor biopsy. 

In this vignette we present short introduction to the blind source spearation techniques, the biological foundation of the problem and finally we walk you through examples on how to use [`deconICA`](https://github.com/UrszulaCzerwinska/DeconICA) R package.

If you are interested only in practical examples of [`deconICA`](https://github.com/UrszulaCzerwinska/DeconICA), skip directly to [Tutorial](#Tut) section.

You can access this documentation on the [DeconICA website](https://urszulaczerwinska.github.io/DeconICA/).

# Background

## Blind source separation

Blind source separation (BSS) is the separation of a set of source signals from a set of mixed signals, without the aid of information (or with very little information) about the source signals or the mixing process. The separation is possible under a variety of conditions. 

A known example of BSS is a **cocktail party problem**, there is a group of people talking at the same time. You have multiple microphones picking up mixed signals, but you want to isolate the speech of a single person. BSS can be used to separate the individual sources by using mixed signals [@Hyvarinen2000]

## Application of BSS to biological data

Through decomposition of the transcriptome matrix into components (aka factors or sources) we hope to recover underlying biological functions and cell types. 

In tumor biopsies is it expected to find a part of Tumor Microenvironment (TME). TME includes tumor cells, fibroblasts, and a diversity of immune cells. Most studies have
focused on individual cell types in model tumor systems, and/or on individual molecules
mediating a crosstalk between two cells. Unraveling the complexity, organization, and mutual
interactions of TME cellular components represents a major challenge.

Several methods have been proposed to estimate the mixing proportions of sources in biological mixtures, such as: least squares regression [@Abbas2009] and more recently, non-negative least squares regression [Qiao2012], quadratic
programming [Gong2011] and supported vector
regression [@Newman2015]. Even though [@Vallania2017] shows that the used algorithm do not impact substantially the results. According to @Vallania2017, what matters are the gene signatures used as an input of aforementioned methods.

**BSS methods do not use pre-defined cell-type signatures**. The transcriptomic matrix is decomposed into a certain number of sources and then the sources are interpreted with available knowledge (gene signatures, cell profiles).

The main argument of using BSS over supervised decomposition techniques is that the obtained result is **unbiased** by *a priori* biological hypothesis (however there are always statistical hypothesis about the nature of data) or knowledge. In addition BSS techniques allow **discovery** of new biological signatures that can extend our available knowledge.

In the case of cell type separation from mix of tumor bulk, supervised techniques as CIBERSORT [@Newman2015], MCP counter [@Becht2016, @MCPcounter], TIMER [@Li2016] etc. are based on optimized blood signatures. With an evidence brought by single cell data, these signatures are not always characterizing immune cells infiltrating tumors [@Schelker2017]. Some methods, like EPIC [@Racle2017], use single-cell based signatures. However, today, the single cell based signatures are limited to few cancer subtypes and often based on small number of patients, incomparable with the heterogeneity that is hidden in the bulk transcriptome cohort studies. 

Therefore, obtaining informative cell-type signature of immune cells infiltrating tumor biopsy samples at high thoughput remains an open question that we attempt to approach with `deconICA` pipeline.

Here is a short overview of BSS or related algorithms that one can potentially use as an input to `deconICA`. At its actual state `deconICA` facilitates starting pipeline with ICA. 

```{r echo = FALSE, fig.align='center', fig.cap='**Graphical representation of dimension reduction & BSS methods.** PCA, ICA, NMF inspired by figures of Andrei Zinovyev, Convex hull: CC BY [@Wang2016]', out.width='976px', out.height='267.5px'}
knitr::include_graphics("./figures-ext/methodsBSS.png")
```

### Independent Components Analysis

Independent Component Analysis (ICA) is a matrix factorization method for data dimension reduction [@Hyvarinen2000]. ICA defines a new coordinate system in the multi-dimensional space such that the distributions of the data point projections on the new axes become as mutually independent as possible. To achieve this, the standard approach is maximizing the non-gaussianity of the data point projection distributions [@Hyvarinen2000]. There is no constraint imposed on the non-negativity (in contrary to NMF) or orthogonality (in contrast to PCA). In our analysis, the negative projections are interpreted in terms of absolute values and only one side of a component is taken into account.

A mathematical way to formalize ICA is the set of equations:

the set of individual source signals $s(t) = (s_{1}(t), ...,  s_{n}(t))^T$ is mixed using a matrix $A = [a_{ij}] \in \mathbb{R}^{m \times n}$ to produce a set of mixed signals, $x(t) = (x_{1}(t), ...,  x_{m}(t))^T$, as follows:

 \[x(t) = A \times s(t)\]

The above equation is effectively 'inverted' as follows. Blind source separation separates the set
of mixed signals $x(t)$, through the determination of an 'unmixing' matrix  to $B = [b_{ij}] \in \mathbb{R}^{m \times n}$
'recover' an approximation of the original signals, $y(t) = (y_{1}(t), ...,  y_{n}(t))^T$.

 \[y(t) = B \times x(t)\]

This algorithm uses higher-order moments for matrix approximation, considering all Gaussian
signals as noise.

Most efficient application of ICA is fastICA [@Hyvarinen2000]. However, the speed comes with a price, the results of the algorithms are not exact. This is why we recommend use of ICA with stabilization (ICASSO [@Himberg2003]) for reproducible results. More about this is the vignette [Running fastICA with icasso stabilisation](Icasso.html). 

For applications in molecular biology, Independent Component Analysis (ICA) models gene expression data as an action of a set of statistically independent hidden factors. 

Here is a small list of ICA application to biological data:

* Independent component analysis uncovers the landscape of the bladder tumor transcriptome and reveals insights into luminal and basal subtypes [@Biton2014]
* Elucidating the altered transcriptional programs in breast cancer using independent component analysis [@Teschendorff2007]
* Principal Manifolds for Data Visualization and Dimension Reduction [@Gorban2008]
* Independent component analysis of microarray data in the study of endometrial cancer [@Saidi2004]
* Blind source separation methods for deconvolution of complex signals in cancer biology [@Zinovyev2013]
* Determining the optimal number of independent components for reproducible transcriptomic data analysis [@Kairov2017]
* Application of Independent Component Analysis to Tumor Transcriptomes Reveals Specific And Reproducible Immune-related Signals [@Czerwinska2018]

### NMF

Non-negative matrix factorization (NMF) is matrix factorization technique assuming that the mixing, source and mixed matrices are all non negative. It is usually written as

\[V = WH\]

Matrix multiplication can be implemented as computing the column vectors of $V$ as linear combinations of the column vectors in $W$ using coefficients supplied by columns of $H$. That is, each column of $V$ can be computed as follows:

\[{v}_{i}={W}{h}_{i}\]

where $v_{i}$ is the $i$-th column vector of the product matrix $V$ and $h_{i}$ is the $i$-th column vector of the matrix $H$.

There are many types of NMF, that differ in implementation, ways the error terms are defined and minimized.

@CellMix1 proposed a framework of semi-supervised NMF for solving cell and tissue mixtures in his R package *CellMix* [@CellMix3]. His work deconvolution was of great inspiration for us and food for thoughts on the advantages and limits of BBS applied for cell type deconvolution.

Renaud Gaujoux is also an author of NMF R package [@NMF1; @NMF2; @NMF3]. 

Work of Cantini et al. (@Cantini2018) showed that ICA produces more reproducible results than NMF when applied to tumor transcriptomes. 

However, we can always imagine the cases where NMF factorisation will be judged more adequate than ICA. Our pipeline is adaptable to interpretation of NMF factors without a need for major adjustments. We will be seen extend vignettes to show application of deconICA for interpretation of NMF components.

Here is a small list of NMF application to biological data:

* Metagenes and molecular pattern discovery using matrix factorization [@Brunet2004]
* Tumor Clustering Using Nonnegative Matrix Factorization With Gene Selection [@Chun-HouZheng2009]
* Semi-supervised Nonnegative Matrix Factorization for gene expression deconvolution: a case study [@CellMix1]
* Post‐modified non‐negative matrix factorization for deconvoluting the gene expression profiles of specific cell types from heterogeneous clinical samples based on RNA‐sequencing data [@Liu2017]
* Virtual microdissection identifies distinct tumor- and stroma-specific subtypes of pancreatic ductal adenocarcinoma [@Moffitt2015]
* A non-negative matrix factorization method for detecting modules in heterogeneous omics multi-modal data [@Yang2015]

### Convex hull methods

An emerging family of BSS methods are convex geometry (CG)-based methods. Here, the "sources" are found by searching the facets of the convex hull spanned by the mapped observations solving a classical convex optimization problem [@Yang2015]. The convex hull-based method does not require the independence assumption, nor the non-correlation assumption which can be interesting in the setup of closely related cell types. @Wang2016 apply their method of convex analysis of mixtures (CAM) to tissue and cell mixtures claiming to provide new signatures. So far the published R-Java package does not allow to extract those signtures and it is not scalable to tumor transcriptomes. Another tool CellDistinguisher [@Newberg2018] provides an [user-friendly R package](https://github.com/GeneralElectric/CellDistinguisher). However, authors do not provide any method for estimation of number of sources. Additionally, quantitative weights are provided only for signature genes which number can vary for different sources. They do not apply their algorithm to complex mixtures as tumor transcriptome.

However, combining convex hull methods and `deconICA`can possibly lead to a meaningful interpretation.

Here is a small list of convex-hull application to biological data:

* Computational de novo discovery of distinguishing genes for biological processes and cell types in complex tissues [Newberg2018]
* Mathematical modelling of transcriptional heterogeneity identifies novel markers and subpopulations in complex tissues [@Wang2016]
* Geometry of the Gene Expression Space of Individual Cells [@Yang2015]
* Applying unmixing to gene expression data for tumor phylogeny inference [@Schwartz2010]
* Inferring biological tasks using Pareto analysis of high-dimensional data [@Hart2015]

### Attractor metagenes

Another way of generating signatures, that can be run in semi-supervised or unsupervised mode is attractor metagenes method proposed by @Cheng2013. Authors describe their rationale as follows:

>We can first define a consensus metagene from the average expression levels of all genes in the cluster, and rank all the individual genes in terms of their association (defined numerically by some form of correlation) with that metagene. We can then replace the member genes of the cluster with an equal number of the top-ranked genes. Some of the original genes may naturally remain as members of the cluster, but some may be replaced, as this process will "attract" some other genes that are more strongly correlated with the cluster. We can now define a new metagene defined by the average expression levels of the genes in the newly defined cluster, and re-rank all the individual genes in terms of their association with that new metagene; and so on. It is intuitively reasonable to expect that this iterative process will eventually converge to a cluster that contains precisely the genes that are most associated with the metagene of the same cluster, so that any other individual genes will be less strongly associated with the metagene. We can think of this particular cluster defined by the convergence of this iterative process as an "attractor" i.e., a module of co-expressed genes to which many other gene sets with close but not identical membership will converge using the same computational methodology.

The produced signatures' weights are non-negative. In the original paper, the generation of tumor signatures leads to three reproducible signatures among different tumor types. Typically with the essential parmeter $\alpha = 5$, they discovered typically approximately 50 to 150 resulting attractors. Although, it is possible by tuninig $\alpha$ obtain more or less signatures that would be possibly interpretable with `deconICA`. 

Attractor metagenes R code is available on [Synapse portal](https://www.synapse.org/#!Synapse:syn1446295).

Literature:

* Biomolecular Events in Cancer Revealed by Attractor Metagenes [@Cheng2013]
* Discovering Genome-Wide Tag SNPs Based on the Mutual Information of the Variants [@Elmas2016]
* Meta-analysis of the global gene expression profile of triple-negative breast cancer identifies genes for the prognostication and treatment of aggressive breast cancer [@Al-Ejeh2014]


# Tutorial {#Tut}

## How to install
You can install `deconICA` from GitHub with:
```{r, eval = FALSE, echo = TRUE}
#install.packages("devtools")
devtools::install_github("UrszulaCzerwinska/DeconICA", build_vignettes = TRUE)
```
or 
```{r, eval = FALSE, echo = TRUE}
install.packages("githubinstall")
githubinstall::githubinstall("DeconICA", build_vignettes = TRUE)
```

[**TO DO**] You can install the stable version from [CRAN](https://cran.r-project.org/web/packages/deconica/index.html)
```{r, eval = FALSE, echo = TRUE}
install.packages('deconica', dependencies = TRUE)
```

Then load package with
```{r}
library(deconica)
```

## Demonstration of DeconICA package {#cases}

### Simulated data

At first, we asses an ability to estimate abundance of cell types in a synthetic cell mixture. Here we use function `simulate_gene_expresssion()`inspired by `CellMix::rmix` function. However, compared to `rmix` function cell profiles distribution follow user defined distribution (uniform in `rmix`) set here by default to negative binomial which approaches the biological reality.

First we create the *mix1* that is a mix of 10 cell types mixed at random proportions. Obatained matrix has 10000 genes and 130 samples. Each cell type has 20 specific markers with 2-fold difference with respect to other genes.
```{r}
set.seed(123)
mix1<-simulate_gene_expresssion(10, 10000, 130,0, markers=20)
```

We can visualize the mixed expression matrix.

```{r}
pheatmap::pheatmap(mix1$expression)
```
We can also visualize the *purified gene expression* of each cell.
```{r}
pheatmap::pheatmap(mix1$basis_matrix)
```

Then we apply ICA (matlab version with stabilisation see [vignette: Running fastICA with icasso stabilisation](Icasso.html)) and we decompose to 11 components.
If you don't have file you can find in `data-vignettes` the file `mix1_ica.RData`.
```{r, eval=FALSE}
mix1_ica <- run_fastica (
  mix1$expression,
  overdecompose = FALSE,
  with.names = FALSE,
  gene.names = row.names(mix1$expression),
  samples = colnames(mix1$expression),
  n.comp = 11,
  R = FALSE
)
```
```{r echo=FALSE}
load(paste0(path.package("deconica", quiet = TRUE),"/data-vignettes/mix1_ica.RData"))
```

Subsequently, we compute correlation between components and the original cell profiles.

```{r}
basis.list <- make_list(mix1$basis_matrix)

mix1_corr.basis <-
  correlate_metagenes (mix1_ica$S, mix1_ica$names, 
                       metagenes = basis.list,
                       orient.long = TRUE,
                       orient.max = FALSE)

```

```{r}
mix1_corr.basis_p <- radar_plot_corr(mix1_corr.basis, size.el.txt = 10, point.size = 2)
```
```{r}
mix1_corr.basis_p$p
```

We automatically assign a component to a cell type.
```{r}
mix1.assign <- assign_metagenes(mix1_corr.basis, exclude_name = NULL)
```
We use top 10 genes as signatures.
```{r}
mix1_ica.10 <-
  generate_markers(mix1_ica, 10,  sel.comp = as.character(mix1.assign[, 2]))
```
It is possible to visualize *the basis matrix* as a heatmap 
```{r}
mix1_ica.10.basis <-
  generate_basis(
    df = mix1_ica,
    sel.comp = as.character(mix1.assign[,2]),
    markers = mix1_ica.10
  )
```
```{r}
pheatmap::pheatmap(mix1_ica.10.basis, fontsize_row = 8)
```

We compute scores which are by default simple mean value of expression of the top genes in the original gene matrix.
```{r}
scores.mix1.ica <- 
  get_scores(mix1_ica$log.counts, mix1_ica.10)
```
We can see on the correlation plot almost prefect correspondence with the original proportions (`mix1$prop`)

```{r eval = FALSE}
scores_corr_plot(scores.mix1.ica,t(mix1$prop),  method="number",  tl.col = "black")
```
```{r echo = FALSE}
corr.plot <- scores_corr_plot(scores.mix1.ica,t(mix1$prop),  method="number",  tl.col = "black")
```


```{r, echo =FALSE}
corrplot::corrplot(
    corr.plot$corr.filtered,
    is.corr = FALSE,
    tl.cex = 1 / graphics::par("cex"),
    cl.cex = 1 / graphics::par("cex"),
    method= "number",
    tl.col = "black"
    )
```

One important feature of the ICA is that even if we do not know the exact number of components, when we overestimate the number of components, the signals are not altered. Our team proposed a method called MSTD that was developed for cancer transcriptomes to estimate optimal number of components [@Kairov2017]. Here we just slightly overestimate the number of components to 20.

Again you can load the decompostion `mix1_ica.20.RData` from `data-vignettes`repository of the package 
```{r, eval=FALSE}
mix1_ica.20 <- run_fastica (
  mix1$expression,
  overdecompose = FALSE,
  with.names = FALSE,
  gene.names = row.names(mix1$expression),
  n.comp = 20,
  R = FALSE
  )
```

```{r echo=FALSE}
load(paste0(path.package("deconica", quiet = TRUE),"/data-vignettes/mix1_ica.20.RData"))
```
Then we repeat the pipeline...
```{r}
mix1_corr.basis.20 <-
  correlate_metagenes (
    mix1_ica.20$S,
    mix1_ica.20$names,
    metagenes = basis.list,
    orient.long = FALSE,
    orient.max = TRUE
  )
mix1_corr.basis.20_p <- radar_plot_corr(mix1_corr.basis.20,
                size.el.txt = 10,
                point.size = 2)

```
```{r}
mix1_corr.basis.20_p$p
```
```{r}
mix1.assign.20 <- assign_metagenes(mix1_corr.basis.20, exclude_name = NULL)
mix1_ica.20$S <- mix1_corr.basis.20$S.max
mix1_ica.20_10 <- generate_markers(mix1_ica.20, 10,  sel.comp = as.character(mix1.assign.20[,2]),  orient.long = FALSE)
scores.mix1.ica.20 <- 
  get_scores(mix1_ica.20$log.counts, mix1_ica.20_10)
```
```{r eval= FALSE}
scores_corr_plot(scores.mix1.ica.20,t(mix1$prop),  method="number",  tl.col = "black")
```

```{r echo = FALSE}
corr.plot.20 <- scores_corr_plot(scores.mix1.ica.20,t(mix1$prop),  method="number",  tl.col = "black")


```

```{r echo = FALSE}
corrplot::corrplot(
    corr.plot.20$corr.filtered,
    is.corr = FALSE,
    tl.cex = 1 / graphics::par("cex"),
    cl.cex = 1 / graphics::par("cex"),
    method= "number",
    tl.col = "black"
    )
```

And we observe that the estimated proportions are correctly estimated.


### *in vitro* mixtures of immune cells

In this demo we use data published in ([@Becht2016])[https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE64385] `GSE54385`that you can download yourself directly from GEO using following lines of code.
```{r, eval = FALSE}
library(Biobase)
library(GEOquery)
library(limma)

# load series and platform data from GEO
GSE64385 <- getGEO("GSE64385", GSEMatrix =TRUE, AnnotGPL=TRUE)[[1]]
```
Or you can load `GSE64385.RData` from `data-vignettes` repo.

This dataset contains 5 immune cell types sorted from 3 healthy donors' peripheral bloods and mixed at different proportions. 
```{r echo=FALSE}
load(paste0(path.package("deconica", quiet = TRUE),"/data-vignettes/GSE64385.RData"))
```
```{r}
head(exprs(GSE64385))
```
Here is the raw matrix of proportions.
```{r}
cell_prop <- pData(GSE64385)[ , c(1,2, 10, 11,12, 13, 14, 15, 16, 17)]
kable(cell_prop, "html", row.names = TRUE) %>%
  kable_styling(font_size = 8)
```
We manually extracted the immune cell proportions from the matrix
```{r}
cell_prop.clean <- cell_prop [,1:2]
cell_prop.clean$NK <- c(0,0, 10, 5, 2.5, 1.3, 0.6, 10, 0.6, 1.3, 2.5, 5)
cell_prop.clean$Bcell <-  c(0,0, 0.6, 10, 5, 2.5, 1.3, 5, 10, 0.6, 1.3,2.5)
cell_prop.clean$Neutrophils <- c(0,0, 0.3, 0.2, 2.5, 1.3, 0.6, 0.6, 1.3, 2.5, 0.2, 0.3 )
cell_prop.clean$Tcell <- c(0,0,2.5, 1.3, 0.6, 10, 5, 1.3, 2.5, 5,10,0.6 )
cell_prop.clean$Monocytes <- c(0,0, 5, 2.5, 1.3, 0.6, 10, 0.6, 1.3, 2.5, 5, 10)
kable(cell_prop.clean, "html", row.names = TRUE) %>%
  kable_styling(font_size = 8)

```
Then we scaled them to relative proportions.
```{r}
rowSums(cell_prop.clean[, 3:7])
cell_prop.clean.scaled <-
  cell_prop.clean[, 3:7] / rowSums(cell_prop.clean[, 3:7])
kable(cell_prop.clean.scaled, "html", row.names = TRUE) %>%
  kable_styling(font_size = 8)
```
Then we performed ICA decomposing into 6 components as we account for the possible junk component.
```{r eval = FALSE}
GSE64385_ica <- run_fastica (
  exprs(GSE64385),
  overdecompose = FALSE,
  with.names = FALSE,
  gene.names = row.names(exprs(GSE64385)),
  samples = colnames(exprs(GSE64385)),
  n.comp = 6,
  R = FALSE
)
```
```{r echo=FALSE}
load(paste0(path.package("deconica", quiet = TRUE),"/data-vignettes/GSE64385_ica.RData"))
```
The components are oriented according to *long tail* even without verifying correlations, once we select the top 10 genes of each component we can generate scores.
```{r}
GSE64385_ica_markers_10 <- generate_markers(GSE64385_ica, 10)
GSE64385.ica_scores <- 
  get_scores(GSE64385_ica$log.counts, GSE64385_ica_markers_10)
```
```{r eval= FALSE}
scores_corr_plot(GSE64385.ica_scores[,1:5],cell_prop.clean.scaled,  method="number",  tl.col = "black")
```
```{r echo = FALSE}
log.corr <- scores_corr_plot(GSE64385.ica_scores[,1:5],cell_prop.clean.scaled,  method="number",  tl.col = "black")
```
```{r echo = FALSE}
corrplot::corrplot(
    log.corr$corr.filtered,
    is.corr = FALSE,
    tl.cex = 1 / graphics::par("cex"),
    cl.cex = 1 / graphics::par("cex"),
    method= "number",
    tl.col = "black"
    )
```

And even better accuracy can be obtained if we compute scores on *un-logged* data, actual counts. 

```{r}
GSE64385.ica_scores_unlog <- 
  get_scores((2^GSE64385_ica$log.counts)-1, GSE64385_ica_markers_10)
```
```{r eval = FALSE}
scores_corr_plot(GSE64385.ica_scores_unlog[,1:5],cell_prop.clean.scaled[1:5],  method="number",  tl.col = "black")
```
```{r echo = FALSE}
unglog.corr <- scores_corr_plot(GSE64385.ica_scores_unlog[,1:5],cell_prop.clean.scaled[1:5],  method="number",  tl.col = "black")
```
```{r echo = FALSE}
corrplot::corrplot(
    unglog.corr$corr.filtered,
    is.corr = FALSE,
    tl.cex = 1 / graphics::par("cex"),
    cl.cex = 1 / graphics::par("cex"),
    method= "number",
    tl.col = "black"
    )
```

### Blood data paired with FACS estimated proportions

Here we will use `SDY420` dataset from [Immport](http://www.immport.org/immport-open/public/home/home) database that @Aran2017 kindly shared with us. They contain PBMC  expression data of 104 healthy patients and paired FACS proportion estimation.
```{r echo = FALSE}
#import expression data
GE_SDY420 <- read.delim(paste0(path.package("deconica", quiet = TRUE),"/data-raw/xCell_ImmPort/GE_SDY420.txt"), row.names=1, stringsAsFactors=FALSE)
```
```{r eval = FALSE}
#import expression data
GE_SDY420 <- read.delim("./data-raw/xCell_ImmPort/GE_SDY420.txt", row.names=1, stringsAsFactors=FALSE)
```
```{r}
dim(GE_SDY420)
```
```{r}
summary(GE_SDY420)[,1:3]
```

In order to define number of components we can perform `doICABatch` that scans an array of decompositions looking for the most reproducible ones. The methodology was published in [@Kairov2017].

The function saves on the disk all the studied decomposition and generates plots of stability.

```{r eval =  FALSE,  out.height='269.5px', out.width='325px'}
GE_SDY420_batch.res<-doICABatch(GE_SDY420,
           seq(2,100,1),
           names = row.names(GE_SDY420),
           samples = colnames(GE_SDY420))
```

```{r echo = FALSE, fig.align='center', out.height='359.3px', out.width='433.3px'}
knitr::include_graphics("./figures-ext/MSTD_estimate.png")
```
```{r echo = FALSE, fig.align='center', out.height='359.3px', out.width='433.3px'}
knitr::include_graphics("./figures-ext/AverageStability.png")
```

The MST = 39 indicates most reproducible number of components. Let's verify if among 39 components we find components associated with the immune cells.

You can load the file with decomposition from `data-vignettes`repo.

```{r eval = FALSE}
GE_SDY420_ica_39 <- run_fastica (
  GE_SDY420,
  gene.names = row.names(GE_SDY420),
  samples = colnames(GE_SDY420),
  overdecompose = FALSE,
  with.names = FALSE,
  n.comp = 39,
  R = FALSE
)

```


```{r echo = FALSE}
load(paste0(path.package("deconica", quiet = TRUE),"/data-vignettes/GE_SDY420_ica_39.RData"))
```

Then we perform our pipeline comparing to pure immune profiles from [@Newman2015].
```{r}
GE_SDY420_ica_39.corr.LM22 <-
  correlate_metagenes (GE_SDY420_ica_39$S, GE_SDY420_ica_39$names, metagenes = LM22.list)

GE_SDY420_ica_39.corr.LM22.p <-
  radar_plot_corr(GE_SDY420_ica_39.corr.LM22,
                  ax.size = 12,
                  size.el.txt = 22,
                  point.size = 7)
```

```{r, fig.dim = c(25, 35)}
GE_SDY420_ica_39.corr.LM22.p$p
```
```{r}
GE_SDY420_ica_39.LM22.reciprocal.corr <-
  assign_metagenes(GE_SDY420_ica_39.corr.LM22, exclude_name = NULL)
```
```{r}
kable(GE_SDY420_ica_39.LM22.reciprocal.corr, "html", row.names = FALSE)
```


```{r}
GE_SDY420_ica_39.LM22.max.corr <- get_max_correlations(GE_SDY420_ica_39.corr.LM22)
```
```{r}
kable(GE_SDY420_ica_39.LM22.max.corr, "html", row.names = FALSE)
```

Both reciprocal and maximal correlations indicate a set of components that can be labelled as immune cells. Let's verify if we can use them to estimate proportions of those cell types.

```{r}
SDY420_markers_10 <-
  generate_markers(
    df = GE_SDY420_ica_39,
    n = 10,
    sel.comp = as.character(unique(GE_SDY420_ica_39.LM22.max.corr$IC)),
    return = "gene.list"
  )
```
```{r}
GE_SDY420_ica_39_scores <-
  get_scores ((2 ^ GE_SDY420_ica_39$log.counts) - 1,
  SDY420_markers_10,
  summary = "mean",
  na.rm = TRUE
  )
```

```{r}
head(GE_SDY420_ica_39_scores)
```

As we have FACS measured proportions we can import them.

```{r echo = FALSE}
#import facs estimated proportions
FACS_SDY420 <-  read.delim(paste0(path.package("deconica", quiet = TRUE),"/data-raw/xCell_ImmPort/FCS_SDY420.txt"), row.names=1, stringsAsFactors=FALSE)
```
```{r eval = FALSE}
#import facs estimated proportions
FACS_SDY420 <- read.delim("../data-raw/xCell_ImmPort/FCS_SDY420.txt", row.names=1, stringsAsFactors=FALSE)
```
```{r}
dim(FACS_SDY420)
```
```{r}
common.samples <- intersect(colnames(FACS_SDY420), colnames(GE_SDY420))
length(common.samples)
```
```{r}
FACS_SDY420.fil <- data.frame(t(FACS_SDY420[, common.samples]))
```
```{r}
head(FACS_SDY420.fil)[,1:4]
```

And we can confront the abundance scores.

```{r eval = FALSE}
scores_corr_plot(GE_SDY420_ica_39_scores,FACS_SDY420.fil,  tl.col = "black")
```
```{r}
knitr::include_graphics("./figures-ext/CorrBlood.jpeg")
```

The estimation of abundance gives remarkable accuracy (Pearson correlation coefficient):

* IC38: B-cells: 0.76 
* CD4 T-cells: 0.629
* CD8 T-cells: 0.63
* NK: 0.43
* Monocytes: 0.42

This results is also highly comparable with (or even better than) results obtained with other tools, published in *xCell* publication by @Aran2017.

```{r, fig.align='center', echo= FALSE}
knitr::include_graphics("figures-ext/SDY420.comp.png")
```

## Overview of functions

In this section we will discuss main functions of the package and their different options on computationally light toy examples. If you are interested in demo application of `deconICA` with biological arguments go straight to section [Demonstration of DeconICA package](#cases)


You can use `run_fastica()` function to decompose a matrix into independent components
```{r}
S <- matrix(runif(10000), 5000, 2)
A <- matrix(c(1, 1, -1, 3), 2, 2, byrow = TRUE)
X <- data.frame(S %*% A)
res <- run_fastica(X = X, row.center = TRUE, n.comp = 2, overdecompose = FALSE)
str(res)
```
`run_fastica` runs `fastica` from `fastica` package. In this trivial example we create sources matrix `S`and mixing matrix `A` that we multiply to obtain `X`. Then we decompose `X` into `n.comp = 2`, with row centering (substracting mean from each row) `row.center = TRUE`. We also checked `overdecompose = FALSE`, `overdecompose = TRUE` would ignore number of components we defined with `ncomp`. It finds its use for more advanced analysis applied to transcriptome [see section [Demonstration of DeconICA package](#cases)]. Other parameters were selected as default. 

Full description of the `run_fastica` parameters can be found in help `?run_fastica`. 

The main differences between `run_fastica` and `fastica` are:

* if column names are provided with the matrix, duplicated names are removed and entries with higher variance are kept
* it transforms data into log2 if data are in row counts 
* it runs a PCA before ICA to denoise the matrix
* it allows running matlab version fastica with *icasso* stabilisation if matlab software is installed on your machine (more about this point in [vignette: Running fastICA with icasso stabilisation](Icasso.html))
* it returns in a `list` 
  + input `data.frame` without duplicated entries and before log transformation: `log.counts`
  + `names` row names vector
  + `samples`sample names vector
  + `A`, `S`, `K` and `W` matrices (see `?run_fastica` for details) if run in `R=TRUE`
  + `A` and `S` matrices if run with `R=FALSE`
* `overdecompose` parameter that selects number of composed needed to perform overdecomposition of the input matrix
  
Therefore, `run_fastica()` performs ICA decomposition of the matrix and provides additional features useful for the downstream analysis. The use of more advanced options will be demonstrated later on in this tutorial. It generates *components* or *sources* to which we will refer later on in the tutorial.

The step naturally following `run_fastica()` is `correlate_metagenes()`.

It is common that after an unsupervised decomposition, components should have attributed an interpretation or a  meaning or a label. We call this process *interpret* components or *identify* sources. A domain knowledge is necessary to interpret components. In the case of transcriptomic data, components can be seen as weighted gene list. 

An efficient way to in interpret a component is to use correlation with some known profile or as we call it a *metagene*. If we dispose of a known weighted list of genes that characterize a biological phenomena or a cell type (a metagene), we can then correlate them with obtained components and verify is some of decomposed sources are close to the known cells/functions. 

In our `deconICA` pipeline `correlate_metagenes()` can be used in order to correlate metagenes (gene lists: knowlege-based or data-driven) with the data-driven components.

```{r}
library(deconica)
data(Example_ds)
#decompose the matrix
set.seed(123)
res_run_ica <- run_fastica (
  Example_ds,
  overdecompose = FALSE,
  with.names = TRUE,
  n.comp = 10
)

#correlate with Biton et al. metagenes
corr <- correlate_metagenes(S = res_run_ica$S,
                          gene.names = res_run_ica$names, 
                          metagenes = Biton.list
                           )
```
In this case we use an example of an extract, of 60 samples, from transcriptomic data from breast cancer [@Bekhouche2011]. At first, dataset is decomposed into an arbitrary number of 10 components. The `correlate_metagenes()` correlates the obtained `S` matrix with *Biton et al.* metagenes [@Biton2014]. This set of 11 metagenes is data-driven and was derived in the article @Biton2014 from pan-cancer transcriptome as the reproducible signals, common to many cancer types. Some of the signatures as BCLAPATHWAYS or UROTHELIALDIF (urothelialdifferentiation) were shared among many datasets, but within bladder cancer. They can bees as a sort of negative control.

However, any set weighted signatures with reasonable size can be used as metagenes. Later in this tutorial we use immune cell profiles optimized for cell type deconvolution from [@Newman2015] `LM22.list`.

One can controle a minimal number of genes used for correlation with `n.genes.intersect` it is set to a magic number in statistics (30) by default.

```{r}
names(corr)
```
The `correlate_metagenes()` returns `n`, `r`, `P` matrices which correspond to [`Hmisc::rcorr`](https://www.rdocumentation.org/packages/Hmisc/versions/4.1-0/topics/rcorr) function output, number of genes on which correlation is based, correlation coefficient and p-value.

The `S.or` stands for `S` matrix that is *oriented*. Why the matrix should be oriented? If you used ICA to decompose gene expression then the positive and negative projections do not have meaning by itself. We developed a methodology orienting data in the direction of the *long tail* of the distribution. Which means highest absolute values of a component weight should be positive. 
However, if the distribution doesn't have *tails* and is close to gaussian, an alternative, can be orienting components through maximal correlation. If we are confident with our metagenes, we can orient the `S` matrix so that the maximal correlation is always positive. We demonstrate the use of `orient.max = FALSE` and obtained `S.or`in [use cases section](#cases). One can decide not to orient the `S` matrix through selecting `orient.long = FALSE` and `orient.max = FALSE`.

If we have a look at `r` matrix, we see, indeed, it contains correlations between components and provided metagenes.
```{r}
head(corr$r)
```

These correlation matrix can be visualized in many ways. We propose a *radar plot* to evaluate quickly if there is a good match between metagenes and components. In this representation, we focus on positive correlations that have red color and placed away from the center of the circle (radar). We can see which component is attributed to the highest correlation for each metagene. 

```{r}
p <- radar_plot_corr(corr, size.el.txt = 10, point.size = 2)
```
```{r,echo =FALSE}
p$p
```

The function  `radar_plot_corr()` returns as well the matrix in long format suitable for `ggplot2` plots in case you want to use a different type of visualization.

```{r}
head(p$df)
```
In order to *zoom in* into a correlation with a specific metagene, one can use a function `lolypop_plot_corr()`

Here we can visualize for example SMOOTH_MUSCLE metagene that seems a bit ambigous.
```{r eval = FALSE}
lolypop_plot_corr(corr$r,"M3_SMOOTH_MUSCLE")
```
```{r echo = FALSE}
p1 <-lolypop_plot_corr(corr$r,"M3_SMOOTH_MUSCLE")
```
```{r, echo = FALSE}
p1
```

We can observe that the highest correlation is with IC4 0.44 but  IC1: 0.28 is close as well. These two components may represent two different types of functions or cells reated to muscles.

However, if we look at IMMUNE metagene, we can select one strongly correlated component: IC9, Pearson correlation coefficient equal to 0.73.
```{r eval = FALSE}
lolypop_plot_corr(corr$r,"M8_IMMUNE")
```
```{r echo=FALSE}
p2 <- lolypop_plot_corr(corr$r,"M8_IMMUNE")

```
```{r echo=FALSE}
p2
```

In case we have many profiles, an automatic extraction of corresponding pairs metagene - component is handy. 

A natural way is consider the component which Pearson correlation coefficient is highest as the one *corresponding* to a metagene.
You can use `get_max_correlations()` to retrieve this information from the correlation matrix. 
```{r}
# retrieve max correlations
 max.corr <- get_max_correlations(corr)
# order
 max.corr.ordered <- max.corr[order(-max.corr$r),]
# show table
kable(max.corr.ordered,"html", row.names = FALSE)
```

`get_max_correlations()` provides pearson correlation `r` column and the p-value `p.val` to help decide if the maximal correlation can be used as labelling. One can decide on minimal threshold, or p-value to take a decision.

Another way to assign metagene to a component can be trough reciprocal correlations. This method was used in our research articles [@Becht2016; RBH_paper]. In brief, given correlations between the set of metagenes $M=\{M_{1},...,M_{m}\}$ and $S$ matrix $S=\{IC_{1},...,IC_{N}\}$, if $S_{i}=argmax_k(corr(M_{j},S_{k}))$ and $A_{j}=argmax_k (corr(S_{i},M_{k}))$, then $S_{i}$ and $M_{j}$ are reciprocal. This approach is useful with assumption that there should be one component corresponding to a metagene and one metagene corresponding to a component.

```{r}
# retrieve reciprocal correlations
reciprocal.corr <- assign_metagenes(corr, exclude_name = NULL)
# show table
kable(reciprocal.corr, "html",row.names = FALSE)
```

Here the corresponding pairs do not follow a specific order. The six of components find a reciprocal match.

### Full pipeline example

#### Decompose data

In order to fully explore a dataset with identification and quantification of immune-related signals we suggest to over decompose the data matrix. 

We recommend to use for this purpose MATLAB implementation of the algorithm (see [vignette: Running fastICA with icasso stabilisation](Icasso.html))

One can run it like this on BEK complete data (WARNING: requires MATLAB and take a few minutes): 
```{r eval = FALSE}
library(deconica)
BEK_ica_overdecompose <- run_fastica (
  BEK,
  isLog = FALSE,
  overdecompose = TRUE,
  with.names = FALSE,
  gene.names = row.names(BEK),
  R = FALSE
)
```

In order to follow this tutorial simply load precomputed ICA decomposition.

```{r}
data(BEK_ica_overdecompose)
```

#### Interpret components

```{r}
# correlate with Biton et al. metagenes
corr_Biton <- correlate_metagenes(S = BEK_ica_overdecompose$S,
                            gene.names = BEK_ica_overdecompose$names,
                            metagenes = Biton.list
                            )
# correlate with LM22 cell profiles
corr_LM22 <- correlate_metagenes(S = BEK_ica_overdecompose$S,
                            gene.names = BEK_ica_overdecompose$names,
                            metagenes = LM22.list
                            )

```

We can illustrate correlations with reference metagenes...
```{r, eval = FALSE}
radar_plot_corr(corr_Biton, size.el.txt = 10, point.size = 2)
```
```{r include=FALSE}
p3 <- radar_plot_corr(corr_Biton, size.el.txt = 10, point.size = 2)
```
```{r, echo = FALSE}
p3$p
```

... and with LM22 immune cell type profiles.
```{r, eval = FALSE}
radar_plot_corr(corr_LM22, size.el.txt = 10, point.size = 2)
```
```{r include=FALSE}
p4 <- radar_plot_corr(corr_LM22, size.el.txt = 30, point.size = 8, ax.size = 10)
```
```{r echo=FALSE, fig.dim = c(25, 35), out.width='130%'}
p4$p
```

And zoom in the *M8_IMMUNE* metagene.
```{r, eval = FALSE}
lolypop_plot_corr(corr_Biton$r,"M8_IMMUNE")
```
```{r echo=FALSE}
p5 <- lolypop_plot_corr(corr_Biton$r,"M8_IMMUNE")
```
```{r, echo=FALSE}
p5
```

As we have several components corresponding to the M8_IMMUNE and a few matches with immune profiles. One can use enrichment test to confirm this results.

First we run reciprocal assignment so that we can exclude from immune signals the confounders as, for example, cell cycle.

```{r}
 reciprocal.corr.Biton <-
   assign_metagenes(corr_Biton, exclude_name = c("M8_IMMUNE", "M2_GC_CONTENT"))
 immune.components <-
   identify_immune_comp(corr_Biton$r[, "M8_IMMUNE"], reciprocal.corr.Biton$component)
```
Five components pass the threshold of >0.1 pearson correlation with the IMMUNE component.
```{r}
 immune.components
```

We can verify to which cells they may correspond through enrichment test (based on Fisher exact test).
```{r}
enrichment.immune <- gene_enrichment_test(
  BEK_ica_overdecompose$S,
  BEK_ica_overdecompose$names,
  immune.ics = names(immune.components),
  alternative = "greater",
  p.adjust.method = "BH",
  p.value.threshold = 0.05
  )
```

```{r}
names(enrichment.immune)
```
Output of `gene_enrichment_test()` is a list of three different elements for each component passing fixed p-value threshold. `metagenes` is a weighted list of top *n* genes, `genes.list` is character vector of top *n* genes and `enrichment`is enrichment test result
```{r}
enrichment.immune$metagenes$IC4[1:10,]
```
```{r}
enrichment.immune$genes.list$IC4[1:10]
```

```{r}
kable(enrichment.immune$enrichment$IC4[1:3,], "html", row.names = FALSE) %>%
  kable_styling(font_size = 8)
```
If you use ImmgenHUGO list of signatures, `cell_voting_immgen()`can be used to summarize results.
```{r}
kable(
  cell_voting_immgen(enrichment.immune$enrichment)$IC4,
  "html",
  row.names = FALSE,
  caption = "IC4"
  )
```
```{r}
kable(
  cell_voting_immgen(enrichment.immune$enrichment)$IC28,
  "html",
  row.names = FALSE,
  caption = paste("IC28")
  )
```
```{r}
kable(
  cell_voting_immgen(enrichment.immune$enrichment)$IC39,
  "html",
  row.names = FALSE,
  caption = paste("IC39")
  )
```
```{r}
kable(
  cell_voting_immgen(enrichment.immune$enrichment)$IC49,
  "html",
  row.names = FALSE,
  caption = paste("IC49")
  )
```
```{r}
kable(
  cell_voting_immgen(enrichment.immune$enrichment)$IC68,
  "html",
  row.names = FALSE,
  caption = paste("IC68")
  )
```

This result is not trivial to interpret.

A different .gmt can be also used to perform the enrichment analysis.

We can import for example signature genes from *TIMER* [@Li2016]

```{r eval = FALSE}
setwd(path.package("deconica", quiet = TRUE))
TIMER <-
  ACSNMineR::format_from_gmt("./data-raw/TIMER_cellTypes.gmt")
```

```{r, echo = FALSE}
TIMER <- deconica::TIMER_cellTypes
```

```{r}
enrichment.immune <- gene_enrichment_test(
  BEK_ica_overdecompose$S,
  BEK_ica_overdecompose$names,
  immune.ics = names(immune.components),
  gmt = TIMER,
  alternative = "greater",
  p.adjust.method = "BH",
  p.value.threshold = 0.05
  )
```
```{r}
kable(enrichment.immune$enrichment$IC4, "html", row.names = FALSE, caption = "IC4") %>%
  kable_styling(font_size = 8)
kable(enrichment.immune$enrichment$IC28, "html", row.names = FALSE, caption = "IC28") %>%
  kable_styling(font_size = 8)
kable(enrichment.immune$enrichment$IC39, "html", row.names = FALSE, caption = "IC39") %>%
  kable_styling(font_size = 8)
kable(enrichment.immune$enrichment$IC49, "html", row.names = FALSE, caption = "IC49") %>%
  kable_styling(font_size = 8)
kable(enrichment.immune$enrichment$IC68, "html", row.names = FALSE, caption = "IC68") %>%
  kable_styling(font_size = 8)
```

Based on this results we can suppose that:

* IC4 is related to B-cells
* IC49 and IC68 are related to Myeloid cells 
* IC39 is related to T-cells
* IC28 can be characterized as stroma

This result can be cross-verified with radar plot and maximal correlations with LM22

```{r}
#retrieve max correlations
 max.corr <- get_max_correlations(corr_LM22)
# order
 max.corr.ordered <- max.corr[order(-max.corr$r),]
# show table
kable(max.corr.ordered, "html",row.names = FALSE)

```

Indeed from correlation with LM22 we learn that

* IC4 is related to B-cells
* IC39 to T-cells
* IC49 to M2 Macrophages, Monocytes and Neutrophils
* IC68 to NK

#### Estimate abundance 

Once we identified components related to the immune cells, we can use them to extract  *specific* markers, best number of markers usually vary from 10 to 30

```{r}
markers_10 <-
  generate_markers(
    df = BEK_ica_overdecompose,
    n = 10,
    sel.comp = names(immune.components),
    return = "gene.list"
  )
```
```{r}
markers_10
```

```{r fig.height=8}
basis_10 <-
  generate_basis(
    df = BEK_ica_overdecompose,
    sel.comp = names(immune.components),
    markers = markers_10
  )
pheatmap::pheatmap(basis_10, fontsize_row = 8)
```
We can observe how it looks for `n = 30`(30 specific genes for each cell type - component)

```{r}
markers_30 <-
  generate_markers(
    df = BEK_ica_overdecompose,
    n = 30,
    sel.comp = names(immune.components),
    return = "gene.list"
  )
```

```{r fig.height=13}
basis_30 <-
  generate_basis(
    df = BEK_ica_overdecompose,
    sel.comp = names(immune.components),
    markers = markers_30
  )
pheatmap::pheatmap(basis_30, fontsize_row = 6)
```

And finally to obtain the abundance scores, a function `get_scores()`can be used.

If you want to use only gene.list 
```{r}
head(get_scores (BEK_ica_overdecompose$counts, markers_10, summary = "mean", na.rm = TRUE))
```
If you want to compute weighted mean (using ICA weights), 
first transform to a list of  data frames:
```{r}
weighted.list_10 <- generate_markers(
  df = BEK_ica_overdecompose,
  n = 10,
  sel.comp = names(immune.components),
  return = "gene.ranked"
)
```

```{r}
weighted.scores <- get_scores (BEK_ica_overdecompose$counts, weighted.list_10, summary = "weighted.mean", na.rm = TRUE)
head(weighted.scores)
```
These scores are linearly correlated with proportion values of each cell type in samples.

We can rename columns according to labels we attributed

```{r}
colnames(weighted.scores) <- c("B-cell", "Stroma", "T-cells", "Myeloid", "NK")
```

Here we present results in a form of a stacked baplot where scores sum up to 1. The scores used for plot should not be used as final scores, there are just more convinient for visualization purposes.
```{r fig.width=10}
stacked_proportions_plot(t(weighted.scores))
```
```{r  fig.width=10}
library(MCPcounter)
MCP.counter.scores <- MCPcounter.estimate(BEK_ica_overdecompose$counts,featuresType="HUGO_symbols")
stacked_proportions_plot(MCP.counter.scores)
```
```{r  fig.width=10}
#without enothelial cells and Fibroblasts
stacked_proportions_plot(MCP.counter.scores[c(-10,-9),])
```

<!-- ## Simulate Data  -->

<!-- TBA -->
<!-- ## Vignette Info -->

<!-- Note the various macros within the `vignette` section of the metadata block above. These are required in order to instruct R how to build the vignette. Note that you should change the `title` field and the `\VignetteIndexEntry` to match the title of your vignette. -->

<!-- ## Styles -->

<!-- The `html_vignette` template includes a basic CSS theme. To override this theme you can specify your own CSS in the document metadata as follows: -->

<!--     output:  -->
<!--       rmarkdown::html_vignette: -->
<!--         css: mystyles.css -->

<!-- ## Figures -->

<!-- The figure sizes have been customised so that you can easily put two images side-by-side.  -->

<!-- ```{r, fig.show='hold'} -->
<!-- plot(1:10) -->
<!-- plot(10:1) -->
<!-- ``` -->

<!-- You can enable figure captions by `fig_caption: yes` in YAML: -->

<!--     output: -->
<!--       rmarkdown::html_vignette: -->
<!--         fig_caption: yes -->

<!-- Then you can use the chunk option `fig.cap = "Your figure caption."` in **knitr**. -->

<!-- ## More Examples -->

<!-- You can write math expressions, e.g. $Y = X\beta + \epsilon$, footnotes^[A footnote here.], and tables, e.g. using `knitr::kable()`. -->

<!-- ```{r, echo=FALSE, results='asis'} -->
<!-- knitr::kable(head(mtcars, 10)) -->
<!-- ``` -->

<!-- Also a quote using `>`: -->

<!-- > "He who gives up [code] safety for [code] speed deserves neither." -->
<!-- ([via](https://twitter.com/hadleywickham/status/504368538874703872)) -->
## References
