1 Introduction

Identifying important transcription factor (TF) motifs, as shown in 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 monaLisa, we implement the randomized lasso stability selection proposed by Meinshausen and Bühlmann (2010) with the improved error bounds introduced by Shah and Samworth (2013). We have modified the stabs::glmnet.lasso function used by stabs::stabsel from the 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 (Meinshausen and Bühlmann 2010).

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.

2 Motif selection with Randomized Lasso Stability Selection

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.

2.1 Load packages

We start by loading the needed packages:

library(monaLisa)
library(JASPAR2020)
library(TFBSTools)
library(BSgenome.Mmusculus.UCSC.mm10)
library(Biostrings)
library(SummarizedExperiment)
library(ComplexHeatmap)
library(circlize)

2.2 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.

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

2.3 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.

# 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]
#>              
#>               NR3C2 Arnt LHX1 SNAI1 MAFG ZSCAN4
#>   peak_51663      0    0    0     0    0      0
#>   peak_57870      0    0    0     0    0      0
#>   peak_2986       2    0    0     0    0      0
#>   peak_124022     0    0    0     0    0      0
#>   peak_29925      0    2    0     1    0      0
#>   peak_95246      0    0    0     0    0      0

# remove TF motifs with 0 binding sites in all regions
zero_TF <- colSums(TFBSmatrix) == 0
sum(zero_TF)
#> [1] 2
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]
#>                fracGC     oeCpG NR3C2 Arnt LHX1 SNAI1
#> peak_51663  0.5155709 0.4079115     0    0    0     0
#> peak_57870  0.4963235 0.3048298     0    0    0     0
#> peak_2986   0.4008264 0.3103806     2    0    0     0
#> peak_124022 0.4572650 0.4429813     0    0    0     0
#> peak_29925  0.4675000 0.3495939     0    2    0     1
#> peak_95246  0.5144509 0.4020976     0    0    0     0

2.4 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.

# # 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
#> class: SummarizedExperiment 
#> dim: 10000 300 
#> metadata(12): stabsel.params.cutoff stabsel.params.selected ...
#>   stabsel.params.call randStabsel.params.weakness
#> assays(1): x
#> rownames(10000): peak_51663 peak_57870 ... peak_98880 peak_67984
#> rowData names(1): y
#> colnames(300): fracGC oeCpG ... CLOCK OLIG2
#> colData names(20): selProb selected ... regStep16 regStep17

# selected TFs
colnames(se)[se$selected]
#>  [1] "NKX2-5"       "GATA1::TAL1"  "HNF1B"        "HNF4A(var.2)" "Nr2f6"       
#>  [6] "ONECUT3"      "MYF5"         "THRB"         "ISL2"         "NR2C2"       
#> [11] "TEAD3"        "TEAD4"        "GATA3"        "RORA"         "NFIC"        
#> [16] "ZEB1"

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

plotStabilityPaths(se)