---
title: "SplicingFactory"
author: "Péter Szikora, Tamás Pór, Endre Sebestyén"
date: "`r Sys.Date()`"
abstract: >
  SplicingFactory is an R package for the analysis of alternative splicing
  isoform diversity, with expression levels generally coming from RNA-Seq data.
  Expression levels are generated by the usual tools, for example a STAR
  alignment followed by featureCounts, RSEM, salmon or kallisto. However, you
  need transcript isoform level data, where the expression estimates of the
  different isoforms for a gene are as specific, as possible. The package
  provides the ability to generate gene-level Shannon-entropy values and other
  diversity measures, such as the Gini- or Simpson-index. These measures can
  quantify transcript isoform diversity within samples or between conditions.
  Additionally, the package analyzes the isoform diversity data, looking for
  significant changes between conditions. A basic task during analysis is to use
  read counts or TPM to calculate a diversity measure and then perform a
  differential analysis to identify genes with a significant change between
  conditions. The vignette explains the use of the package with a typical
  workflow.
  
output:
  rmarkdown::html_document:
    highlight: pygments
    toc: true
    fig_width: 6
    fig_height: 4
vignette: >
  %\VignetteIndexEntry{SplicingFactory}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Standard workflow

SplicingFactory package version: `r packageVersion("SplicingFactory")`

## Input data structure

As input, the SplicingFactory package can use 5 different data types. Besides
matrices and data frames, you can use the output of tximport, DGELists and
SummarizedExperiments. The tximport R package is able to import files containing
transcript abundance estimates generated by tools such as Salmon, Kallisto or
RSEM. The tximport datasets are stored in a list by default, and they can serve
as input for SplicingFactory package.

DGELists are a list-based S4 class for storing read counts and associated
information, while SummarizedExperiment objects contain and describe
high-throughput expression level assays. SplicingFactory automatically uses the
necessary functions for all of these data types, on condition that they were not
previously modified (e.g. list elements renamed).

## Data arrangement

Besides the format, the input table (matrix, data.frame or the tabular
expression data extracted from other data structures) needs to be arranged
properly. Every row identifies a distinct transcript in your dataset, while
every column belongs to a distinct sample. The table can contain only numeric
values. Supplementing your table, the package will need a gene and a sample
vector, identifying the genes and the sample conditions in your analysis. The
gene vector assigns genes to every row in your data, that will be used to
aggregate transcript level expression information at the gene level. It is
important that the columns and the rows of the table (genes and samples) are in
the same order, as the gene and the sample vectors. In case of
`SummarizedExperiment`, the genes vector can be empty (but not exclusively)
as the package can automatically extract the necessary vector from the object.

## Example dataset

The package contains an example dataset called `tcga_brca_luma_dataset`. The
data was downloaded from The Cancer Genome Atlas on 12th of April, 2020. It
contains transcript level read counts for 300 pre-selected genes of 40 patients
with Luminal A type breast cancer (primary tumor and solid normal samples).
Transcript level expression was estimated with RSEM.

## Splicing diversity analysis

### Importing example data

```{r setup}
library("SplicingFactory")
library("SummarizedExperiment")

# Load dataset
data(tcga_brca_luma_dataset)

# Extract gene names
genes <- tcga_brca_luma_dataset[, 1]

# Extract read count data without gene names
readcounts <- tcga_brca_luma_dataset[, -1]

# Check read count dataset
dim(readcounts)

head(readcounts[, 1:5])
```

As a first step, before doing the diversity calculation, you might want to
filter out genes with a low expression.

### Transcript diversity calculation

We are going to use the `calculate_diversity` function to calculate two different
types of transcript diversity.There are several mandatory and optional pamaters
for the function.

* x - Input data, in various formats, discussed in more detail in the Input data
structure and Data arrangement section.
* genes - A vector with gene names used for aggregating the transcript level data.
* method - Method to use for splicing diversity calculation.
* norm - If set to `TRUE`, the entropy values are normalized to the number of
transcripts for each gene.
* tpm - In the case of a tximport list, you might want set the `tpm` argument.
As the default option is `FALSE`, the raw  read counts will be extracted from
your input data. Set it to `TRUE` if you want to use TPM values.
* assayno - An optional argument is `assayno`, which is useful if you are
planning to analyze a `SummarizedExperiment` input, containing multiple assays.
`assayno` is a numeric value, specifying the assay to be analyzed.
* verbose - Set it to `TRUE` if you want more detailed diagnostic messages.

To calculate Laplace entropy, where values are normalized between 0 and 1, use:

```{r laplace}
laplace_entropy <- calculate_diversity(readcounts, genes,  method = "laplace",
                                       norm = TRUE, verbose = TRUE)

head(assay(laplace_entropy)[, 1:5])
```

To calculate Gini index, you don't need to specify the `norm` argument, as the
Gini index is by definition ranges between 0 (complete equality) and 1 (complete
inequality).

```{r gini}
gini_index <- calculate_diversity(readcounts, genes, method = "gini",
                                  verbose = TRUE)

head(assay(gini_index)[, 1:5])
```

Both for the Laplace-entropy and Gini index calculation, the package returns a
`SummarizedExperiment` object, that you can investigate further with the
`assay` function.

The package automatically filters out genes with a single isoform, as splicing
diversity values can only be calculated for genes with at least 2 splicing
isoforms.

Some genes might show `NA` diversity values. This means that the expression was
zero for all transcript isoforms of the gene in these samples and the package
could not calculate any diversity value as there is no meaningful diversity for
genes which did not show any expression in your experiment. Lack of expression
might also be the result of technical issues.

To further analyze and visualize your data, you might do the following. To see
the distribution and density of the splicing diversity data, you can visualize
it by using ggplot2 from the tidyverse package collection.

```{r divplots}
library("tidyr")
library("ggplot2")

# Construct data.frame from SummarizedExperiment result
laplace_data <- cbind(assay(laplace_entropy),
                      Gene = rowData(laplace_entropy)$genes)

# Reshape data.frame
laplace_data <- pivot_longer(laplace_data, -Gene, names_to = "sample",
                             values_to = "entropy")

# Add sample type information
laplace_data$sample_type <- apply(laplace_data[, 2], 1,
                                  function(x) ifelse(grepl("_N", x),
                                                     "Normal", "Tumor"))

# Filter genes with NA entropy values
laplace_data <- drop_na(laplace_data)

# Update gene names and add diversity type column
laplace_data$Gene <- paste0(laplace_data$Gene, "_", laplace_data$sample_type)
laplace_data$diversity <-  "Normalized Laplace entropy"

# Construct data.frame from SummarizedExperiment result
gini_data <- cbind(assay(gini_index), Gene = rowData(gini_index)$genes)

# Reshape data.frame
gini_data <- pivot_longer(gini_data, -Gene, names_to = "sample",
                          values_to = "gini")

# Add sample type information
gini_data$sample_type <- apply(gini_data[, 2], 1,
                               function(x) ifelse(grepl("_N", x),
                                                  "Normal", "Tumor"))

# Filter genes with NA gini values
gini_data <- drop_na(gini_data)

# Update gene names and add diversity type column
gini_data$Gene <- paste0(gini_data$Gene, "_", gini_data$sample_type)
gini_data$diversity <-  "Gini index"

# Plot diversity data
ggplot() +
  geom_density(data = laplace_data, alpha = 0.3,
               aes(x = entropy, group = sample, color = diversity)) +
  geom_density(data = gini_data, alpha = 0.3,
               aes(x = gini, group = sample, color = diversity)) +
  facet_grid(. ~ sample_type) +
  scale_color_manual(values = c("black", "darkorchid4")) +
  guides(color = FALSE) +
  theme_minimal() +
  labs(x = "Diversity values",
       y = "Density")

# Mean entropy calculation across samples for each gene/sample type combination
laplace_entropy_mean <- aggregate(laplace_data$entropy,
                                  by = list(laplace_data$Gene), mean)
colnames(laplace_entropy_mean)[2] <- "mean_entropy"
laplace_entropy_mean <- as_tibble(laplace_entropy_mean)

# Add sample type information
laplace_entropy_mean$sample_type <- apply(laplace_entropy_mean[, 1], 1,
                                          function(x) ifelse(grepl("_Normal", x),
                                                             "Normal", "Tumor"))

# Add diversity type column
laplace_entropy_mean$diversity <-  "Normalized Laplace entropy"

# Mean gini calculation across samples for each gene/sample type combination
gini_mean <- aggregate(gini_data$gini, by = list(gini_data$Gene), mean)
colnames(gini_mean)[2] <- "mean_gini"
gini_mean <- as_tibble(gini_mean)

# Add sample type information
gini_mean$sample_type <- apply(gini_mean[, 1], 1,
                               function(x) ifelse(grepl("_Normal", x),
                                                  "Normal", "Tumor"))

# Add diversity type column
gini_mean$diversity <-  "Gini index"

ggplot() +
  geom_violin(data = laplace_entropy_mean, aes(x = sample_type, y = mean_entropy,
                                               fill = diversity),
              alpha = 0.6) +
  geom_violin(data = gini_mean, aes(x = sample_type, y = mean_gini,
                                    fill = diversity),
              alpha = 0.6) +
  scale_fill_viridis_d(name = "Diversity") +
  coord_flip() +
  theme_minimal() +
  labs(x = "Samples",
       y = "Diversity")
```

The two methods are calculating different results. Genes with a single dominant
isoform have a near-zero entropy, while they have a Gini index close to 1. The
overall distribution of the data is similar between the Normal and Tumor
conditions.

### Differential analysis

To further analyze the data, the steps of a differential diversity analysis are
implemented in the `calculate_difference` function, aiming to identify genes
with significant changes in splicing diversity. The result table will contain
the mean or median values of the diversity across sample categories, the
difference of these values, the log2 fold change of the two different
conditions, p-values and adjusted p-values for each genes.

There are several mandatory and optional arguments for this function:

* x - Output of the previous function, which is a data.frame with
gene names as the first column and splicing diversity values for
each sample in additional columns.
* samples - The previously defined vector with sample
conditions specifying the category of each sample in the correct order.
* control - Name of the control sample category, defined in the samples
vector, e.g. control = "Normal" or control = "WT".
* method - Method to use for calculating the average splicing diversity
value in a condition. Can be "mean" or "median".
* test - Method to use for p-value calculation. There are two statistical tests
implemented at the moment - Wilcoxon rank sum test and label shuffling test.
* randomizations - Optional parameter for the label shuffling test. You
can specify the number of random shuffles (default = 100).
* pcorr - P-value correction method to use for the test results. Benjamini &
Hochberg by default.
* assayno - An optional argument is `assayno`, which is useful if you are
planning to analyze a `SummarizedExperiment` input, containing multiple assays.
`assayno` is a numeric value, specifying the assay to be analyzed.
* verbose - Set it to `TRUE` if you want more detailed diagnostic messages.

To analyze the previously calculated normalized Laplace entropy values stored in
a `SummarizedExperiment` object with a Wilcoxon sum rank test, use the
`calculate_difference` function as follows:

```{r}
# Update the SummarizedExperiment object with a new sample metadata column for
# sample types, as the the object returned by calculate_diversity does not
# contain this information.
colData(laplace_entropy) <- cbind(colData(laplace_entropy),
                                      sample_type = ifelse(grepl("_N", laplace_entropy$samples),
                                                           "Normal", "Tumor"))

# Calculate significant entropy changes
entropy_significance <- calculate_difference(x = laplace_entropy, samples = "sample_type",
                                             control = "Normal",
                                             method = "mean", test = "wilcoxon",
                                             verbose = TRUE)

head(entropy_significance)
```

The package sends a note about 16 genes, with low sample size, excluded from the
analysis. as these genes had several NA diversity values, the result of 0
expression values. You need at least 3 samples in each sample category and a
total of 8 samples for a Wilcoxon test and at least 5 samples in each sample
category for the label shuffling.

Genes with a significant change in splicing diversity can be further analyzed
and visualized by using e.g. MA-plots, where the log2 fold changes are shown on
the y axis, and the mean diversity values on the x axis. Dots are colored red if
the adjusted p-value is smaller than 0.05. It is recommended to filter for an
absolute log2 fold change larger than 0.5 besides the adjusted p-value.

```{r}
entropy_significance$label <- apply(entropy_significance[, c(5, 7)], 1,
                                    function(x) ifelse(abs(x[1]) >= 0.5 & x[2] < 0.05,
                                                       "significant", "non-significant"))

entropy_significance$mean <- apply(entropy_significance[, c(2, 3)], 1,
                                   function(x) (x[1] + x[2]) / 2)


ggplot(entropy_significance, aes(x = mean, y = log2_fold_change)) +
  geom_point(color = "lightgrey", size = 1) +
  geom_point(data = entropy_significance[entropy_significance$label == "significant", ],
             color = "red", size = 1) +
  theme_minimal() +
  labs(title = "Normalized Laplace entropy",
       subtitle = "Wilcoxon signed rank test",
       x = "Mean entropy",
       y = "Log2(fold change)")
```

The analyzed genes also can be visualized on a Volcano-plot, which shows the
adjusted p-values of the genes on a logarithmic scale on the y axis and the log2
fold change values between the two conditions on the x axis. We used a cutoff of
0.5 and -0.5 for the log2 fold changes.

```{r}
ggplot(entropy_significance, aes(x = log2_fold_change,
                                 y = -log10(adjusted_p_values),
                                 color = label)) +
  geom_point() +
  scale_color_manual(values = c("grey", "red"), guide = "none") +
  geom_hline(yintercept = -log10(0.05), color = "red") +
  geom_vline(xintercept = c(0.5, -0.5)) +
  theme_minimal() +
  labs(title = "Normalized Laplace entropy",
       subtitle = "Wilcoxon signed rank test",
       x = "Log2 fold change of mean entropy values",
       y = "-Log10(adjusted p-value")
```

## Session info
```{r}

sessionInfo()
```