---
title: "L-shaped selection from expression and methylation data"
author: "Berta Miro Cau and Alex Sánchez-Pla"
date: "`r Sys.Date()`"
output:
  rmarkdown::html_vignette: default
vignette: > 
  %\VignetteIndexEntry{L-shaped selection from methylation and expression}
  %\VignetteEncoding{UTF-8}  
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

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

# Introduction
This vignette provides examples of how to use the `LHeuristic` package, 
designed to identify genes with an L-shaped pattern in scatterplots 
comparing gene expression and methylation.

The package enables users to:

- Data Input

    - Load gene expression and methylation values from CSV files
    
- Data Preprocessing

    - Set parameters for data usage, choosing either all genes or those meeting 
    specific filters (e.g., $\mbox{corr}< 0$)
    - Select the L-shape detection method and configure its parameters
    
- Data Processing

    - Run computations to identify L-shaped genes
    - Generate scatterplots for all genes or filtered selections
    
- Data Output

    - Save the list of identified genes and scatterplots


# Input data

The data is a subset of methylation and expression arrays from the 
colon adenocarcinoma (COAD) dataset in The Cancer Genome Atlas (TCGA). 
The full dataset is available on the TCGA website.

Methylation and expression data are in two separate files. 
Columns represent patient samples, and rows represent genes.

```{r loadLibs}
library(Lheuristic)

# library(kableExtra)
# library(VennDiagram)

data("TCGAexpression")
data("TCGAmethylation")
```

## Using `lhCreateMAE` based on `MultiAssayExperiment` as data container

To effectively integrate and manage omics data,
we utilize the `MultiAssayExperiment` data structure.
This approach ensures consistency between datasets and 
facilitates downstream analysis.

The `lhCreateMAE` function constructs a MultiAssayExperiment
object using two datasets—typically methylation and
expression data—while enforcing consistency in sample
and feature names. It ensures:

The datasets have matching row and column names.

Experiments are labeled appropriately
(e.g., "methylation" and "expression").

Sample metadata can be incorporated
via the colData argument.

The function includes an internal validation step,
which verifies that both datasets are aligned correctly.
If discrepancies exist in row or column names, the
function will return an error, prompting users to
correct inconsistencies before proceeding.


```{r createMAE, warning=FALSE, message=FALSE}
library(MultiAssayExperiment)

colDat <- DataFrame(sampleID = colnames(TCGAmethylation))
rownames(colDat) <- colDat$sampleID


# Construct MultiAssayExperiment
mae1 <- lhCreateMAE(xDat = TCGAmethylation, yDat = TCGAexpression,
    xName = "methylation", yName = "expression", colData = colDat)

# Display dataset names
print(names(experiments(mae1)))
```

# Data analysis

The `LHeuristic` package implements the heuristic algorithm by 
Sanchez et al. (2019) to detect L-shaped scatterplots. However, since many 
researchers rely on negative correlation to identify genes "regulated by 
methylation," the package includes both methods for completeness.

Each method has various parameters optimized for selecting L-shaped patterns 
and detecting negative correlations in scatterplot data.

The functions also support parameter tuning, allowing for other forms of 
correlation selection.

## Correlation method

The *correlation method* is a simple approach that identifies genes with a
significant negative correlation between their expression and methylation 
values.

```{r corMethod}
cl <- correlationSelection(mae1, type = "Spearman",
    pValCutoff = 0.25, rCutoff = -0.5, adj = TRUE)

correlationL <- cl[!is.na(cl$SigNegCorr) & 
    cl$SigNegCorr, ]
correlationNoL <- cl[!is.na(cl$SigNegCorr)& 
    !cl$SigNegCorr, ]

message(
    "The number of genes selected 
    with the correlation method is: ", sum(correlationL$SigNegCorr),
    "\n"
)
```

We can observe the distribution of the correlation coefficients:

```{r corrDistri,fig.width = 6, fig.height=4}
d <- density(correlationL[, 1])
x2 <- data.frame(x = d$x, y = d$y)

library(ggplot2)
ggplot(x2, aes(x,y)) + geom_line() +
    labs(
    title = "Significant correlations in the TCGA dataset",
    x = "Correlation",
    y = "Density"
    )+
    theme_minimal()

```

The number of genes retained depends on the *p-value* and *r* cutoff.

The result is a table listing "selected" and "unselected" genes, with columns
for the *r coefficient*, *p-value*, *adjusted p-value*, 
*distance correlation*, and a Boolean indicating whether the gene was selected
based on a significantly negative correlation.

```{r showCorrs}
head(correlationL)
```

The selected genes can be plotted and saved as a PDF (leaving the file name 
empty displays the plot on screen). The figure below shows the first four 
selected genes.

```{r plotGenesFromCorr1, fig.height=6, fig.width=6}
genes2plot <- rownames(correlationL)[1:3]

plotGenesMat(mae1, geneNames = genes2plot,
    text4Title = correlationL[rownames(correlationL), ""],
    saveToPDF = FALSE
)
```


To plot scatterplots showing the relationship between methylation and expression
for a list of genes, use the function `plotGenesMat`.

## Heuristic method

The heuristic method identifies L-shaped scatterplots by overlaying a grid 
on the graph and defining specific cells that must contain a minimum 
(or maximum) percentage of points for the scatterplot to qualify as L-shaped.

It computes a score where points in the designated L region receive 
positive values, while points outside this region receive negative values.
Properly adjusting scores and weights results in positive scores for 
L-shaped scatterplots and negative scores for non-L-shaped ones.

This concept can be clarified with a "three-band rule":

1. Overlay a \(3 \times 3\) grid on the scatterplot.

2. Classify the scatterplot as **"L" or "non-L"** based on a few conditions    
    2.1. A *minimum* number of points must be in the upper-left (cell (1,1))
    and lower-right (cell (3,3)) corners.
    2.2. A *maximum* number of points must be in the upper-right (cell (1,3)), 
    as points here indicate hypermethylation and hyperexpression, which are 
    contrary to our goal.
    2.3. A minimum of points in cell (3,1) is usually *not required* 
    unless we strictly want an L-shape; we can also accept diagonals 
    that reflect a negative correlation.

3. Score points in each sub-grid as follows:
    3.1. Points in permitted regions (left margin: cells (1,1), (2,2), (3,1),
    (3,2), (3,3)) score positively if classified as L, or zero 
    if classified as non-L.
    3.2. Points in non-desired regions (outer band: cells (1,2), (1,3), 
    (2,3)) always score negatively.
    3.3. Some regions, like cell (2,2), may be neutral and not score.

4. *Tune scoring parameters either manually (based on experience and dataset*
*characteristics) or through cross-validation* 
(**if a set of positive and negative L-shaped genes is available**).

The scheme can be summarized using the following equation:


$$
S(X) = W_L \circ X \times \mathbf{1}_L(X) + W_{L^C} \circ X \times 
\mathbf{1}_{L^c}(X)
$$ 
where:


- \(X\) is the matrix of *counts*, representing the number of counts in each 
cell of the grid.
- \(W_L\) is the matrix of scores per cell and point if the scatterplot is 
classified as \(L\).
- \(W_{L^c}\) is the matrix of scores per cell and point if the scatterplot is 
classified as non-\(L\) (\(L^c\)). 
The symbol \(\circ\) denotes the Hadamard product of the two
matrices \(W_{L/L^c}\) (i.e., element-wise multiplication), 
and \(\mathbf{1}_{L/L^c}()\) is
the indicator function for \(L\) or \(L^c\).

The classification of the scatterplot as \(L\) or \(L^c\) can also be 
represented by the Hadamard product of three matrices:

\[
\mathbf{1}_L(X) = \bigwedge_{i,j} X \circ C \circ 
\left( mMP \times \sum_{i,j} x_{ij}
\right),
\]

where:
- \(X\) is the matrix of *counts*, indicating the number of counts 
in each grid cell.

- \(C\) is the matrix of conditions to verify if the scatterplot 
should be classified as \(L\).
- \(mMP\) is the matrix of minimum and maximum percentages of points 
required in each cell for classification as \(L\).
- \(\circ\) represents point-wise logical operations, allowing the product
of the three matrices to act as a logical operation.
- \(\bigwedge_{i,j}\) denotes a logical "AND" operation across all cells;
if all cells are TRUE, the result is assigned to \(L\); if any cell fails,
it is assigned to \(L^c\).


```{r setWeights}
sampleSize <- dim(mae1[[2]])[2]
numGenes <- dim(mae1[[2]])[1]

reqPercentages <- matrix(c(2, 20, 5, 5, 40, 20, 3, 3, 2), 
    nrow = 3, byrow = TRUE)
sum(reqPercentages)
(maxminCounts <- toReqMat(sampleSize, reqPercentages))

(theWeightMifL <- matrix(c(2, -2, -sampleSize / 5, 1, 0, -2, 1, 1, 2), 
    nrow = 3, byrow = TRUE))
(theWeightMifNonL <- matrix(c(0, -2, -sampleSize / 5, 0, 0, -2, 0, 0, 0),
    nrow = 3,
    byrow = TRUE
))

heur <- scoreGenesMat(mae1,
    aReqPercentsMat = reqPercentages, aWeightMifL = theWeightMifL,
    aWeightMifNonL = theWeightMifNonL
)

message("Number of scatterplots scored  : ", dim(heur)[1], "\n")
message("Number of L-shape scatterplots : ", sum(heur[, 1]), "\n")

heurL <- heur[heur$logicSc, ]
heurNoL <- heur[!heur$logicSc, ]
```

We can check the results in the following table, which includes a logical value
indicating whether the gene exhibits an L-shape based on our criteria, 
along with the *numerSc* score.


```{r showTrueLs}
knitr::kable(heurL)
```

Next, we can visualize the scatterplots of the selected genes or save them
on a PDF file.

```{r plotTrueLs, fig.height=6, fig.width=6}

genes2plot2 <-   rownames(heurL)[1:4] #  

plotGenesMat(mae1, geneNames = genes2plot2,
    fileName = NULL, text4Title = heurL[genes2plot2, "numeriSc"],
    saveToPDF = FALSE
)

```


In summary, the Heuristic method allows us to select genes with an L-shaped 
scatterplot as follows:

1. **Select datasets:** A pair of row-column matched matrices for expression and
methylation.

2. **Set parameters:**
    2.1. Grid definition
    2.2. Binary scoring
    2.3. Numerical scoring

3. **Score the selected data:** Return scores, the classification 
(L-Shape = TRUE / non-L-shape = FALSE), and plots for each gene.

# Comparison of Selected Gene Lists

After identifying genes using the two methods, we can create the intersection 
of all lists and visualize the results with a Venn Diagram.

We may choose genes selected by two or more methods to ensure greater 
consistency in the selection.


```{r checkCommon}
inCommonL <- intersect(rownames(correlationL), 
    rownames(heurL))
inCorrelationLOnly <- setdiff(rownames(correlationL),
    inCommonL)
inheurLLOnly <- setdiff(rownames(heurL), inCommonL)
```


We can also plot selected genes.
As an example we will plot the genes 1, 2, 3 and 5 by index.

```{r plotSelected, fig.height=6, fig.width=6}
par(mfrow = c(2, 2))
myGene1 <- inCommonL[1]
titleT <- paste(myGene1, "(May be GRM)")
plotGeneSel(mae1,myGene1,
    titleText = titleT,
    x1 = 1/3, x2 = 2/3)

myGene2 <- inCommonL[2]
titleT <- paste(myGene2, "(May be GRM)")
plotGeneSel(mae1, myGene2,
    titleText = titleT,
    x1 = 1/3, x2 = 2/3)


myGene3 <- inCommonL[3]
titleT <- paste(myGene3, "(May be GRM)")
plotGeneSel(mae1,myGene3,
    titleText = titleT,
    x1 = 1/3, x2 = 2/3)

myGene5 <- inCommonL[5]
titleT <- paste(myGene5, "(May be GRM)")
plotGeneSel(mae1, myGene5,
    titleText = titleT,
    x1 = 1/3, x2 = 2/3)
```

# References

Sanchez-Pla A, Miró B, Carmona F et al. 
A heuristic algorithm to select genes potentially regulated by methylation
[version 1;  not peer reviewed]. 
F1000Research 2019, 8:1017.
(https://doi.org/10.7490/f1000research.1116986.1)

# Session info

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


