---
title: "Regression Based Approach for Motif Selection"
author: "Dania Machlab, Lukas Burger, Charlotte Soneson, Michael Stadler"
date: "`r Sys.Date()`"
bibliography: monaLisa-refs.bib
output: 
  BiocStyle::html_document:
    toc_float: true
vignette: >
  %\VignetteIndexEntry{selecting_motifs_with_randLassoStabSel}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 7,
    fig.height = 5,
    out.width = "80%",
    fig.align = "center",
    crop = NULL # suppress "The magick package is required to crop" issue
)
library(BiocStyle)
```

<img src="monaLisa_logo_v1.png" style="display: block; margin-left: auto;
margin-right: auto; width: 50%; border: 0" />

# Introduction

Identifying important transcription factor (TF) motifs, as shown in 
`r Biocpkg("monaLisa", vignette="monaLisa.html", label="the main vignette")`, could also be done using a regression based approach, 
where motifs are selected and have to compete against each other for selection. 
In this framework, the response vector can be the observed experimental measure 
of interest, e.g. log-fold changes of accessibility for a set of regions, and 
the predictors consist of the TF motif hits across those regions. In 
`r Biocpkg("monaLisa")`, we implement the randomized lasso stability selection 
proposed by @StabSel with the improved error bounds introduced by @compStabs. 
We have modified the `stabs::glmnet.lasso` function used by `stabs::stabsel` 
from the `r CRANpkg("stabs")` package to implement the randomized lasso.

Lasso stability selection implements the lasso regression multiple times on 
subsamples of the data, and returns a selection probability for each predictor 
(number of times selected divided by number of regressions done). With the 
randomized lasso, a weakness parameter is additionally used  to vary the lasso 
penalty term $\lambda$ to a randomly chosen value between 
[$\lambda$, $\lambda$/weakness] for each predictor. Although the main appeal of 
stability selection is in cases where the number of predictors exceeds the 
number of observations, it also performs better variable selection with noisy 
data [@StabSel].

With this approach, TF motifs compete against each other to explain the 
response vector, and we can also include additional predictors like GC content 
to compete against the TF motifs for selection. This is especially useful if 
the response is biased by sequence composition, for example if regions with 
higher GC content tend to have higher response values.


# Motif selection with Randomized Lasso Stability Selection{#stabsel}

In the example below, we select for TF motifs explaining log-fold changes in 
chromatin accessibility (ATAC-seq) across the enhancers between mouse liver and 
lung tissue at P0, but this can be applied to other data types as well 
(ChIP-seq, RNA-seq, methylation etc.). Positive log2-fold changes indicate 
more accessibility in the liver tissue, whereas negative values indicate 
more accessibility in the lung tissue.

## Load packages

We start by loading the needed packages: 

```{r libs, message=FALSE}
library(monaLisa)
library(JASPAR2020)
library(TFBSTools)
library(BSgenome.Mmusculus.UCSC.mm10)
library(Biostrings)
library(SummarizedExperiment)
library(ComplexHeatmap)
library(circlize)
```

## Load dataset

In our example dataset, we have quantified ATAC-seq reads on enhancers in 
mouse P0 lung and liver tissues. The log2-fold change (our response vector in 
this example) is for liver vs lung chromatin accessibility. We are using a set 
of 10,000 randomly sampled enhancers to illustrate how randomized lasso 
stability selection can be used to select TFs.

```{r loadData}
# load GRanges object with logFC and peaks
gr_path <- system.file("extdata", "atac_liver_vs_lung.rds", 
                       package = "monaLisa")
gr <- readRDS(gr_path)
```


## Get TFBS per motif and peak

We will now construct the transcription factor binding site (TFBS) matrix for 
known motifs (from a database like JASPAR) in peak regions. We use the 
`monaLisa::findMotifHits` function to scan for TF motif hits. This matrix will
be the predictor matrix in our regression. This step may take a while, and it 
may be useful to parallelize it using the `BPPARAM` argument (e.g. to run on 
`n` parallel threads using the multi-core backend, you can use: 
`findMotifHits(..., BPPARAM = BiocParallel::MulticoreParam(n))`).

As mentioned, this framework offers the flexibility to add additional 
predictors to compete against the TF motifs for selection. Here, we add the 
fraction of G+C and CpG observed/expected ratio as predictors, to ensure that 
selected TF motifs are not just detecting a simple trend in GC or CpG 
composition.

```{r predictorMatrix, warning=FALSE}
# get PFMs (vertebrate TFs from Jaspar)
pfms <- getMatrixSet(JASPAR2020, list(matrixtype = "PFM", 
                                      tax_group = "vertebrates"))

# randomly sample 300 PFMs for illustration purposes (for quick runtime)
set.seed(4563)
pfms <- pfms[sample(length(pfms), size = 300)]

# convert PFMs to PWMs
pwms <- toPWM(pfms)

# get TFBS on given GRanges (peaks)
# suppress warnings generated by matchPWM due to the presence of Ns 
# in the sequences
peakSeq <- getSeq(BSgenome.Mmusculus.UCSC.mm10, gr)
suppressWarnings({
  hits <- findMotifHits(query = pwms, subject = peakSeq, min.score = 10.0,
                        BPPARAM = BiocParallel::SerialParam())
})

# get TFBS matrix
TFBSmatrix <- unclass(table(factor(seqnames(hits), levels = seqlevels(hits)),
                            factor(hits$pwmname, levels = name(pwms))))
TFBSmatrix[1:6, 1:6]

# remove TF motifs with 0 binding sites in all regions
zero_TF <- colSums(TFBSmatrix) == 0
sum(zero_TF)
TFBSmatrix <- TFBSmatrix[, !zero_TF]

# calculate G+C and CpG obs/expected
fMono <- oligonucleotideFrequency(peakSeq, width = 1L, as.prob = TRUE)
fDi <- oligonucleotideFrequency(peakSeq, width = 2L, as.prob = TRUE)
fracGC <- fMono[, "G"] + fMono[, "C"]
oeCpG <- (fDi[, "CG"] + 0.01) / (fMono[, "G"] * fMono[, "C"] + 0.01)

# add GC and oeCpG to predictor matrix
TFBSmatrix <- cbind(fracGC, oeCpG, TFBSmatrix)
TFBSmatrix[1:6, 1:6]
```


## Identify important TFs

We can now run randomized lasso stability selection to identify TFs that are 
likely to explain the log-fold changes in accessibility. 

```{r stabSelTFs}
# # randLassoStabSel() is stochastic, so we set a seed to reproduce 
# # ... a parallel run
# RNGkind("L'Ecuyer-CMRG")
# set.seed(123)
# se <- randLassoStabSel(x = TFBSmatrix, y = gr$logFC_liver_vs_lung, 
#                        cutoff = 0.8, mc.preschedule = TRUE, 
#                        mc.set.seed = TRUE, mc.cores = 2L)

# if not running in parallel mode, it is enough to use set.seed() before 
# ... using the function to ensure reproducibility (with 1 core)
set.seed(123)
se <- randLassoStabSel(x = TFBSmatrix, y = gr$logFC_liver_vs_lung, 
                       cutoff = 0.8)
se

# selected TFs
colnames(se)[se$selected]
```


The stability paths visualize how predictors get selected, decreasing 
regularization stringency (from left to right): 

```{r plotStabilityPaths}
plotStabilityPaths(se)
```

Each line corresponds to a predictor, and we can see the selection 
probabilities as a function of the regularization steps, corresponding to 
decreasing values for the lambda regularization parameter in lasso. The 
predictor (TF motif) selection happens at the last step, given the specified 
minimum probability.

We can also visualize the selection probabilities of the selected TF motifs, 
optionally multiplied by the sign of the correlation to the response vector, 
to know how the TF relates to the change of accessibility (`directional` 
parameter). Note that although one can vary the `selProbMinPlot` argument 
which sets the selection probability cutoff, it is best to re-run randomized 
lasso stability selection with the new cutoff, as this influences other 
parameter choices the model uses internally. See @StabSel for more details.

```{r plotSelProbs, fig.width=15, fig.height=9}
plotSelectionProb(se, directional = TRUE)
```

Next, we visualize the correlation structure of the TFBS matrix and selected TF 
motifs. While the colinearity of predictors is generally an issue in 
regression-based approaches where variables are selected, randomized lasso 
stability selection generally does better at co-selecting correlated 
predictors. In practice, we still see it select predictors with correlations as 
high as 0.9. However, it is good to keep in mind that this can be an issue,
and that predictors that are extremely correlated with each other might not 
end up being co-selected. In a bad case example, if we have two identical 
predictors, neither of them could end up being selected if half of the times 
the regressions are done, one would be randomly selected, and the other one 
in the other half, giving each a selection probability of 0.5. If many cases 
of high correlations exist or this is a concern, one may consider selecting a 
representative set of predictors to use. 

```{r TFBScor_selected, fig.width=10, fig.height=7}
# subset the selected TFs
sel <- colnames(se)[se$selected]
se_sub <- se[, sel]

# exclude oeCpG and fracGC
excl <- colnames(se_sub) %in% c("oeCpG", "fracGC")
se_sub <- se_sub[, !excl]

# correlation matrix 
TFBSmatrixCorSel <- cor(TFBSmatrix[, colnames(se_sub)], method = "pearson")

# heatmap
pfmsSel <- pfms[match(colnames(TFBSmatrixCorSel), name(pfms))]
maxwidth <- max(sapply(TFBSTools::Matrix(pfmsSel), ncol))
seqlogoGrobs <- lapply(pfmsSel, seqLogoGrob, xmax = maxwidth)

hmSeqlogo <- rowAnnotation(logo = annoSeqlogo(seqlogoGrobs, which = "row"),
                           annotation_width = unit(2, "inch"), 
                           show_annotation_name = FALSE
)

colAnn <- HeatmapAnnotation(AUC = se_sub$selAUC, selProb = se_sub$selProb,
                            show_legend = TRUE, 
                            show_annotation_name = TRUE,
                            col = list(
                              AUC = colorRamp2(c(0, 1), 
                                               c("white", "brown")),
                              selProb = colorRamp2(c(0, 1), 
                                                   c("white", "steelblue")))
)

Heatmap(TFBSmatrixCorSel, 
        show_row_names = TRUE, 
        show_column_names = TRUE, 
        name = "Pear. Cor.", column_title = "Selected TFs",
        col = colorRamp2(c(-1, 0, 1), c("blue", "white", "red")), 
        right_annotation = hmSeqlogo,
        top_annotation = colAnn)
```


```{r TFBScor_all}
# exclude oeCpG and fracGC
excl <- colnames(se) %in% c("oeCpG", "fracGC")
se_sub <- se[, !excl]

# correlation matrix 
TFBSmatrixCorSel <- cor(TFBSmatrix[, colnames(se_sub)], method = "pearson")

# heatmap annotations
ann <- HeatmapAnnotation(
    selected = as.character(se_sub$selected), 
    AUC = se_sub$selAUC, 
    selProb = se_sub$selProb, 
    col = list(selected = c("TRUE"="forestgreen", 
                            "FALSE"="grey"), 
               AUC = colorRamp2(c(0, 1), 
                                c("white", "brown")), 
               selProb = colorRamp2(c(0, 1), 
                                    c("white", "steelblue")))
)

# heatmap
Heatmap(TFBSmatrixCorSel, show_row_names = FALSE, show_column_names = FALSE, 
        name = "Pear. Cor.", column_title = "All TFs",
        col = colorRamp2(c(-1, 0, 1), c("blue", "white", "red")), 
        top_annotation = ann)
```

\  

We can examine the peaks that have hits for a selected TF motif of interest, 
ordered by absolute accessibility changes.

```{r topPeaks}
TF <- sel[2]
TF

i <- which(assay(se, "x")[, TF] > 0) # peaks that contain TF hits...
nm <- names(sort(abs(gr$logFC_liver_vs_lung[i]), 
                 decreasing = TRUE)) # ... order by |logFC|

head(gr[nm])
```


# Session info and logo
The monaLisa logo uses a drawing that was obtained from 
http://vectorish.com/lisa-simpson.html
under the Creative Commons attribution - non-commercial 3.0 license: 
https://creativecommons.org/licenses/by-nc/3.0/.

This vignette was built using:  
```{r, session}
sessionInfo()
```

# References
