---
title: "VSClust with Bioconductor objects"
author: "Veit Schwämmle"
package: vsclust
abstract: >
  Here, we describe the workflow to run variance-sensitive clustering on data 
  stored in a SummarizedExperiment, QFeatures or MultiAssayExperiment object. 
  This vignette is distributed under a CC BY-SA license.
vignette: >
  %\VignetteIndexEntry{VSClust on Bioconductor object}
  %\VignetteEngine{knitr::rmarkdown}
  %%\VignetteKeywords{Genomics, Transcriptomics, Proteomics, Metabolomics, Clustering, Quantitative, Statistics }
  %\VignetteEncoding{UTF-8}
output:
  BiocStyle::html_document:
    toc_float: true
---

```{r style, echo = FALSE, results = 'asis'}
BiocStyle::markdown()
```

```{r env, message = FALSE, warning = FALSE, echo = FALSE}
library("vsclust")
library("MultiAssayExperiment")
```

# Introduction

For a more detailed explanation of the VSClust function and the workflow, please
take a look on the vignette for running the VSClust workflow. 

Here, we present an example script to integrate the clustering with data object 
from Bioconductor, such as `QFeatures`, `SummarizedExperiment` and
`MultiAssayExperiment`. 

# Installation and additional packages

Use the common Bioconductor commands for installation:

```{r eval=FALSE}
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("vsclust")
```
The full functionality can be obtained by additionally installing and loading the
packages `yaml`, `shiny`, `clusterProfiler`, and `matrixStats`.


# Initialization

Here, we define the different parameters for the data set `RNASeq2GeneNorm` from
the `miniACC` object. 

The number of replicates and experimental conditions will be retrieved 
automatically by specifying the metadata for the grouping. 
```{r}
#### Input parameters, only read when now parameter file was provided #####
## All principal parameters for running VSClust can be defined as in the 
## shiny app at computproteomics.bmb.sdu.dk/Apps/VSClust
# name of study
Experiment <- "miniACC" 
# Paired or unpaired statistical tests when carrying out LIMMA for 
# statistical testing
isPaired <- FALSE
# Number of threads to accelerate the calculation (use 1 in doubt)
cores <- 1 

# If 0 (default), then automatically estimate the cluster number for the 
# vsclust run from the Minimum Centroid Distance
PreSetNumClustVSClust <- 0 
# If 0 (default), then automatically estimate the cluster number for the 
# original fuzzy c-means from the Minimum Centroid Distance
PreSetNumClustStand <- 0 

# max. number of clusters when estimating the number of clusters. 
# Higher numbers can drastically extend the computation time.
maxClust <- 10 
```

# Statistics and data preprocessing

At first, we load will log-transform the original data and normalize it to the 
median. Statistical testing will be applied on the resulting object. 
After estimating the standard deviations, the matrix consists of the averaged 
quantitative feature values and a last column for the standard deviations of the 
features.

We will separate the samples according to their `OncoSign`.

```{r fig.width = 12 }
data(miniACC, package="MultiAssayExperiment")
# log-transformation and remove of -Inf values
logminiACC <- log2(assays(miniACC)$RNASeq2GeneNorm)
logminiACC[!is.finite(logminiACC)] <- NA
# normalize to median
logminiACC <- t(t(logminiACC) - apply(logminiACC, 2, median, na.rm=TRUE))

miniACC2 <- c(miniACC, log2rnaseq = logminiACC, mapFrom=1L)

boxplot(logminiACC)


#### running statistical analysis and estimation of individual variances
statOut <- PrepareSEForVSClust(miniACC2, "log2rnaseq", 
                               coldatname = "OncoSign", 
                               isPaired=isPaired, isStat=TRUE)

```

We can see that there is no good separation of cancer signatures on the PCA 
plot.

# Estimation of cluster number

There is no simple way to find the optimal number of clusters in a data set. For 
obtaining this number, we run the clustering for different cluster numbers and 
evaluate them via so-called validity indices, which provide information about 
suitable cluster numbers. VSClust uses mainly the "Maximum centroid distances" 
that denotes the shortest distance between any of the centroids. Alternatively, 
one can inspect the Xie Beni index.

The output of `estimClustNum` contains the suggestion for the number of 
clusters.

We further visualize the outcome.
```{r fig.width = 12}

#### Estimate number of clusters with maxClust as maximum number clusters to run 
#### the estimation with
ClustInd <- estimClustNum(statOut$dat, maxClust=maxClust, cores=cores)

#### Use estimate cluster number or use own
if (PreSetNumClustVSClust == 0)
  PreSetNumClustVSClust <- optimalClustNum(ClustInd)
if (PreSetNumClustStand == 0)
  PreSetNumClustStand <- optimalClustNum(ClustInd, method="FCM")
#### Visualize
  estimClust.plot(ClustInd)



```
Both validity indices agree with each other and suggest 7 as the most reasonable
estimate for the cluster number. However, we can also see that this decreases 
the number of clustered features quite drastically from over 150 to about 90.


# Run final clustering

Now we run the clustering again with the optimal parameters from the estimation.
One can 
take alternative numbers of clusters corresponding to large decays in the 
Minimum Centroid Distance or low values of the Xie Beni
index.

First, we carry out the variance-sensitive method

``` {r}


#### Run clustering (VSClust and standard fcm clustering
ClustOut <- runClustWrapper(statOut$dat, 
                            PreSetNumClustVSClust, NULL, 
                            VSClust=TRUE,
                            cores=cores)
Bestcl <- ClustOut$Bestcl
VSClust_cl <- Bestcl

```

We see how different groups of genes show distinctive pattern of their 
expression in different oncological signatures. 


```{r}
sessionInfo()
```