---
title: "Introduction to OutSeekR"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to OutSeekR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup}
library(OutSeekR);
```

# Introduction

The **OutSeekR** package implements a statistical approach to detecting outliers in RNA-seq and related data.

# Methodology

An analysis using **OutSeekR** should start with a matrix or data frame of normalized RNA-seq data, e.g., fragments per kilobase of transcript per million fragments mapped (FPKM), or similar data; in the case of RNA-seq data, we require normalized values rather than the counts themselves because the statistics we calculate assume continuous data. For simplicity, we'll refer to the rows of the input data set as 'transcripts' and the columns of the input data set as 'samples'. Transcript identifiers should be stored as the row names of the matrix or data frame, and sample identifiers should be stored as the column names.

The statistical approach of **OutSeekR** centers around the use of five statistics calculated for each transcript in the observed data:

1. the range of standard scores (z-scores) computed using the mean and standard deviation;
2. the range of z-scores computed using the median and median absolute deviation (MAD);
3. the range of z-scores computed using the 5%-trimmed mean and 5%-trimmed standard deviation;
4. the fraction of observations assigned to the smaller cluster based on K-means clustering with K = 2 clusters; and
5. the cosine similarity between the most extreme observed value and the largest quantile of a representative theoretical distribution (see [Simulating null data]).

Specifically, it uses the five statistics calculated on the observed data and compares the distributions of these statistics with counterparts calculated using simulated null data. Observed data yielding statistics more extreme than those of the null data suggest the presence of outliers.

Note: Input data with low variation (e.g., where 50% or more of the values are identical) may cause issues in calculating the MAD or performing K-means clustering. It is recommended to exclude such data from the analysis to ensure reliable results.
**OutSeekR** first identifies which transcripts/features contain outlier values.  For a given transcript, the null hypothesis is there are no outlier values (*i.e.* all values come from some unknown population distribution) while the alternative hypothesis is that the transcript contains 1 or more outlier values from another distribution.  A small p-value or FDR-value gives evidence for the presence of 1 or more outliers.

After determining which transcripts contain outliers, an iterative procedure is used to identify the exact number of outliers each transcript has.
## Overview

As an overview, the **OutSeekR** algorithm does the following:

- Simulates null data designed to share the characteristics of the observed data but without outlying values.
- Calculates the five statistics on each transcript in the observed data and each simulated transcript in the null data.
- For each transcript:
  - Combines the five statistics for the observed transcript with the five statistics from the null data.
  - Ranks each of the five vectors of statistics individually.
  - Computes the rank product for the observed transcript and each null transcript, defined as $$\sqrt[k]{\prod_{i=1}^k Rank_i},$$ where $k$ is the number of nonmissing ranks across the five statistics.
  - Computes a p-value as the proportion of rank products among the null transcripts greater than the rank product for the observed transcript. This p-value corresponds to the outlier test for the sample with the **most extreme value** in the observed transcript. If the p-value indicates this sample is an outlier, we exclude this sample, recalculate the outlier statistics in the remaining samples, combine these recalculated statistics with those computed from the null transcripts, recompute the rank product, and compute the new p-value. We continue in this fashion until we find a nonsignificant p-value (or, more precisely, a p-value greater than a user-specified threshold); each subsequent sample will be less outlying than the current one, so no more outliers are present, and we continue with the next transcript.
- Combines results across transcripts and performs p-value adjustment using the false discovery rate (FDR).

## Simulating null data

We hypothesize that the underlying data can be represented by one of four theoretical distributions: the normal, log-normal, exponential, or gamma distributions. For each transcript, we fit a generalized additive model for location, scale, and shape (GAMLSS) to the observed data based on each distribution using `gamlss()` from the [**gamlss**](https://www.gamlss.com/) package; the model with the smallest Bayesian Information Criterion corresponds to the optimal theoretical distribution for the transcript. To account for additional variability, we also calculate 'residuals' for each transcript as the differences between the empirical quantiles of the observed data and the quantiles of the optimal theoretical distribution. Using the same GAMLSS-based method, we identify the optimal distribution of these residuals from among the normal, log-normal, exponential, or gamma. (Note that the optimal distributions of a transcript and its residuals need not be the same.)

With the optimal distribution of each observed transcript and its residuals, we can create the simulated data. We iterate a user-specified number of times---1000 by default. At each step, we randomly select a transcript from the observed data and generate as many values as there are samples in the observed data according to its optimal distribution. We then add noise based on the optimal distribution of its residuals. The final null transcript is the sum of the simulated values and the simulated residuals. The final data set of simulated transcripts will contain as many rows as transcripts selected and the same number of columns as the observed data. While 1000 iterations are sufficient for exploratory analysis, we recommend using at least 10,000 iterations for publication-level results, with 100,000 or even one million iterations providing more robust estimates.

The null data do not contain outliers by construction but are still in some sense representative of the observed data. By calculating the five outlier statistics on these transcripts, we can get a sense of their distributions in data similar to that observed but without outliers; this logic forms the basis for the p-value of the outlier test.

## Parallelization in **OutSeekR**

The statistical approach to outlier detection implemented by **OutSeekR** is time-consuming. To reduce runtime, the code in **OutSeekR** has been written to allow parallelization. Parallelization is supported through the use of the package [**future.apply**](https://future.apply.futureverse.org/). **future.apply** is built on top of the future framework implemented by the [**future**](https://future.futureverse.org/) package, which provides a uniform way to parallelize code independent of operating system or computing environment. In particular, **OutSeekR** uses the `future_*apply()` functions defined in **future.apply** rather than their base R equivalents in order to take advantage of a parallelization strategy set with `future::plan()` by the user.

**Note:** Depending on the size of the input data set and the number of null transcripts requested, the objects created by **OutSeekR** may exceed a limit defined in **future** to prevent exporting very large objects, which will trigger an error with a message referring to `future.globals.maxSize`. This can be avoided by setting `future.globals.maxSize` to a large value, e.g., `options(future.globals.maxSize = 8000 * 1024^2);  # = 8 GB`. See [the documentation of `future.options`](https://search.r-project.org/CRAN/refmans/future/html/future.options.html) for further details.

## Output

The output of the main function in **OutSeekR**, `detect.outliers()`, is a named list with the following components:

- `p.values`: a matrix of unadjusted p-values for the outlier test run on each transcript in the observed data.
- `fdr`: a matrix of FDR-adjusted p-values for the outlier test run on each transcript in the observed data.
- `num.outliers`: a vector giving the number of outliers detected for each transcript based on the threshold.
- `outlier.test.results.list`: a list of length `max(num.outliers) + 1` containing entries `roundN`, where `N` is between one and `max(num.outliers) + 1`.  `roundN` is the data frame of results for the outlier test after excluding the $(N-1)$th outlier sample, with `round1` being for the original data set (i.e., before excluding any outlier samples).
- `distributions`: a numeric vector indicating the optimal distribution for each transcript. Possible values are 1 (normal), 2 (log-normal), 3 (exponential), and 4 (gamma).
- `initial.screen.method`: Specifies the statistical criterion for initial feature selection. Valid options are 'p-value' and 'FDR'.

# Example

**OutSeekR** includes a sample data set called `outliers` that we shall use as an example. It is a data frame consisting of 50 samples across 500 transcripts. Note that transcript identifiers and sample identifiers are stored as `rownames(outliers)` and `colnames(outliers)`, respectively.

```{r view-data}
str(outliers, list.len = 5);
outliers[1:6, 1:6];
outliers[495:500, 45:50];
```

An analysis in **OutSeekR** is run using the function `detect.outliers()`. With default settings, it can be run by simply passing the input data set via the `data` parameter. Other parameters include the following:

- `num.null`: the number of transcripts to generate when simulating from null distributions; default is 1000.
- `initial.screen.method`: the statistical criterion for initial gene selection.
- `p.value.threshold`: the p-value threshold for the outlier test; default is 0.05.
- `fdr.threshold`: the false discovery rate (FDR)-adjusted p-value threshold for determining the final count of outliers; default is 0.01.
- `kmeans.nstart`: the number of random starts when computing k-means fraction; default is 1. See `?stats::kmeans` for further details.

In our first example, we pass the `outliers` data frame and use the defaults for most other parameters. To keep runtime manageable, we limit the number of null transcripts generated to 1,000. Prior to running `detect.outliers()`, we set up parallel processing. See [Parallelization in **OutSeekR**] for further details.

```{r run-1}
# Set random seed for reproducibility.
set.seed(371892);

# Set up parallel processing.
future::plan(future::multisession);

outlier.test.run.1 <- detect.outliers(
    data = outliers,
    num.null = 1e3
    );

str(outlier.test.run.1, max.level = 2);

# Restore sequential processing.
future::plan(future::sequential);
```

Unadjusted p-values for the outlier test on each transcript are contained in the `p.values` entry of the object returned by `detect.outliers()`; FDR-adjusted p-values are contained in the `fdr` entry. 

```{r examine-p-values}
head(outlier.test.run.1$p.values);
head(outlier.test.run.1$fdr);
```

As explained in [Overview], an iterative procedure is used to determine the exact number of outliers per transcript.  The iterations stop after the p-value or FDR-value exceeds the specified threshold.

The entries `num.outliers` give transcript-specific counts of the number of outliers based on the threshold of specified `initial.screen.method`.

```{r transcript-level-outlier-counts}
head(outlier.test.run.1$num.outliers);
```

We can get a quick view of the distribution of the number of outliers across transcripts using `table()`:

```{r table}
table(outlier.test.run.1$num.outliers);
```

More granular results can be found in the `outlier.test.results.list` entry of the object returned by `detect.outliers()`. It is a named list of data frames, where each reflects the results of a different 'round' of outlier testing; by 'round', we mean the process of excluding the most outlying sample (for all rounds except the first), calculating the five outlier statistics, ranking them alongside the statistics of the null data, computing the rank product, and computing the unadjusted p-value. Each data frame in `outlier.test.results.list` contains the transcript and sample identifiers, the five statistics calculated on the observed transcript, the rank product, and the unadjusted and FDR-adjusted p-value for the outlier test.

```{r view-outlier-test-results-list}
str(outlier.test.run.1$outlier.test.results.list);
```

It may be of use to combine the results across all rounds of outlier testing into a single data frame. One way to do this (while retaining the round information) using `lapply()` and `rbind()` is shown below:

```{r collapse-rounds}
outlier.test.results.combined <- lapply(
    X = seq_along(outlier.test.run.1$outlier.test.results.list),
    FUN = function(i) {
        df <- outlier.test.run.1$outlier.test.results.list[[i]];
        df$round <- i;
        df <- df[, c(
            'round',
            colnames(outlier.test.run.1$outlier.test.results.list[[i]])
            )];
        }
    );
outlier.test.results.combined <- do.call(
    what = 'rbind',
    args = outlier.test.results.combined
    );
# Combining the data frames produces duplicates in the row names.  R
# will de-duplicate them, but as all the necessary information is
# included in the columns of the data frame (specifically, 'round' and
# 'transcript'), we'll simply discard the row names.
rownames(outlier.test.results.combined) <- NULL;
head(outlier.test.results.combined);
```

We can observe the optimal theoretical distribution for each transcript and the frequencies across all transcripts.

```{r distributions}
head(outlier.test.run.1$distributions);
table(outlier.test.run.1$distributions);
```

The results are numeric codes, with 1 corresponding to the normal distribution, 2 the log-normal distribution, 3 the exponential distribution, and 4 the gamma distribution.

Lastly, we show a second run of the algorithm where we adjust the p-value and FDR thresholds:

```{r run-2}
# Set up parallel processing.
future::plan(future::multisession);

outlier.test.run.2 <- detect.outliers(
    data = outliers,
    num.null = 1e3,
    initial.screen.method = 'fdr',
    p.value.threshold = 0.25,
    fdr.threshold = 0.05
    );

# Restore sequential processing.
future::plan(future::sequential);

str(outlier.test.run.2, max.level = 2);

# Examine p-value and FDR matrices.
head(outlier.test.run.2$p.values);
head(outlier.test.run.2$fdr);

# Check the distribution of number of outliers detected.
table(outlier.test.run.2$num.outliers);
```