---
title: "Transcriptomic Classification using Mixture of Normal Distributions"
output:
  BiocStyle::html_document:
    toc: true
bibliography: splitTypeR.bibtex
vignette: >
  %\VignetteIndexEntry{Transcriptomic Classification using Mixture of Normal Distributions}
  %\VignettePackage{splitTypeR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r style, echo = FALSE, results = 'asis', message = FALSE}
BiocStyle::markdown()
suppressPackageStartupMessages({
    library(knitr)
    library(SummarizedExperiment)
})
```

<br />

**Package**: `r Rpackage("splitTypeR")`<br />
**Authors**: `r packageDescription("splitTypeR")[["Author"]]`<br />
**Version**: `r packageDescription("splitTypeR")$Version`<br />
**Compiled date**: `r Sys.Date()`<br />
**License**: `r packageDescription("splitTypeR")[["License"]]`<br />

# Licensing and citing

This package and the underlying __`r Rpackage("splitTypeR")`__ code are 
distributed under the Artistic license 2.0. You are free to use and 
redistribute this software.

If you use this package for a publication, we would ask you to cite the
following:

```{r citation, echo=TRUE}
## Citation information for the splitTypeR package
citation("splitTypeR")
```

<br />

# Introduction

Cancer classification, using RNA sequencing, is now part of modern precision 
oncology. For some cancers, RNA classification provides deep insights into 
tumor aggressiveness and growth rates. For example, breast cancer subtypes 
(like Luminal-A, Luminal-B, HER2-enriched, or Basal-like) help forecast 
patient outcomes [@Cascianelli2020]. The Luminal-A is the most 
common molecular subtype of breast cancer, making up to 60% of all breast 
cancer cases. It is also associated with good prognosis, while 
Luminal-B tumors have a more aggressive phenotype and a worse 
prognosis [@Yersal2014]. Once widely accepted by major oncology
organizations, gene signatures for RNA subtyping are used to classify cancer
subtypes, predict patient survival, and guide treatment.

From a research perspective, the classification of new samples using their 
transcriptional profiles can be done using different tools. For some cancers, 
machine learning models and dedicated software are available such as 
the Bioconductor Genefu package for breast cancer molecular 
subtyping [@Gendoo2016], and the CMScaller package for subtyping of 
colorectal cancer [@Eide2017].

In the absence of dedicated software, unsupervised clustering with a 
heatmap visualization to classify transcriptomic samples with a 
gene signature is often employed. This method is relatively simple and 
offers instant visual support. However, this method has a few drawbacks:

* _Lack of statistical metrics_: It doesn't 
provide a statistical probability score for a sample's class assignment; it 
mainly acts as a visual summary. 
* _Metric sensitivity_: Changing the clustering metric (e.g., from "Euclidean" to 
"Pearson" distance) can alter the resulting clusters. 
* _Outlier distortion_: Extreme outliers can squash the color scale, and 
artificially create skewed subsets.

The __`r Rpackage("splitTypeR")`__ package resolves these issues by providing 
an automated statistical framework to classify heterogeneous biological 
samples based on gene signature lists, effectively isolating the 
signature-positive samples. 

This classification method is designed for bulk transcriptomic datasets.

<br />

# Installation

As with any R package, the __`r Rpackage("splitTypeR")`__ package should first 
be loaded with the following command:

```{r loadingPackage, warning=FALSE, message=FALSE}
library(splitTypeR)
```

<br />

# Quick start

The __`r Rpackage("splitTypeR")`__ package includes a small demonstration 
dataset: two gene signatures
(`signaturesDemo`) and a matrix of normalized expected counts for 30
patients (`expNormalCountsDemo`). The whole subtyping method can be run
on it in a few lines:

```{r quickStart, echo=TRUE}
## Load the demo signatures and the demo normalized expected counts
data(signaturesDemo)
data(expNormalCountsDemo)

## Fix seed for reproducible results
set.seed(121)

## Run the subtyping method
## The normalized expected counts used as input are log2-transformed 
quickResult <- runSubtypingBimodal(geneLists=signaturesDemo,
    expectedCountsMatrixLog=log2(expNormalCountsDemo + 1),
    permRatio=0.75, permNbr=20, upscaleNbr=10)

## Display a summary of the results
quickResult
```

Each step of the method is detailed in the next section, followed by a
complete demonstration on a published dataset.

<br />

# Method

The following figure represents the statistical method employed to turn a noisy, 
single-number enrichment score into a confident binary classification of 
each sample. In summary, the method estimates the uncertainty around each
enrichment score, models the cohort as a mixture of two biological states, and 
then statistically tests each sample against that mixture.


```{r pipelineSchema, echo=FALSE, fig.align='center', out.width="100%", fig.cap="General workflow."}
knitr::include_graphics("enrichment_pipeline_300dpi_crop2.jpg")
```
<br />

The statistical method is split into 5 steps. Each of those steps is 
performed for one specific signature at a time. A signature is typically a 
predefined set of features (transcriptomic expression in this situation) that 
is associated with some biological state, pathway, or phenotype.


## Step 1. Enrichment score

The enrichment scores are first computed for every single sample 
with the Bioconductor `r BiocStyle::Biocpkg("GSVA")` 
package [@Hanzelmann2013]. Those scores represent the real score for each 
sample. Later steps treat this real score as the anchor point around which 
uncertainty gets modeled.

## Step 2. Permutation

Random subsets of samples, of fixed size, from the full dataset are 
repeatedly drawn. For each of these resampled subsets, the enrichment score 
for the samples present is recomputed.

This step enables a per-sample uncertainty estimation (standard deviation).

## Step 3. Upscaling 

For each sample, a fixed number of simulated enrichment scores are drawn 
from a normal distribution centered on the real score, using the standard 
deviation estimated from permutations.

This step enables the upscaling of the total number of values that will 
be used in the next step. 

## Step 4. Mixture model

The synthetic values across all samples are pooled to fit two overlapping 
normal distributions to them. The two normal distributions are extracted using 
a mixture model approach implemented in the `r BiocStyle::CRANpkg("mixtools")`
package [@JSSv032i06]. The assumption is that the whole cohort is a
blend of two underlying populations (a "low-mean" population and a
"high-mean" population).


## Step 5. Classification

Each sample is classified with a one-sided hypothesis test based on the 
"low-mean" distribution of the mixture model (from Step 4):

* **Null hypothesis (H₀)**: the sample belongs to the "low-mean" distribution
* **Alternative hypothesis (H₁)**: the sample belongs to the "high-mean" distribution

The 95th percentile of the "low-mean" distribution is used as the
classification threshold. Each sample's real enrichment score (from Step 1)
is compared to this threshold. A score above the threshold has less than a
5% probability of being produced by the "low-mean" distribution (one-sided
p-value < 0.05). In that case, H₀ is rejected.

```{r step5Schematic, echo=FALSE, message=FALSE, fig.width=7, fig.height=4, fig.align='center', fig.cap="Classification of the samples (schematic). The curve represents the \"low-mean\" distribution of the mixture model. The dashed line is the classification threshold (95th percentile) and the shaded area represents the upper 5% of the distribution. Sample A, with a score below the threshold, remains unclassified while sample B, with a score above the threshold, is assigned to the signature."}
library(ggplot2)

## Illustrative "low-mean" distribution (values chosen for the schematic)
lowMean <- -0.4
lowSD <- 0.2
threshold <- qnorm(p=0.95, mean=lowMean, sd=lowSD)

curve <- data.frame(score=seq(-1, 0.6, length.out=500))
curve$density <- dnorm(curve$score, mean=lowMean, sd=lowSD)

## Two example samples, one on each side of the threshold
examples <- data.frame(score=c(-0.3, 0.25),
    label=c("Sample A\nunclassified", "Sample B\nassigned to\nthe signature"),
    hjust=c(0.5, 0), nudgeX=c(0, -0.02))

ggplot(curve, aes(x=score, y=density)) +
    geom_area(data=curve[curve$score >= threshold, ], fill="firebrick",
        alpha=0.4) +
    geom_line(linewidth=1) +
    geom_vline(xintercept=threshold, linetype="dashed", color="gray30") +
    annotate("text", x=threshold + 0.02, y=max(curve$density) * 1.05,
        label="Threshold\n(95th percentile)", hjust=0, vjust=1, size=3.5) +
    annotate("text", x=threshold + 0.035, y=0.12, label="5%",
        color="firebrick4", size=3.5, fontface="bold") +
    geom_point(data=examples, aes(x=score, y=0), size=3.5, shape=17) +
    geom_text(data=examples, aes(x=score + nudgeX, y=0.1, label=label,
        hjust=hjust), vjust=0, size=3.2, lineheight=0.9) +
    labs(x="Enrichment score", y="Density",
        title="\"Low-mean\" distribution of the mixture model") +
    theme_bw()
```

When H₀ is rejected (p-value < 0.05), the sample is assigned to the tested 
signature. Otherwise, the sample is labeled as "unclassified" for that 
signature.

Since each signature is tested independently, a sample can be assigned 
to one signature, to several signatures, or to none of them. 

<br />


# Available signatures

The __`r Rpackage("splitTypeR")`__ package provides transcriptomic signatures.

The _getGeneSignaturesNames()_ function lists the names of all the 
available signatures.

```{r signaturesNames, echo=TRUE}
## List the names of the available signatures
getGeneSignaturesNames()
```

<br />

The _getGenesInSignature()_ function returns the genes of one signature, 
while the _getGeneSignatures()_ function returns a named list containing one 
or more signatures. When no name is given, _getGeneSignatures()_ returns all 
the available signatures.

```{r signaturesContent, echo=TRUE}
## Get the genes of the basal-like signature
getGenesInSignature(name="Tiriac_PDAC_PDO_Basal_like")

## Get all the available signatures and count the genes in each one
allSignatures <- getGeneSignatures()
lengths(allSignatures)
```

<br />

The current signatures are:

* **Tiriac_PDAC_PDO_Basal_like**: the basal-like signature 
(26 genes) defined for pancreatic ductal adenocarcinoma (PDAC) 
patient-derived organoids (PDOs) in Tiriac et al. [@Tiriac2018].
* **Tiriac_PDAC_PDO_Classical**: the classical signature 
(62 genes) defined for PDAC PDOs in Tiriac et al. [@Tiriac2018].
* **Moffitt_Basal_like**: the original basal-like signature 
(25 genes) defined for PDAC in Moffitt et al. [@Moffitt2015].
* **Moffitt_Classical**: the original classical signature 
(25 genes) defined for PDAC in Moffitt et al. [@Moffitt2015].
* **Miyabayashi_Fast**: the fast signature (231 genes), updated from the 
Ensembl IDs (GENCODE v48/org.Hs.eg.db), defined 
for PDAC xenografts in Miyabayashi et al. [@Miyabayashi2020].
* **Miyabayashi_Slow**: the slow signature (333 genes), updated from the 
Ensembl IDs (GENCODE v48/org.Hs.eg.db), defined 
for PDAC xenografts in Miyabayashi et al. [@Miyabayashi2020].

<br />

## Using your own signatures

The method is not limited to the signatures included in the 
__`r Rpackage("splitTypeR")`__
package. Any named `list` of at least two gene vectors can be passed to the 
_runSubtypingBimodal()_ function. The names of the list are used as the 
signature names in the results. The gene identifiers must use the same 
format as the row names of the expression matrix (e.g., gene names, Ensembl 
IDs). Genes absent from the expression matrix are ignored. It is recommended 
to use signatures with 10 genes or more to ensure statistical stability and 
reliability of the enrichment scores.

```{r ownSignatures, echo=TRUE, eval=FALSE}
## Example of custom signatures
mySignatures <- list(
    SignatureA=c("GENE1", "GENE2", "GENE3", "GENE4", "GENE5", "GENE6",
            "GENE7", "GENE8", "GENE9", "GENE10", "GENE11"),
    SignatureB=c("GENE4", "GENE9", "GENE10", "GENE12", "GENE13", "GENE15",
            "GENE16", "GENE17", "GENE18", "GENE19", "GENE20", "GENE21"))
```

<br />


# Demonstration

For this demonstration, we are going to use the published RNA-seq dataset of 
pancreatic ductal adenocarcinoma (PDAC) patient-derived organoids (PDOs) 
published in Tiriac et al. [@Tiriac2018]. The RNA-seq has been aligned with 
GENCODE v48 on GRCh38.p14 genome with STAR v2.7.9a [@Dobin2013]. Only one
PDO per patient was retained. The expected counts were processed with 
`r BiocStyle::Biocpkg("DESeq2")` v1.50.2 [@Love2014] on genes expressed in at 
least 20% of the PDOs. The normalized expression is available in a 
_SummarizedExperiment_ object. 

In PDAC, the basal-like and classical subtypes, as defined in Moffitt et 
al. [@Moffitt2015], are the two consistently employed molecular subtypes. The 
classical signature is associated with better prognosis and better response to 
standard chemotherapy, while the basal-like subtype is associated 
with worse overall prognosis [@Moffitt2015; @OKane2020]. These signatures were
later redefined for patient-derived organoids [@Tiriac2018].

<br />

## Prepare the signatures

The basal-like and classical gene signatures used in the demonstration are 
the ones defined for PDOs in Tiriac et al. [@Tiriac2018]. 

```{r geneSignatures, echo=TRUE}

## Required package
library(splitTypeR)

## Prepare the PDOs basal-like and classical signatures
signaturesList <- list()
signaturesList[["Basal-like"]] <- 
        getGenesInSignature("Tiriac_PDAC_PDO_Basal_like")
signaturesList[["Classical"]] <- 
        getGenesInSignature("Tiriac_PDAC_PDO_Classical")

```

<br />

## Prepare the log2-transformed normalized expected counts

A matrix of normalized expected counts is required. The rows are the genes, 
while the patients (samples) are the columns.

The row annotation must correspond to the signature format. As an example, if 
the signatures are composed of gene names, the row annotation for the 
normalized expected counts must also be gene names.

For GSVA analysis, the normalized expected counts should be log-transformed to 
enable the use of the Gaussian distribution function. The current 
behavior of GSVA is to automatically choose one of the available distribution 
functions. 


```{r expCount, echo=TRUE}

## Required package
library(SummarizedExperiment)

## Load the SummarizedExperiment object
path <- system.file("extdata", "Tiriac_NormExpCounts.RDS", 
                        package="splitTypeR")
se <- readRDS(path)

## Extract the normalized counts and assign the gene names as row names
## The row names should correspond to the information present in the 
## gene signatures. In this case, the gene names are used. 
normCounts <- assays(se)[["normalizedCounts"]]
rownames(normCounts) <- rowRanges(se)$gene_name

## Log2-transform the data, adding a small pseudo-count to 
## prevent log2(0) errors 
log2normCounts <- log2(normCounts + 1)

```

<br />

## Run the subtyping method

The _runSubtypingBimodal_ function runs all the steps. Three parameters
control the permutation and upscaling steps:

| Parameter | Description | Default | Trade-off |
|:----------|:------------|:-------:|:----------|
| `permRatio` | Proportion of the samples randomly selected in each permutation (Step 2). Must be between 0 and 1, and must not select all the samples. | 0.75 | A lower value creates more variation between permutations, but each GSVA score is computed on fewer samples. |
| `permNbr` | Number of permutations (Step 2). Must be 5 or higher. | 20 | More permutations give a more stable standard deviation for each sample, but take longer to run since GSVA is run once per permutation. |
| `upscaleNbr` | Number of values drawn from the normal distribution of each sample (Step 3). Must be 2 or higher. | 10 | More values give more data to fit the mixture model (Step 4), with little impact on the run time. |

<br />

Each sample is present, on average, in `permRatio` × `permNbr` permutations
(15 with the default values), which is the number of values used to estimate
its standard deviation.

<br />

```{r subtyping, echo=TRUE}

## Fix seed for reproducible results
set.seed(81)

## Run the subtyping method using the log2-transformed data
result <- runSubtypingBimodal(geneLists=signaturesList, 
    expectedCountsMatrixLog=log2normCounts,
    permRatio=0.75,
    permNbr=20,
    upscaleNbr=12)

## Summary of the classification results
result

```

<br />

The resulting object contains information from each step of the process. 
The information can be accessed using the appropriate getter function.

<br />

```{r signaturesGetter, echo=TRUE}
## Get the list of genes in the classical signature
signatures(result)[["Classical"]]
```

<br />

```{r gsvaResultsGetter, echo=TRUE}
## Get the GSVA scores for all samples in each signature
gsvaResults(result)
```

<br />

```{r permutationGetter, echo=TRUE}
## Get the GSVA scores for the first 5 samples in the first 6 permutations 
## (columns) done for the classical signature
permutations(result)[['Classical']][seq_len(5), seq_len(6)]
```

<br />

```{r standardDeviationGetter, echo=TRUE}
## Get the standard deviation of the GSVA scores for all samples for the 
## classical signature
## The standard deviation is inferred from the permutation step
standardDeviation(result)$Classical
```

<br />

```{r upscalingGetter, echo=TRUE}
## Get the first 6 upscaling GSVA scores (columns) for the first 5 samples 
## for the classical signature
upscaling(result)[['Classical']][seq_len(5), seq_len(6)]
```

<br />

```{r modelGetter, echo=TRUE}
## Get the mixture model for the classical signature
## This is a mixEM object
classicalModel <- model(result)[['Classical']]

## Get the two distribution means
classicalModel$mu

## Get the two distribution standard deviations
classicalModel$sigma

## Get the mixing proportions (the fraction of the values assigned to each
## distribution)
classicalModel$lambda
```

<br />

```{r thresholdGetter, echo=TRUE}
## Get the classification threshold of each signature
## The threshold is the 95th percentile of the "low-mean" distribution
threshold(result)
```

<br />

```{r classificationGetter, echo=TRUE}
## Get the final classification of the first samples 
## for the classical signature
head(classification(result)[['Classical']])
```

<br />

## Visual representation of the results at the permutation step

The _plotPermutationSamplesDistribution()_ function enables the visualization 
of the GSVA results obtained during the permutation step for a specific 
signature and specific samples. 

```{r graphPermutation, echo=TRUE}

## The GSVA results obtained for three samples during the permutation step 
## for the classical signature
## The seed is fixed to ensure reproducibility
plotPermutationSamplesDistribution(x=result, signature="Classical", 
    pointColor="darkorange", alpha=1, samples=c("hF23", "hT82", "hT89"), 
    seed=121)

```

<br />

## Visual representation of the results at the upscaling step

From the permutation step, the standard deviation is calculated for each 
sample. The upscaling is done with the real GSVA score as mean and the 
calculated standard deviation.

The _plotPermutationSamplesNormal()_ function shows the GSVA scores sampled 
during the upscaling step. The mean is shown as a vertical dotted line, while 
the normal distribution is shown as a dashed line. The sampled GSVA scores 
are shown as dots in this graph.

```{r graphUp, echo=TRUE}

## Fix seed for reproducibility
set.seed(121)

## The GSVA scores sampled for three samples during the upscaling step 
## for the classical signature
plotPermutationSamplesNormal(result, signature="Classical", 
    samples=c("hF23", "hT82", "hT89"), pointColor="darkviolet")
```

<br />

## Visual representation of the results at the mixture model step

The _plotDensityWithUpscalingSamplesBimodal()_ function enables the 
visualization of the two normal distributions obtained for a specific 
signature. In this graph, the GSVA values are those generated during the 
upscaling step.

```{r graphDensity, echo=TRUE}

plotDensityWithUpscalingSamplesBimodal(breaks = 26L,
    signature="Basal-like", result)

plotDensityWithUpscalingSamplesBimodal(breaks = 26L,
    signature="Classical", result)

```

<br />

## Classification of the samples

Let's compare the classification obtained here with the one from 
Tiriac et al. [@Tiriac2018].


```{r classificationTiriac, echo=TRUE}

## Assignments from Tiriac et al.
pathAssign <- system.file("extdata", "Moffitt_Tiriac_2018_Assignation.txt",
                        package="splitTypeR")
tiriacAssign <- read.table(pathAssign, sep="\t", header=TRUE)

## Assignments from the mixture model are extracted from the result
## By default, the samples are considered unclassified
## The samples are classified only if they are assigned to one of the signatures
## When assigned to both signatures, the samples remain unclassified
assignation <- data.frame(
        Samples=rownames(classification(result)[['Classical']]), 
        newAssignation="unclassified")
assignation$Classical <- classification(result)[['Classical']]$classification
assignation$Basal <- classification(result)[['Basal-like']]$classification
assignation$newAssignation[which(assignation$Classical == "Classical" & 
                    assignation$Basal == "unclassified")] <- "Classical"
assignation$newAssignation[which(assignation$Basal == "Basal-like" & 
                    assignation$Classical == "unclassified")] <- "Basal-like"

## Merge the current assignments with the Tiriac et al. assignments into 
## one data frame
assignation <- merge(assignation[, c("Samples", "newAssignation")], 
                        tiriacAssign, by.x="Samples", 
                        by.y="Line", all.y=FALSE, all.x=FALSE)

table(assignation$newAssignation, assignation$Moffitt)
```

<br />

About
`r round(sum(assignation$newAssignation == assignation$Moffitt)/nrow(assignation)*100)`%
of
the samples (`r sum(assignation$newAssignation == assignation$Moffitt)` out of
`r nrow(assignation)`) are assigned the same subtype.

The following graph shows the GSVA scores of each sample for both
signatures. The dashed lines represent the classification threshold of each
signature (the 95th percentile of the "low-mean" distribution). A sample
located to the right of the vertical line is assigned to the classical
signature, while a sample located above the horizontal line is assigned to
the basal-like signature. The samples for which the classification differs
from Tiriac et al. [@Tiriac2018] are labeled.

```{r compared, echo=TRUE, fig.width=7, fig.height=5.5}

## Required package
library(ggplot2)

## Get the GSVA score of each sample for both signatures
gsva <- data.frame(
    Samples=rownames(classification(result)[['Classical']]),
    Classical=classification(result)[['Classical']]$GSVA_score,
    Basal=classification(result)[['Basal-like']]$GSVA_score)

## Add the assignments and identify the samples with a different
## classification than Tiriac et al.
gsva <- merge(gsva, assignation, by="Samples")
gsva$Agreement <- ifelse(gsva$newAssignation == gsva$Moffitt,
    "Same classification", "Different classification")

## Classification threshold of each signature
thresholdClassical <- threshold(result)[['Classical']]
thresholdBasal <- threshold(result)[['Basal-like']]

## Scatter plot of the GSVA scores with the classification thresholds
## The vertical dashed line is the classical threshold and the horizontal
## dashed line is the basal-like threshold
ggplot(gsva, aes(x=Classical, y=Basal)) +
    geom_vline(xintercept=thresholdClassical, linetype="dashed",
        color="gray40") +
    geom_hline(yintercept=thresholdBasal, linetype="dashed",
        color="gray40") +
    geom_point(aes(color=Agreement, shape=Moffitt), size=3) +
    geom_text(data=gsva[gsva$Agreement == "Different classification", ],
        aes(label=Samples), vjust=-1, size=3.5) +
    scale_color_manual(values=c("Same classification"="gray60",
        "Different classification"="firebrick")) +
    labs(x="GSVA score - Classical signature",
        y="GSVA score - Basal-like signature",
        color="Compared to Tiriac et al.", shape="Tiriac et al.") +
    theme_bw()

```

<br />

Most samples that are not assigned to the same subtype have intermediate
scores and are located where the classical and basal-like samples meet.
Sample hF50 is located in the region where both signatures are assigned and
therefore remains unclassified.

<br />


# Acknowledgment

D.A. Tuveson is a distinguished scholar of the Lustgarten Foundation and 
Director of the Lustgarten Foundation-designated Laboratory of Pancreatic 
Cancer Research. This project is supported by the Lustgarten Foundation.

Claude Sonnet 5 was used for improving the workflow figure. It was also used 
for proofreading and copyediting sections of this text.

Claude Opus 5.5, through Claude Code, was used to assist with code 
maintenance, including renaming the demonstration signatures, updating the 
associated documentation and unit tests, removing unused internal data, and 
adding the Quick start section to this vignette. All changes were reviewed 
and validated by the package maintainer.

<br />

# Session info

Here is the output of `sessionInfo()` on the system on which this document was 
compiled:

```{r sessionInfo, echo=FALSE}
sessionInfo()
```

<br />