Neural Networks beyond explainability: Selective inference for sequence motifs
Authors: Antoine Villié
Philippe Veber
Yohann De Castro
Laurent Jacob
Venue: Accepted by TMLR ()

Abstract: Over the past decade, neural networks have been successful at making predictions from biological sequences, especially in the context of regulatory genomics. As in other fields of deep learning, tools have been devised to extract features such as sequence motifs that can explain the predictions made by a trained network. Here we intend to go beyond explainable machine learning and introduce SEISM, a selective inference procedure to test the association between these extracted features and the predicted phenotype. In particular, we discuss how training a one-layer convolutional network is formally equivalent to selecting motifs maximizing some association score. We adapt existing sampling-based selective inference procedures by quantizing this selection over an infinite set to a large but finite grid. Finally,we show that sampling under a specific choice of parameters is sufficient to characterize the composite null hypothesis typically used for selective inference - a result that goes well beyond our particular framework. We illustrate the behavior of our method in terms of calibration, power and speed and discuss its power/speed trade-off with a simpler data-split strategy. SEISM paves the way to an easier analysis of neural networks used in regulatory genomics, and to more powerful methods for genome wide association studies (GWAS).

Published in Transactions on Machine Learning Research (06/2023)
Neural Networks beyond explainability: Selective inference for sequence motifs
Antoine Villié
antoine, villie@unio-lyon. Université de Lyon, Université Lyon 1, CNRS, VetAgro Sup, Laboratoire de Biométrie et Biologie Evolutive, UMR5558, Villeurbonne, France
phillppe..eke wtw.:pe..eeeerruni---eee Wwnio-lyonl..f Université de Lyon, Université Lyon 1, CNRS, VetAgro Sup, Laboratoire de Biométrie Biologic Evolutive, UMR5558, Villeurbonne, France
Philippe Veber
Yohann De Castro
yoharn. he-.:am..dd.-ustroeeeoo....oo
Institut Camille Jordan, École Centrale Lyon, CNRS UMR 5208 Instifut universitaire de Prance (TUF)
Laurent Jacob
Sturent.jacobeeny
Sorbonne Université, CNRS, IBPS, Laboratory of Computational and Quantitative Biology (LCQB), UMR 7238, Paris 75005, France
Reviewed on OpenReview: https: // openreview. net/ forum? id= addEHTSngg
Abstract
Over the past decade, neural networks have been successful at making predictions from biological sequences. As in other fields of deep learning, tools have been devised to extract features such as sequence motifs that can explain the predictions made by trained net- work. Here we intend to go beyond explainable machine learning and introduce SEISM, a selective inference procedure to test the association between these extracted features and the predicted phenotype. In particular, we discuss how training a onc-layer convolutional network is formally equivalent to selecting motifs maximizing some association score. We adapt existing sampling- based selective inference procedures by quantizing this selection over an infinite set to large but finite grid. Finally, we show that sampling under specific choice of parameters is sufficient to characterize the composite mull hypothesis typically used for selective inference- a result that goes well beyond our particular framework. We illus- trate the behavior of our method in terms of calibration, power and speed and discuss its power power/speed trade-off with a simpler data-split strategy. SEISM paves the way an casier analysis of neural networks used in regulatory genomics, and to more powerful methods for genome wide association studies (GWAS).
1 Introduction
In the recent years, neural networks have been successfully used for making predictions from biological se- quences. In particular, they have brought significant improvements in regulatory genomies, e.g. to predict cell-type specific transcription factor binding, gene expression, chromatin accessibility or histone modifica- tions from DNA sequence (Zhou & Troyanskaya, 2015; Kelley et al., 2018; Avsec et al., 202 These tasks are expected to be good proxy for predicting the functional effect of non-coding variants, and help us in turn make better sense of the observed human genetic variation and its effect on various phenotypical traits including diseases. Most successful models have used convolutional neural networks (CNNs, LeCun & Bengio, 1998) and more recent approaches have explored self-attention mechanisms (Vaswani al., 2017). These models have been trained from experimental data obtained from ChIP-seq. ATAC-seq. DNaso-soq, CAGE assays, that provide examples where both the DNA sequence and the outcome of interest are known.

---

Published in Transactions on Machine Learning Research (06/2023)
A commonly outlined limitation of neural networks is their lack of explainability or black box aspect, i.e., the contrast. between their excellent prediction accuracy and the possibility to explain these in intuitive mechanistic terms (Ras et al., 2022; Molnar, 2022). Elementary one-layer CNNs don't face this issue, as their trained filters have a straightforward interpretation as position weight matrices (PWMs, Harr et al., 1983; Schneider & Stephens, 1990), historical and basic element of regulatory genomics. Nonetheless, these simple models are notoriously too simple to capture the complexity of the regulatory code which requires to account not only for individual motif presence but for their long range sequence context and mutual interactions (Avsec et al., 2021b). Multi-layer CNNs and self-attention mechanisms model this additional complexity but are less straightforward to interpret. Tools inspired from the explainable deep learning literature have been adapted to extract features beyond PWMs and one-layer CNNs to explain the predicted regulatory behavior (Novakóvsky et al., 2022). It is therefore often possible lo explain the predictions of a trained neural network for biological sequences, either directly through estimates of its parameters or through features extractod post hoc.
Unfortunately, finding features somewhat associated to an outcome is often not enough, as an observed non- zero association can be spurious. In experimental science, it is actually common to quantify the significance of this association, e.g., by testing the hypothesis that it is zero. Genome wide association studies (GWAS, Visscher et al., 2017) for example find genetic variants correlated with a trait by building a linear model explaining this trait by each variant and testing the hypothesis that the weight is zero. Statistical significance has its own limitations (Wasserstein relevant features. Quantifying the significance of associations between interpretable features and predicted outcome is equally important in the context of neural networks but has received little attention to our knowledge.
Lazar, 2016), but often provides an intuitive scale for identifying
Empárical mull distribution with data split Oz
Motifs space meshing
Motifs selection (NN training)
Selection event definition and sampling Motif epaie
Empirical conditional mull distribution p. p-value computation
Input dataset
Phenotopee y' 3.2
V- Nini
Sesuenore AGCA...
CE CAPIACEA
CGTA..C
0.1
Pacnotrpes
ATGC...A
6.7
.... 2.3/
cheAliTAA
CAGA..T
A(A (
o in Al
selection event QERECtAtiON SOOTE
*solo
Amas
Figure 1: Overview of our SEISM procedure. (a) The input is set of sequences and corresponding pheno- types in some space y (b) It trains a convolutional neural networks to predict phenotype from sequences, which leads to the selection of sequence motifs. (c) Then SEISM partitions the space of motifs to quantize the selection. The selection event is the set of phenotype vectors that would lead to selecting an element in the same mesh. (d) Using a sampling strategy, SEISM builds a null distribution for the test statistic, conditional to the selection event. The p-values associated with a selected motif is the quantile of its score under this distribution.
Here, we set out to go beyond explainable machine learning by introducing SElective Inference for Sequence Motifs (SEISM), depicted in Figure 1, valid statistical inference procedure for these features. In order to do so, cast commonly used CNNs in feature selection framework, and show that it achieves similar selection performances as existing bioinformatics algorithms on de novo motifs discovery tasks. This selection needs to be accounted for when iesting the association of the features with the predicted trait. This problem has been discussed and addressed in the growing literature on post-selection inference over the past few years, using e.g. data-split (Wasserman & Roeder, 2009) or selective inference strategies (Taylor & Tibshirani, 2015; Reid et al., 2018; Slim et al., 2019). The former split the data into two parts, performing selection on one and inference on the other. They produce valid inference but necessarily result in reduction of the sample
2

---

Published in Transactions on Machine Learning Research (06/2023)
size, which is unsatisfying when the original sample size is limited. By contrast, the latter condition the mull distribution on the selecting event, which generally provides more power but can prove more computationally intensive.
Our contributions are as follows:
We formally cast one-layer CNNs into motif discovery tool, reaching similar performances de-novo motils discovery tools from the bio-informatics literature (Section 3).
We define post-selection inference framework for the features selected by the neural network, using either data-split or selective inference (Section 4), cach being more appropriate in given sample size regime.
Both strategies require sampling under a normal mull hypothosis which is composite- mean vectors define the same null-and depends on an unknown parameters. We provide invariance results suggesting a practical procedure that works around these issues (Section 4.6). To our knowledge, they were blind spot in sampling-base basodd post-sclection inference approaches beyond our specific context.
Existing selective inference methods are only defined for selections over a finite set. We work around this issue by quantizing our selection to very large but finite space, making it amenable to existing sampling strategies. We show that the resulting procedure is well calibrated.
We provide PyTorch implementation of SEISM at:
https://gitlab. in2p3 .ntp3.fr/amtei....e.e.e. httt:/attiieeeiiieeeeeeeeeee wwwlieieseiisseeeeeeeeeeeeeeeeee
In this paper, we restrict our presentation to simple one-layer CNNs and sequence motifs. The procedure we introduce here, however, is not limited to this framework. It can be applied to more expressive features proposed in the explainable machine learning literature, but may require some further work depending on the feature considered.
2 A short overview of our SEISM procedure
SEISM aims to detect sequence motils associated with a biological outcome, and to test the statistical significance of this association. Here we briefly describe the selective inference version of SEISM in order to give the reader an overview of the procedure. It is summarized in Algorithm 1, and more details will be given in the following sections.
(i) SEISM takes input biological sequences X associated with phenotype y. The user must also specify the number of motifs to find, as well as parameter controlling the meshing of the motif space, that is the precision with which the found motifs will be tested.
(ii) The motif selection step corresponds to the maximisation of a so-called association score s(,,), +), which depends on the phenotype and on the motifs = through their activation patterns in the biological sequences x This step is formally equivalent to training a one hidden layer CNN. We implement a greedy procedure, optimizing each new filter over the residuals of the previously entered ones, using a gradient descent method initialized at the k-mer with the best score. To that end, we enumerates the k-mers contained in X using the DSK software (Rizk et al., 2013) and compute their scores s(;)))
(iii) SEISM splits the set of sequence motifs into meshes according to the input parameter. This step leads to the definition of set of null hypotheses and of selection event E, f.c. the set of outcomes y that would have led to the selection of motifs within the same meshes as the ones selected in (ii), namely M.,). Formally, the selection event reads
the sequence of meshes (M, AMM,,..,MM
F : y E ::: Vj E [
arg argmax s(2, P,y') M,)
(1)
ZEZ
for some projection matrix P,, to be defined later.
3

---

Published in Transactions on Machine Learning Research (06/2023)
(iv) It approximates the conditional mull distribution of the test statistics by sampling biological outcomes y' under the null, conditionally to the selection event. This sampling is performed using hit-and-run strategy (according to Algorithm 2), by building a discrete time Markov chain on E whose distribution converges to the uniform one.
(v) SEISM finally computes the p-values for the mull hypotheses defined (iit), associated with the selected motifs in ii, using the empirical distribution of the test statistics, and returns the motifs with their association p-values. Given these p-valucs, one can adjust the number of selected motifs discarding the ones with non-significant p-values. This multiple-testing issue has not been investigated in this paper, but the practitioner can use for instance a Bonferroni bound to select the number of motifs.
The data-split version SEISM applies the same (i)-(ii) steps on a fraction of the data and simply compares the scores of the selected motifs on the remaining data to the distribution of scores for the same motif with data sampled under the mull distribution-ass opposed to the selective mull generated by (iii)-(v). Sampling is much faster under the null than under the selective null, because it does not involve rejection step. Both samplings will need our results in Section 4.6 to avoid depending on particular value of the mean and variance parameters.
Algorithm SEISM algorithm (general formulation)
# Description: SEISM selects a sel of sequence molifs (zi,..., ;g)) based on an associalion score s(,,), and evaluate their p-volues based on a partition = M,.
Response y E V C RR, sequence samples X, feature function z E Z + xp E R", association score 8: 2 x y + IR, number of selected molifs q 2 1, meshes Z = Mi, sampling algoritlm HR.
Inputs:
==1
Result: ((p1, .......
... (pg. 2n)), sequence of p-values and sequence motifs.
# Selection step: Selection of the sequence motifs (Z1,..., 2g) and the sequence of meshes (Miy--.., My.)
for j = .,...,9 do
2j arg max szz P,y)
11 Pk orthogonal projection onto Span {
262
tej // the mesh Mi, is selected
in f i s.t. 2j E Mi ;
end
# Inference step: SEISM provides p-value ps the statistical inftuence of the selected sequence motifs z conditional on the selection event (1) of observations y' that would have led to same selection of the sequence of meshes AM,,,..,,M,,
y (1) (N) + HR(y, (if....,,,.)).
// Sampling outcomes under the selected mull
for j =
do null: .....
F (1)
....y y ((N)) +- empirical cumulative distribution function of s(r H₁y') under the selected 11 rM;, is motif representing Mi, and Iī the orthogonal projection onto
Span { [x=
Pj - P3 (s( 1 .... (N))
// output the jth p-value
end
3 One hidden layer CNNs select sequence motifs maximizing an association score
One-layer CNNs have been at the core of the rising popularity of deep learning over the past decade, by enabling major improvements in computer vision tasks (Krizhevsky et al., 2012). Although they are formally a specialized fully connected foodforward networks with additional constraints on the weights, CNNs are equivalent to, and more often thought of as, of convolutions of the vectorial input with some smaller vectors referred to filters. When applying the network, dot products taken between each of them and
4

---

Published in Transactions on Machine Learning Research (06/2023)
successive windows of the vectorial input followed by some non-linear operation, producing an activation profile for each filter. In one-layer networks, these activations are pooled across the windows into a single scalar for each filter and these scalars are combined typically through a linear or regular fully connected network-to provide a prediction for the input. Because convolution filters are homogeneous to the input, they easily lend themselves to interpretation: as small image patches for image inpuls, and as sequence motifs for appropriately encoded biological sequence imputs. Accordingly, activation profiles reflect how much each piece of the input is similar the filter in the sense of the dot product and applying a one- layer CNN amounss to applying predictive function to modified representation of the original data by these similarity profiles. Because convolution filters are jointly optimized with the parameterization of the predictive function, CNNs are often described as a strategy to jointly learn data representation and a function acting on this represenlation, both being oplimized for a prediction objective. In computer vision, the optimized filters of the first layer typically learn to detect edges with different orientations. In biological soquenees, they learn short sequences whose presence anywhere in the input is predictive of the output phenotype used for training.
A
00.11
0.27
0.01
0.99
0.67
0.77
0.00
0.00
I ITAAG GGG
C
0.00
0.48
0.00
0.00
0.25
0.01
0.60
0.00
G
0.01
0.17
0.13
0.00
0.05
0.21
0.40
0.00
T
(0.88
0.08
0.86
0.01
0.03
0.01
0.00
1.00
Figure 2: A motif represented by its position weight matrix and corresponding sequence logo. The total height of the letters indicates the information content of the position (in bits), closcly related to the Shannon entropy.
Unlike imput sequences that are formed by discrete succession of letters in some alphabet, trained filters are continuous and therefore account for possible variation in the predictive short sequence, e.g., T mostly followed by C but sometimes an A or G and so on (Figure 2). These probabilistic objects have also been used for a long time in the bioinformatics literature and referred to as position weight matrices (PWMs). Inferring PWMs either according to their frequency in a set of sequences (Bailey et al., 2006) or their discriminating power between two sets (Bailey, 2021) has been a major theme over the past thirty years. Here we formalize the training a one-layer CNN as equivalent to the selection of set of sequence motifs that are optimal for some association score. This formalization will be instrumental in the definition of our hypothesis lesting procedure in Section 4.
Notations Let X represent a data set of n one-hot encoded sequence samples {æ1, 22.......... In), in set X of biological sequences assumed to be over an alphabet A A-for DNA sequences, A= = (A,C, encoding maps each letter in AA to a vector in {0, 1014,, with all-zero entries except for a single 1 at the coordinate corresponding to the order of the letter in A for DNA sequences, A is encoded as (1,0, 0, 0). Every 24 is therefore encoded as a matrix in {0, 1)1A|x [ although in practice, encoded sequences are often padded with dummy columns to have the same lengths. We denote Y/ y the measurement of a biological property associated with sequence æi, and y E ym the corresponding vector of outcomes. We consider one- layer CNNs with Gaussian non- lineariy with scale w, max global pooling and a linear prediction function. These CNNs parameterize a function f X + y by filters of length k namely Z = {II...., 2n} E z", where 2 is a subsct of A|x given by the simplex in this paper:
One-hot
IA [ Σ Zij =1 =1]
f=x2
2 2
Z
(2)
and g weights B E IR9.
5

---

Published in Transactions on Machine Learning Research (06/2023)
More precisely, define f(x))
Cn where Cn In 1 is the centering operator, In the identity matrix, 1 the all-one vector in R and
B), with up x E Rexe defined 2
||= ull?
f 2.X
max uEieile
exp
2
where [as]le denotes the set of l consecutive entries of the vector as (and of its reverse-complemmnt coun- terpart), and is bandwidth hyperparameter whose impact and tuning is studied in Appendix A. This model differs with typical CNN in two ways. First, it uses Gaussian activation function instead of an exponential one; second the use of the centering operator that sets the average of the activation to zero. These adjustments were made to improve the SEISM algorithm's selection performances.
3.1 From empirical risk minimization to association scores
The function f is learned in a classical penalized empirical risk minimization framework, using the data {X,y):
min (Z. (Z)8)
n 1 22 + 111
(4)
for some À > 0. Equation (4) formalizes the idea that learning a one-layer CNN on one-hot encoded sequences amounts to learning a data-representation 2 of the sequences parameterized by a set Z of filters corresponding to PWMs and a linear function with weights 3 acting on this representation. Noting that there exists a unique explicit optimal / for Eq. (4), follows immediately that:
arg min
{ min (n "ly Z. B1 + {([l11a)
arg max z
ridge (Z,
(5)
where gridge defines a particular quadratic association score between an outcome y and set of filters Z:
ridge 2 (2.80Z y) = y² z [1,2.X)TT + AmI,] ((?2.X"Tp
It formalizes the training of a CNN as the selection of a set of filters whose association with y in the sense ridge is maximal. Of note, one has of 6x
ridge lim An x x (Z.y) = (x22X)) = e HSIC ( (Z. y) X
so for large values of the regularization hyperparameter selecting filters by learning a CNN is equivalent to selecting filters with the classical HSIC score (Song et al., 2012), because already includes a centering operator. In addition to connecting ridee with HSIC, we observed that the centering in the definition of yp 2 led the selection of better sequence motifs our experiments. Observe that the centering matrix is an orthogonal projection matrix onto E == Range(C,, the orthogonal of the vector line generated by the vector 1, and it holds
y z.X 3 = ||C.y s0 z 11. + |1 C.yle.
The solution of (4) is unchanged if y is replaced by Cay, and generality loss. Furthermore, this shows that can work with skewed data in classification context, since imbalanced classes will have no effect on the result.
assume that y & without any
3.2 Greedy optimization
It is common to solve (4) by stochastic gradient descent (SGD) jointly over the filters. More generally, this approach for training a neural network with single, large hidden layer is known to find a global optimizer at the large limit under some assumptions (Soltanolkolabi el al., 2019). Our objective here is slightly different: we do not necessarily aim at approximating a continuuus measure with a large number of particules, but we aim at selecting small number of particules lending themsclves to a biological interpretation. Furthermore, the number of relevant motifs on a given dataset is generally unknown. In this context, known that jointly optimizing the convolution filters leads to irrelevant PWMs, with some actual iot split across several filters

---

Published in Transactions on Machine Learning Research (06/2023)
and other duplicated (Koo & Eddy, 2019). A possible strategy is to forego filter-level interpretation, train overparameterized network- with a much larger q than the expected number of motifs-aai use attribution methods to extract relevant motifs or other interpretable features from the trained network (Shrikumar et al., 2018). Here we adopt a different strategy using a forward stepwise procedure, where we iteratively optimize each of the convolution filters over the residual error left by the previous ones.
More precisely at each of the q steps, we select 2j such that:
2j arg max mww.stiitaeee.... stidge (z, Pzy),
(8)
#EZ
where Pj is the projection operator onto the orthogonal of the subspace Span {1, 11, see line 2 of i<j
Algorithm This is how 2; is optimized over the residuals of the previous filters. The vector 1 enforces that we project on a subspace of E, in particular P1 C.. Without this projection, iterating (8) would return the same Of note, joint optimization proced of the q filters don't face this issue, and forward selection procedures over finite sets of features work around the problem by iteratively removing the selected elements from the set over which selection if performed (Slim et al., 2019). This sequential strategy combined with the testing procedure introduces in Section 4 provides a data-driven mean to choose the number q of relevant motifs.
In practice, we solve (8) with a slandard gradient descent algorithm, initialized at the k-mer with the best association score. The k-mer list is obtained using the DSK software (Rizk al., 2013). The length k first varies according to uscr-defined range, and the optimal value is chosen by SEISM, as described in Appendix A. We work on a less constrained set than Z (2) and don't 't enforce the positivity constraint during optimization. We project the optimized motifs onto the full Z at the end of the process. Our procedure also requires to choose a motif length k. We proceed adaptively by choosing the length leading to the highest score, within a user-specified range.
With the one-layer CNNs training formally cast as the successive selection of q sequence motifs optimizing an association score, we now turn to the problem of testing the significance of these associations. Of note, what follows is only based on the definition of an association score and could be applied to perform inference on heher features coming from the training step of any algorithm, as long as one can define an association score between the feature and the outcome.
4 Post-selection testing of the association between the outcome and trained convolution filters
We now turn to the problem of testing the association between the selected motils z and the trait y. Tn order to do so, we need to solve three interrelated problems. First, the motifs were specifically selected for their association with the trait, which leads to the well known post-selection inference problem. Any inference procedure that disregards that the hypothesis was constructed using the same data used for testing is likely invalid and produces deflated p-values. Second, we deal with a continuous selection event, because (8) is performed over a continuous set Z. By contrast, existing solutions for post-selection inference address selections over finite sets. Third, the null hypothesis commonly used for similar post-selection inference problems is composite, i.e., it corresponds to several values of the parameters. Existing methods work around this issue by fixing theses paramcters to arbitrary values, thereby limiting the scope under which they are calibrated. Here we present our solutions lo these three problems.
Consider the Gaussian model:
y = μ + dC
(9)
where μ the target deterministic signal, and c ~ NO, C.,) the standard Gaussian distribution on E.
4.1 Selective null hypothesis
We follow Yamada al. (2018) and test the association of motif z through the following mull hypothesis:
Ho "s(2, 11) 0"
(10)

---

Published in Transactions on Machine Learning Research (06/2023)
for some association score For z chosen independently of the data, Hg could be tested by sampling y' under the corresponding distribution, and using the quantile of the s(2,y/) scores corresponding to s(z, y) a p-value- i.e., the probability when sampling under Ho to observe score as extreme as s(z, y). In our case, however, the motifs & in the trained convolution filters were specifically selected for their strong association with y, and this procedure would not produce calibrated p-values. This problem known as post-selection inference, and has been discussed and addressod a growing literature over the past few years. Data-split strategies lead to valid inference but necessarily result in a reduction of the sample size, which is unsatisfying when the original sample size is limited. Alternatively, selective inference frameworks were developed in the recent years to address these issues. We refer to (Hastie et al., 2015, Chapter 6) and references therein for a general presentation. Taylor et al. (2014) and Lee et al. (2016) address scenarios where the selection event, e. the set of data outputs that would result in the selection of the same set of features, is polyhedral determined by the finite intersection of linear constraints. Reid & Tibshirani (2013), and later Reid et al. (2015) extend this selection to clusters or groups of features, still in the linear framework. Yamada et al. (2018) extended post-selection inference to the non-linear framework, by proposing kernel-based approach, where the selection is performed through the HSIC criterion. Slim et al. (2019) generalize this work, by allowing the selection to be carried out with wider range of tools, making use of quadratic association scores.
To our knowledge, the selective inference literature only addresses the problem of selecting features from a discrete collection and does not provide a solution for selections from a continuous set like our 2. Hence, testing (10) directly is not feasible and we resort to the quantization of the motif space to address this problem.
In addition to that, we push the analysis of the statistical model further, in order to be able to apply it with weaker assumptions on the data distribution.
4.2 Selective inference over continuous set of features
Formally, our selection event Evont. is the set of outcomes y' that would have led to the selection the same set of motifs Z = (zI...., 2g} than the one selected using y from the real dataset, when applying the same selection procodure:
Econt. : E : Vj .l......hh.... arg max s(2, P,y') = =z;}
(1I)
2E2
where P, is the orthogonal projection onto Span {1, x=
Kij
A simple rejection approach to sample from the null (10) conditioned to Feont. would be to sample y in E under (9, 10) and only retain those in Econt. Unfortunately, Econt. belongs to a strictly lower-dinensional vector space of R" and is therefore a mull set for the Lebesgue mcasure on R For s MSIC and sridge, and noting that a maximum also a critical point, we indeed obtain:
y' E Fcout. = Vj E (1,..,,4)) Pzy' E Span { [V=
For Q = and assuming that the different directions of the gradient. are independent, this spans is a vector subspace with dimension n 4 k. We empirically observed that sampling from this subspace produced non-zero proportion of y' in Econt.. Nonetheless, choosing a sampling distribution that leads to the correct conditional distribution is not straightforward-and and may not even be possible-ss discussed in Supplementary Material B. Moreover, relying on conditional probability with respect to a null set is not well defined and may lead to the Bole--Kolmogorov paradox (Bungert & Wacker, 2022), which further complicates its use.
We choose to circumvent this issue using a partition of the space of motifs spaces, over which our selec- tion (8) operates, into very large but finite set of meshes: Z MM, M,. As depicted in Figure 3, we consider
regular partition of each coordinates into m bins:
8

---

Published in Transactions on Machine Learning Research (06/2023)
A
0.75
0.25
e
A
0.5
0.5
C
0.75
IC
0.75
0.5
0.25
G
Figure 3: Discretization of the 3-letters alphabet simplex {A,C, G), with a binning parameter for the meshes m 4.
Based on this partition into meshes, we define a quantized selection event E as follows. First, given an outcome y we define the sequence of the q selected meshes (M,,..., M.,) as
Vj E f1. arg max s(z, P,y) My;
#E2
Second, the selection event is given by:
E(i in) = y e:: Vj (I....,.
s(, E Mr.
(12)
ZE2
the set of outcomes y' that would have led to the selection of motifs within the same meshes as the selected ones (M,,...,MM,)).
We now show how quantization (12) of the selection problem make enables the definition of a valid inference procedure. We start with the simplest case where we select a single motif (q 1).
4.3
Test with only one motif 1, / and fixed
In this section, considering the motif was chosen by the SEISM selection procedure, selection event (12) boils down to:
B(i { y/ y
arg max s(2,8') E Mi,
(13)
#EZ
We use this simplified case introduce our mull hypotheses and test statistics attached to this selection event, and consider options:
A first option consists in representing the mesh M,, by its center e₁ Then the corresponding mull hypothesis is the following:
H², : "s(c1, 4) = 0", It can be tested using statistic V = s(c1. y).
(14)
A second possibility is to represent Mi, by the motil with the highest association score within. In this case, the null hypothesis becomes:
H
"V Mi s(2, 0",
(15)
We test it using statistic Vi"
max s(2, y)
zē M;
9

---

Published in Transactions on Machine Learning Research (06/2023)
Algorithm 2 Hypersphere Directions hit-and-run sampler
/+ Description:
The Hypersphere Directions hit-and-run sampler creates a discrete-time Markov chain on an open and bounded region and is used to approximate a uniform distribution on the selection event F.
Inputs: Response W E C R", B and R the numbers of burn-in iterations and replicates.
Result: y y((+1) y (B+R) E EECR" R" the replicates sampled under the conditional mull distribution
10 y (0) + L(y); /* L is the cumulative distribution function of O²
11 for t= ..... B+ + do
Sample uniformly ( from (0 E R",
||0| 1}:
12
(1-1)
y (t-1)
(1) - max
13
max
max
~0
0,
0 <0
02
(t-1)
1 y (
14
600)
max
min
min
/* Sampling (() from (a(t)661) ensures that
0 <D
00/00 >0
Or
g((t-1) + ((00 1
while y'th (00 & E do
15
/* This loop parallelized on several cores until one of them discovers replicates in the selection event.
Sample uniformly X from (09))(0)) [;0);
16
j'(t) ÿj (1(+11 + (00)
17
y (t) + - L-
18
19
end
20 end
In both cases, we reject the mull hypothesis the test statistics are greater than a threshold, determined
by their cumulative distributions under the nulls (14), (15) conditionally to E(i)) : F Fl. and F. FT. In practice, there is no closed form for these conditional cumulative distributions, and we rely on empirical version that we build using a hit-and-run sampler algorithm, as described in Section 4.4.
Hypotheses (14) and (15) lead to very similar results when the meshes are small enough, which is easily the case in practice. (14) gives us insights on one specifie motif of the mesh the center, but (15) tells us about illustrate the difference, let us
whether there exists a motif within Mi associated with the phenotype. consider a meshing with only one bin per coordinate, that is the meshing with only one mesh, containing all the motifs:
Testing the center-based null hypothesis (14) boils down to testing the association of μs with the motif e with the same probabilities for each letter of A at every position, and produces a p-value of 1, regardless of the data, since for any k-mer u, ||c u||2= k (0.75² + x 0.252), which leads to XX X, x = according to the centering step, and to zero score for any y E
By contrast, one can obtain strictly less than value for (15), because different z E can lead lo different scores, which means that there may exist motif in 2 associated with y- -ut does not inform us which motif it is.
4.4 Sampling from the conditional null distribution with the Hit- fit--and-Ru Run algorithm
Even after reducing our selection to a finite set (Section 4.2), a rejection sampling strategy that would draw y' from either (9, 16) or (9, 17) and only retain those leading to the selection of the same mesh as y is not tractable as the rejection rate is empirically too low. Following Slim et al. (2019), we resort a Hypersphere Direction strategy (Algorithm 2).
10

---

Published in Transactions on Machine Learning Research (06/2023)
The hit-and-run algorithm produces uniform samples from open and bounded acceptance region corresponding, in our case, to the selection event. It starts from any point in the acceptance region, draws a random direction from this point and samples along this direction until it finds one elements that also falls in the acceptance region. then follows the same procedure from this new starting point. The hit-and-run sampler therefore also relies on rejection but it does so along single dimension rather than from R It explores the selection event step by step, starting from a point that belongs to this event, which guarantoes higher acceptance rate. To speed up the procedure, we parallelize the rejection step across several cores. Because each point sampled by the hit-and-run procedure depends the previous one, it is impossible to parallelize the whole sampling process. By contrast, the rejection step used for computing single replicate, once & sampling direction has been fixed, can be parallelized. We draw several distances to the initial point independently, oplimizing new independent points, until one of them belongs to the selection event. This parallelization provides a significant time saving, as discussed in Section 5.3. Algorithm 2 produces uniform samples from an open and bounded acceptance region. The boundedness assumption does not hold in our case the arg max over of the score only depends on the direction of y and not on its norm. The openness requirement is ensured by the definition of the meshes. Following Slim et al. (2019) again, the reparameterization and ou denotes the cumulative distribution function of (M,, o²C,). Sampling uniform ÿ from the open bounded space j0, 1(th indirectly provides normal samples from ((,,o
use
Lyy), where L R" }00, 1["h is defined as L(y); L (3%) for i = I,...,
Combining this sampling strategy with the quantization of the selection event introduced in Section 4.2 and the selective mull hypotheses attached to this event introduced in Section 4.3 provides a selective inference procedure for selected motif ø1 (g = 1) and null defined by given pair (n, of parameters. Our next two steps to handle the selection of multiple motifs, and the general case where several μ describe the same null hypothesis and a is not specified.
4.5 Dealing with the selection of several motifs (g > 1)
We now consider that we selected > motils with the SEISM procedure, leading to the general (12) selection event E(i Generalizing our single-motif strategy of Section 4.3, we propose two options for defining null hypotheses (and test statistics) related to this selection event:
The first one relies the centers of the selected meshes:
Ho, "s(c,, [T 0"
(16)
where IT is the orthogonal projector onto Span (xt In other words, it expresses that the center of the mesh M., is associated with μ after removing its component carried by the span of the centers of the meshes corresponding the the q - other motifs.
And the second one takes advantages of the best motifs in each mesh:
Ho "V(zi))ed (Mi))ess Vz Ei,, s(2, ((zt))+s) p)) -0",
(17)
with II" ((=t) being the projection onto Span { x
Generalizing what we introduced for g = 1 (Section 4.3), we test those hypotheses using V = s(,, IT,,, and V"" = max s(2, I) To that end, we rely on their cumulative distributions under the nulls (16), (17) conditionally to E(ii,..., ig) respectively Fy (i....( and IF" (.....((i empirically approximated with Algorithm 2.
Following the work of Loftus & Taylor (2015) in the finite case, both versions of our mull hypothesis are joint across the q motifs: each of them considers the association between the j-th selected motif and after projecting onto the span of all others, not just the ones that were selected before using IT' and II". This is to contrasted to our sequential selection process, which adjusts at each steps for the previously selected motifs using P.
In order to give more insights these null hypotheses, we derive the following proposition:
11

---

Published in Transactions on Machine Learning Research (06/2023)
Proposition 4.1 (Description of the selective mills), Let Z = {z1...... zg) be q sequence motifs. Let s(:,)) be score such that 'nullily implies orthogonalily" (for instance slii or sridge):
(A1) Nullity implies orthogonality: If {s(z, y) = 0f then {(0*,x yy = =0, for every (y. z) E 2, and for some function z - + E E.
Let μ E and decompose μ as
ayp= + 1
(18)
with μt E & orthogonal Span((?2.x)
It holds that "s( [1)pe) = 0" is equivalent to "a== 0 for some decomposition (18).
If Rank(y = q then the decomposition (18) is unique, and the greedy selection procedure described in Section 3 enforces this situation. We interpret this as follows: we look a motif žạ and would like test its signíficance; in view of property (A1), can climinate the effects that captured by the other motifs by using the orthogonal projection onto the orthogonal of Span(y>") ), given by II (using II, II'
or I1, = =II" ((zt))+s)) and consider y to test the association "s(z). [I) = 0"; equivalent to testing "aj- = 0" by the above proposition.
4.6 Sampling under selective multiple hypotheses with known a
The sampling strategy described in Section 4.4 builds a conditional mull distribution- selective inference procedure-for a given μ and a. In practice, a is not known, and several values of μ can describe the selective null hypotheses (16) or (17) for a given motif selection. Of note, this issue is paapiingg-asdddd not specific to our selective inference procedure. It will arise in any sampling- based post-selection inference strategy including data-split: even if the latter samples from non-sclective mull hypothesis, it still needs specific values for Ju and a. For μt, the issue is that any fixed value can not represent the whole set of possible values, which would modify the null hypothesis actually tested. For the variance parameter, it may be fixed by the user, but this may lead to a non-valid procedure if the chosen value is different from the real one.
offering a
We leave aside the choice of a for now, and describe how we can sample from any null distribution (16) or (17) using p = for given a. Our results holds for scores verifying the following assumption- this includes both HSIC and gridge;
(Az) Nullity implies trauslation-invariant:: If sf v) = then Vy E E s(z,v') s y + y'), for every (y. E E x 2;
Under this assumption, the following proposition ensures that using the quantile of the empirical distribution of scores sampled under μ = 0 leads to calibrated test procedure:
Proposition 4.2. Let be an association score such that (A2) holds. Let Vi = s(c,, IT,y) and V = ma 8(z, y) formed from y sampled from (9) with any mean la such that s(z', = 0, any known variance > 0, and such that z' arg magxaxzzz s(2,8) y). The conditional rall distributions F' (.(.....)))) E= (.(11....0 with mean 0 and variance o verify:
and
F (V) ~ Unif(0, 1) and Fy (V,)) ~ Uni 1)
Proof. Assumption (A2) under the Gaussian model (9) implies the following property:
V(2, A, E Z x A x & such that y = + oe,
(19)
"s(z, A/i) = 0" - "s(z, As( = s(2, a AE)",
12

---

Published in Transactions on Machine Learning Research (06/2023)
which implics that, for composite mull hypothesis of the form Hy "s(2, Ap) = 0", the distribution of s(z Ay) does not depend on the mean μ that satisfies Ho. Hence, even if the hypothesis Ho corresponds to a set of probability distributions of that may depend on the distribution of the statistic Ay) does not depend on under this hypothesis. We can then conclude that if a is known, as it is assumed to be the case in this section, then a lest statistic of the form V = =(z, Hy) has the same distribution as s(2,olle).
4.7 Sampling under selective multiple hypotheses with unknown o
In practice, a often unknown. To address this issue, we rely on the normalized versions of the test statistics V and V" introduced in Section 4.3, defined by
2, 12, ((E))))))))
s(c,, |y/l2
Tý
and Ty
max ZEM;
(20)
|/
where z5 arg 2 szz Pay). We will denote G Ci and G 5,65)))))))) their cumulative distribution functions under the null, conditionally to E(ii....,, (g).
We will also make use of third assumption, here again fulfilled by HSIC and stidge
(As) Two-homogeneity: It holds that s(z, ty) s2 y) for all (y, z) & x Z and all t
Of note, normalizing the association score with respect to the labels does not affect the selection:
s(z. y)
Vy E V, arg magmax s(2, y) arg max =EZ
(21)
y
ZE2
If 4 14=00 O, the distribution of the normalized statistics does not depend on 0, and the empirical cumulative distribution functions of normalized scores obtained by sampling under μa 0 and any a still provide valid inference procedure:
Proposition 4.3. Let s be an association score such that (Az) and (As) hold. Let T2 = s(c,, /T /lll || and T = max s(z, y)/llll formed from y sampled from (9) with mean μ 0, and any variance > 0. Then for all of > 0, their conditional null distributions G (. and G" (.66).... with mean 0 and variance o' verify:
G' ()( ~ Unif(0, 1) and G (" ~ (0,
Proof. Let us consider two different normal models as defined in (9) under the global null hypothesis "p 0" and given by
()) = (1) and y (2) = (2) _2))
Then
s(c,, 1y |1
s(c,, s(c,IIy,' y (2)
and s(z, 11" (z)))ts) ((
s(2, II" ((z))es) y121)
|1
21
The proof directly follows assumption (Aa) applied with t 112. Proposition 4.3 is complementary to Proposition 4.2 and provides a selective inference procedure when o is unknown, under the special mull hypothesis fa = 0.
Our final result investigates the lesting procedures for the general null hypotheses (16) and (17) not re- stricted to μ = 0- with an unknown a. Recall that the decision rule is to reject the null hypothesis the observed value of the statistic is greater than a given threshold We show that choosing t to be a quantile for the global mull hypothesis (μ = 0) leads to calibrated (for the type I error) non-sclective procedure, see (22).
13

---

Published in Transactions on Machine Learning Research (06/2023)
Proposition 4.4 (Global mull achieves lowest observed significance). Let Z (21,...,,2 ................................ ,,) be q sequence molifs. Let -) ) be score such that (A1) and (A2) hold. Let pa E & be such that
Ho "s 1 1e) = =0
Then
8 Z 11 + oE) p
8 Z E)
v4 > 0, sup P
2t
=P
(22)
aell
He
We provide proof in Appendix C. This proof makes an ad-hoc use of Anderson's 's theorem on symmetric convex cone (whereas it is usually devoted to symmetric convex bodies).
Proposition 4.4 shows that data-split produces a calibrated procedure for testing the general mull hypothe- ses (16) and (17) when sampling under the global null (n 0) the test statistics (20). We could not prove equivalent statement for conditional null hypotheses, and Proposition 4.4 therefore does not guarantee the validity of selective inference procedure sampling under the global null (μ = 0). Yet, we used it heuristic justification of SEISM and we observed that it leads to empirically calibrated procedures, see Section 5.2.
In view of Proposition 4.4 and its proof, one can see that the alternatives i such that ||Pa is large have small power. As the selection procedure described in Section 3 achieves good results (Section 5.1), the chosen motifs Z should capture the principal components of μr, and therefore such that ||///1//lll) should be small.
5 Results
5.1 SEISM performs as well as state-of-the-art de novo motif discovery methods
In order to compare the accuracy of our selection step with existing motif discovery algorithms, we use the 40 ENCODE Transcription Factors ChIP-seq datasets from K562 cells (ENCODE Project Consortium, 2004), each of which contains a known TF motif, denoted m², derived using completely independent assays (Jolma et al., 2013). STREME (Bailey, 2021) and MEME (Bailey et al., 2006) are state-of-art bioinformatics methods for de-novo motifs discovery tasks. STREME identifies motifs that maximize a Fisher score of association between the presence of the motif and the binary class of sequences. By looking for maximum likelihood estimates of the parameters of mixture model made up of background distribution and a model for generating k-mers al some positions that may have produced particular dataset using an expectation maximisation technique, MEME finds enriched motifs in this dataset. Finally CKN-seq (Chen et al., 2017) is a one-layer CNN tailored to small scale datascts. We set up STREME, MEME and SEISM to select 5 sequence motifs. SEISM is run with a regularization parameter X = 0.01. CKN-seq jointly optimizes its filters, which notoriously leads to poor performances when few filters are used. We train it consequently over 128 filters. We measure these accuracy of all methods by comparing the motifs they discover with the known motif corresponding to the transcription factor m.. We rely on the Tomtom method (Gupta et al., 2007), which quantifies the probability that the euclidean distance between random motif and m is lower than the distance between the discovered motif and m More precisely for each method we use the lowest. Tomtom p-value between the known TF motif m and any of those discovered by the method. The Tomtom score is then defined as - loggo of the Tomtom p-value. We define the accuracy of the method as the proportion of experiments where the Tomtom score between best match and the true TF motif was higher than some threshold.
Figure 4 (left panel) demonstrates that SEISM is just as good as, if not superior to, state-of-the-art bioinfor- matics algorithms at detecting sequence motifs when thresholding Tomtom p-values at 0.01. The one-layer CNN with jointly optimized filters performs poorly in this experiments, emphasising the importance of greedy optimization for selecting the right motif.
Figure 4 (right panel) shows that SEISM performs slightly worse than STREME and MEME for high thresholds on the Tomtom scores. This suggests that the matrix that SEISM identifies is close enough the PWM matrix of the true motif, but farther away than the matrices identified by STREME or MEME.
14

---

Published in Transactions on Machine Learning Research (06/2023)
100
SEISM
STREME
CKN
(N) TOMTOM
Succes Patte 1961
60
20
20-
SEISM
STREME
accuracy score (llog0 TOMTOM puslless)
Figure 4: Left: Proportion of datasets where the true motif was detected by the designated algorithm. A true motif said to detected its highest Tomtom score with the discovered motifs greater than Right: Accuracy motif discovery algorithms on ENCODE TF ChIP-soq datasots. The curves display the proportion of ChIP-seq datasets were the best motif identified by the designated algorithm has Tomtom score greater than
This discrepancy reflects a different usage of z to parameterize a distribution k-mers. In practice, we observe that on a given dataset, the p-values of the best motils discovered by SEISM and STREME are not separated by more than 2 orders of magnitude, which leads minor differences in the motifs, as illustrated in Table 1.
Reference motif (ATF4 DBD)
SEISM motif p-value= 10-6
STREME motif p-value= 3 x 10-8
ITG
TCL 8 T A
CACA
GATCA
CA
Table 1: Comparison between two discovered sequence motifs by SEISM STREME, and the true motif (ATF4 DBD)
Both SEISM and MEME/STREME exploit a distribution of k-mers at the transcription factor binding site. MEME and STREME maximize the likelihood of a categorical model, whereby the matrix z directly defines the probability observe each letter each of the sites:
Vu, z E Leas 2))
II u
(23)
SEISM on the other hand is based on a Gaussian model. Through representation (3), z is meant to maximize the Gaussian likelihood of a set of k-mers, i.e.
Y(u, z ZE 2) Luaus(u; C II
(24)
where is a constant such that the sum of probabilities over B4xk equals 1. If we consider binary y match the setting of MEME/STREMF, this set made of one k-mer for each positive sequence. Importantly,
15

---

Published in Transactions on Machine Learning Research (06/2023)
the true TF motifs from (Jolma et al., 2013) that we use to assess selection accuracies are also defined through the maximum likelihood in categorical model, which can explain why the obtained with MEME/STTRME are closer to the true PWM than the one obtained with SEISM.
We now illustrate on a simple example how the same distribution of k-mers is parameterized by different matrices under the two models. To build an easy example, we focus on k-mers of length 1, with
P(A) 0.3, P(C) = 0.4, P(G) 0.1, P(T) 0.2
(25)
The matrix 2 0.4, 0.2)* used with the categorical model trivially constructs such a distribution. But using the same matrix in a Gaussian model with parameter fixed as described in Appendix A leads to slightly different distribution:
P(A) 0.28, P(C) 0.43, P(G) 0.11, P(T) 0.18
(26)
A distribution closer to Equation (25) can be constructed with a Gaussian model parameterized by 22 = (0.315, 0.38, 0.08, 0.225)".
To clarify the relationships between those two motifs, we will rewrite (23) considering u is one hot encoded. That is, for each position i, it has only one for letter j(i) and 0's elsewhere:
Yuu, z ZE 2), Leat 2) = II 2i,(()
(27)
Assuming that the columns of are normalized and 1, we can modify (24):
Vuu, 2 ZE 2). Lgaus (u: 2) = C II e 1-
C2 II xF= C2 II FU00
(28)
With the Gaussian model and a few assumptions, the motifs can the be seen as defining the log probability to observe each letter at cach of the & sites. This gives us new interpretation for the filters learned by CNNs and suggests that in this framework it might be interesting to constrain e* ratherr than to be in 2.
We used a Gaussian activation function since it closer to typical CNNs approaches. Our framework is generic enough to allow other activation functions basod on the categorical model, or more rcalistic variants (Ruan & Stormo, 2017).
5.2 Statistical validity and performances
In order to assess the statistical validity and of the SEISM procedure with the different strategies, we simulate datasets under the null hypothesis. To that end, we draw one sequence motif with length k == 8 for each simulated datasct using uniform distribution 10 bits. Then, we draw a set of 30 biological sequences X as follows: all sites are generated according to uniform distribution over A, C, T. G for all sequences, and for half of the sequences one k-mer is drawn according to the categorical model parameterized by ž. The phenotypes y are drawn from ((O,)"G.) to generate data under the null hypothesis for calibration experiments, and from N(yxE. ,o'Cn) to generate data under the alternative for experiments on statistical power, with a = 0.1 in both cases. We then run the SEISM procodure to seloct and test two sequence motifs. For both the data-split strategy and the hypersphere direction sampling one, the distribution from which the replicates are drawn uses the empirical variance from y as variance parameter. Although any choice for this parameter leads to valid procedure, as described in Section 4.6, we make this choice for numerical stability considerations. For the data-split strategy, we sample 1000 replicates under the null hypothesis to compute the p-value. For SEISM, we sample 50, 000 replicates under the conditional null hypothesis using the hypersphere direction sampler, after 10, 000 burn-in iterations.
Z restricted to motifs with information level fixed
Figure 5 (top) shows the Q-Q plot of the distribution of quantiles of the uniform distribution against the p-values obtained across 1000 datasets under the null hypothesis for the data-split strategy and 100 datasets
16

---

Published in Transactions on Machine Learning Research (06/2023)
for the hypersphere direction sampling one. All the data points well-aligned with the diagonal, which confirms the correct calibration of both the data-split and hypersphere direction sampling strategies, either considering the best motif or the center of the mesh and regardless of the size parameter.
Figure 5 (bottom) shows the same Q-Q plot on data generated under the alternative hypothesis. From this figure, we observe that on small datasets, the post-selection strategy more powerful than the data- split one, regardless of the size of the mesh considered or the choice concerning the definition of the null hypothesis. The variance observed on the curves associated with the selective inference procedure is due to the presence of weak residual signal after the first motif result of imperfect selection step. Testing it with the best motif in the mesh captures this signal, resulting in curves under the diagonal. By contrast, focusing on the center of the meshes leads to testing motifs that do not capture this signal, placing us in the conservative situation, described at the end of Section 4.7. The residual signal is not well explained by the mesh's centers, and thus its component on the orthogonal of the span of the activation vector of the second motif important. The larger the mesh, the farther its center is to the selected motif and thus the less signal it captures, which explains the differences between the two
Motif 1
Motif 2
1.0
1.0
me
0.8
mw0.3
motil
Calibration
0.4
0.4
0.2
p-values
1.0
1.0
0.8
Power
0.4
0.4
Empirical quantiles
Figure 5: Q-Q plots obtained by applying data-split and different hit-and-run sampling strategies to select two motifs and test their association with an outcome. Top: data simulated under the null hypothesis. The proximity between the quantiles of the obtained p-values and those the uniform distribution confirms that all SEISM strategies presented in this article are correctly calibrated. Bottom: data simulated under an alternative hypothesis, where the outcome depends on the activation p of a single motif in the sequence. The distribut ions of the p-values computed with the post-selection inference (PST) strategies have a larger deviation to the uniform distribution than the distributions of the p-values computed with the data-split strategy (purple).
17

---

Published in Transactions on Machine Learning Research (06/2023)
5.3 Computation costs
The section serves as an overview of how various user-specified parameters impact the computation time roquired by the post-seloction inference procedure.
As discussed in 4.6, the hit-and-run algorithm is actually a rejection sampler. Its overall computation cost depends mainly on two characteristics: the cost of the selection step, that is the cost of selecting g motifs for given phenotype y, and the acceptance rate. Although some paramcters affect the selection cost, the acceptance rate is primarily responsible for determining if a user-specified combination of parameters results in a tractable configuration for the post-selection method in a reasonable amount of time. This rate is high compared with a naive rejection sampler over &, as the hit-and-run strategy reduces the dimension over which the rejection step performed: from n with a naive sampler to 1. Nonetheless some parameters may have a major impact on the rejection rate. To clarify it, we studied in Figure the impact of several user- specified parameters the number of motifs to be discovered, the precision of the meshes, the regularization parameter of the ridge score and the nunber of computation cores allowed during the rejection step of the hit-and-run sampler.
Although the number of motifs to be found by SEISM undoubtedly affects the selection cost, we can roughly consider that this relationship is linear. The upper left figure in Figure 6, however, demonstrates that the influence on the overall computing cost superlinear, in line with the exponential growth of the number of distinct selection events one may describe with a fixed mesh size. As result, the post-sclection process quickly becomes intraciable for discovering and test more than few motifs.
We make a similar observation for mesh precision: computation time grows exponentially with the number of bins used to define the meshing. This can be explained by the exponential relationship between the number of bins and the number of different meshes (and thus the rejection rate). Of note, mesh precision has no impact on the selection time, and therefore the computation time entirely explained by the acceptance rate.
We observe that the greater the regularization parameter X, the lower the computation time. This can be explained by detailing its impact on the rejection rate. To understand it, it is necessary to note that the motifs are not selected over Z, but over a less constrained set as described in 3. They are only projected onlo 2 at the end of the whole procedure, to ease their interpretation. The meshes are then defined over a vectorial space, leading to an infinite number of meshes. Compared to a small regularization parameter, a higher favors motifs resulting in a x=*x with a higher norm. With regard to the activation function, such motifs are located closer from the k-mers, and thus from Z. À has then no effect on the nunber of existing meshes, but impacts the number of acceptable ones, in the sense that they have a reasonable probability to be selected. A lower À leads to better selection performances, but to higher number of acceptable meshes, and thus to a lower acceptance rate. We empirically set X = 0.01 to provide a good trade-off.
Finally, the rejection sampling step can be parallelized over several computation cores, which accelerates the whole procedure, as described Section 4.4. As long as the acceptance rate is small enough, using j cores to parallelize the rejection step should roughly divide the computation time by j.
We can clearly identify limitations inherent to the use of the selective inference procedure. Although it is more powerful than the data-split approach, it can not be used in every situation. This latter approach does indeed not include any rejection step, and the only factor influencing its overall computation time is the selection time, only marginally influenced by the aforementioned parameters.
5.4 Using SEISM on empirical data
The experiments on simulated data that we present in Sections 5.2 and 5.3 are useful to analysc the cali- bration, power and computational behavior of SEISM in a controlled environment where the ground truth is known. Here we use an empirical dataset to evaluate two other critical aspects: the robustness of SEISM to the Gaussian assumption, and its ability to select the correct number of filters. We rely on the ChIP-seq dataset from Chatagnon et al. (2015). As part of this study, the authors investigate the mechanisms under- lying cell differentiation and are particularly interested in the retinoic acid receptor (RAR), a transcription
18

---

Published in Transactions on Machine Learning Research (06/2023)
2000
1500
1000
200
Namber of molils be selected
Number of birs the meshing
time
Kidge rig.laaiaatiinnnnnoooooo parameter
Figure 6: Impact of different parameters on the computation time for 100 replicates for the post-selection inference procedure. Upper left: Impact of the number of motifs be discovered Upper right: Impact of the number of bins defining the meshes. Bottom left: Impact the ridge regularization parameter. Bottom right: (Log scale) Impact of the number of threads over which the hil-and-run sampling paral- lelized.
19

---

Published in Transactions on Machine Learning Research (06/2023)
factor. They perform a ChIP-seq experiment, which leads to dataset containing 131,895 sequences of length 500 and their associated phenotypes: the - -oggo of p-values resulting from test to determine whether the sequence is associated with high number of bindings with RAR. We then derive smaller dataset contain- ing only 1,000 sequences, enriched with sequences significantly associated with RAR, order to increase the signal to noise ratio and speed up the computations.
The Gaussian assumption is strongly challenged on classification labels, but can also be questioned for continuous phenotypes. Although this problem is not limited to SEISM, we propose here an approach to assess whether the Gaussian assumption is valid for given dataset (X,z) for the SEISM procedure. This method follows 3 steps:
1. Create N datasets (X, xxyy (0) derived from the original one. The sequences are unchanged, but the labels are randomly permuted versions of y. This permutation ensures that these new datasets are under the null hypothesis, while maintaining the original probability distribution of y.
2. Run the whole SEISM procedure each of those permuted datasets, and collect the p-values.
Draw a Q-Q plot: if the distribution is close to the uniform, then validates the use of the Gaussian model for this datasct.
Our analysis in Section 5.3 suggests that the selective inference version of SEISM would be too computa- tionally intensive and would bring little improvement compared to the data-split version this dataset, and we therefore apply the above procedure with data-split. We use the ridge association score with a penlaty paramter À = 0.1. The resulting empirical distribution of the labels is far from a Gaussian one (see Figure 7), but the p-values obtained on the permuted datasets with the aforementioned methodology are uniformly distributed between and 1, as shown Figure 8, which validates the use of SEISM on this particular dataset.
0.06
Dansty
0.04
0.82
0.00
20
40
80
100
Phenotypes
Figure 7: Empirical probability density of the phenotypes in Chatagnon et al. (2015).
We then select and test four motifs using the same data-split version of SEISM on the non-permuted dataset. The resulting motifs are represented in Table with their respective p-values.
The first motif found- the one with the lowest p-value p-value-recoves recovers a known motif for RAR (Balmer & Blomhoff, 2005), see Figure On the other hand, only the first discovered motif associated with a significant value aligning with the current literature for RAR. This confirms the capability of SEISM infer a posterior the right number of feature in the model.
20

---

Published in Transactions on Machine Learning Research (06/2023)
1.0
data-spiit
0.6
pvelues
0.4
0.0
0.2
0.4
0.6
0.8
10
emplriccl quantlies
Q-Q plot obtained by applying the SEISM procedure permuted versions of Chatagnon et al. (2015).
Figure
Motif 10-3
Motif 2 p= 0.013
A TCA
T CAAETAc
Fee CCc
Motif 3 p= 0.47
Motif 4 = 0.23
ACC
CAAA A GGA
C
SAC
Table Motifs and p-values obtained with SEISM on the dataset from Chatagnon al. (2015).
6 Discussion and future works
We have introduced a procedure to test the association between features learned by a neural network and the outcome predicted by this network. We did so by relying on the post-selection inference framework and formalizing the network training as fcature selection step. Along the way, addressed general problems related to selective inference over composite hypotheses, which has implications beyond testing of features extracted by trained neural networks. In particular to our knowledge, all previous procedures had to work under the assumption that the variance was known. Our strategy to normalize the statistic to make it scale-free could casily transferred to kernelPSI for testing the association of kernels with a trait, or to previous selective inference frameworks for testing groups of variables using sampling strategies (Slim et al., 2019; Reid et al., 2018).
From neural network perspective, SEISM provides principled way to select the number of filters of CNN, with a different Sicettiv--sgnnifiaanee. than the usual prediction-oriented cross validation pro-
TGACCI
Figure 9: Known binding motif for RAR (Balmer & Blomhoff, 2005).
21

---

Published in Transactions on Machine Learning Research (06/2023)
cedures. Through the SEISM procedure, we are also drawing connections between neural networks for biological sequences and two related fields: sequence motif detection, and GWAS.
Sequence motif detection has been a major theme in bioinformatics for the past 30 years and many methods have been proposed to identify motifs that are vver-represented in a set of sequence compared to some control class or background distribution. The earliest CNNs for regulatory genomics Alipanahi et al. (2015); Zhou & Troyanskaya (2015) already exploited the fact that trained convolution filters of the first layer could be interpreted as PWMs, and more recent work have sought to extract PWMs from entire multi-layer trained networks through attribution methods. The selection step of our procedure merely formalizes that training a one-layer CNN is equivalent to selecting a finite set of PWMs that have a maximal association to the outcome for some particular score. This formalization also highlights the specific way by which CNNs with exponential activation functions parameterize the distribution of k-mers at a binding site. Although the PWM returned by most bioinformatics models represents a categorical each letter at each site, trained convolution matrices paramcterize Gaussian distribution. In practice, this difference leads to discrepancies belween the trained convolution filters and PWMs learned using categorical likelihoods- -ie--diiclllngg those offered by databases and often used as ground truth. This observation also suggests alternative sets of constraints for convolution filters-e.g. e.g., cach column of the pointwisc exponential of the filter should belong to the simplex.
to draw
By providing an inference procedure for features extracted by the trained model, our work also connects neural networks for genomic sequences to GWAS. The good predictive performances of these neural networks is often explained by their ability to jointly learn an appropriate data representation and regular prediction function acting on this representation. Nonetheless, the space from which these representations are learned is seldom formalized and to our knowledge the association of the extracted features with the predicted outcome is never tested. GWAS on the other hand relies on hypothesis testing, but commonly relies on relatively simple genomic variants such as single nucleotide polymorphisms (SNPs) or k-mer presence (Jaillard et al., 2018; Roux de Bézicux et al., 2022). Our framework paves the way to GWAS over richer sets of variants, e.g. capluring the presence of entire polymorphic genes through large convolution filters, or the interaction of simpler variants through multilayer or self- attention networks (Avsec et al., 2021a). This will require scaling to entire genomes inputs, and making more complex networks, such multi-laycr CNNs and networks using attention mechanisms, amenable to inference. The most important step in achieving this goal is to formulate the training of these networks as a feature selection problem and formalize the association between these features and the phenotype. The inference framework might then be directly derived from this present work. For instance, we may test motif interactions derived from convolutional-atteniion networks (Ullah & Ben-Hur, 2021), or a (motif, position) couple as selected in (Ditz el al., 2022). Regarding multi-layer CNNs, several strategics are conccivable. One solution would be to build on TFModISco, that aims at extracting motifs summarizing the features captured by a trained multi-layer CNN (in particular, accounting for possible interactions). These extracted motifs could be tested using the SEISM framework: the selection event is the set of simulated phenotypes that would lead to the construction of TFModISco motifs (Shrikumar et al., 2018) within the same meshes. Of note in this strategy, the motifs would not be selected using the same score used as a statistic for testing, but this is not an issue. A second, more integrated possibility would be to test the filters of the first layer within a deeper network. Deeper layers indeed model interactions between the motifs, and with such an architccture the filters of the first layer would then be optimizod while taking into account those interactions. The selection procedure would still be a greedy one, and the architecture of the network would vary from step to step: the first layer would contain more and more filters, but the previously entered filters would be fixed. This selection procedure could then be translated into a selection event, and the inference framework could then be applied accordingly. Finally, we could test the deeper feat ures themselves, but this would only make sense for an appropriate architecture that makes these features interpretable. A simple option would be, for example, to apply a global (or large enough) pooling over a first layer with few filters, and test filters of the second layer that would represent motif combinations. Alternatively, the second layer could be an even simpler set of pairwise interactions between motifs (i.e., special filters with only two non-zero entries). Tn practice, this could be done by optimizing the residual errors for the successive filters of this second layer. Admittedly, a few practical problems may arise. First, the hit-and-run sampler requires the selection method to be stable, that is, running the selection method twice on the same input will lead to the same selection on features. This property is required to guarantee the
22

---

Published in Transactions on Machine Learning Research (06/2023)
theoretical convergence of the the algorithm but may not be necessary in practice. Second, some attention may be required to avoid the that computational cost become prohibitive, in particular depending the regularity properties of the selection event leading to a higher rejection probability or to a higher number of replicates required. Granted that these technical challenges can be addressed, we are confident that exlending SEISM to more general netwsrks and corresponding features will benefit. both the fields currently using these networks such as regulatory genomics- GWAS.
23

---

Published in Transactions on Machine Learning Research (06/2023)
7 Acknowledgements
This work has been supported by ANR grants (FAST-BIG project ANR-IT-CE23-001-00 and PIECES project ANR-20-CE45-0011) and was performed using the computation facilities of the LBBE/PRABI.
We thank François Gindraud, Jean-Philippe Rasigade, Lotfi Slim and Dexiong Chen for the insightful dis- cussions and support.
References
Babak Alipanahi, Andrew Delong, Matthew T. Weirauch, and Brendan Frey. Predicting the sequence specificities of DNA- and RNA-binding proteins by deep learning. Nature Biolechnology, 33(8):831 838, July 2015, ISSN 1087-0156. doi: 10.1038/nbt.3300, URL http://åx.ddea doxx www.xxooi.oog/1000 10. 1ww.//0t..nn3. 3300.
Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré. Gradient flows: in metric spaces and in the space of probability measures. Springer Science & Business Media, 2005.
Ziga Avsec, Vikram Agarwal, Daniel Visentin, Joseph R. Ledsam, Agnieszka Grabska-Barwinska. Kyle R. Taylor, Yannis Assacl, John Jumper, Pushmeet Kohli, and David R. Kelley. Effective gene expres- sion prediction from sequence by integrating long-range interactions. Nature Melhods, 18(10):1196 1203, Oct 2021a. ISSN 1548-7105. doi: 10.1038/s41550 022-0-22222 URL https://doi http:://ddo..org10...t org/10. 941592-021-01252-xx
Ziga Avsec, Melanie Weilert, Avanti Shrikumar, Sabrina Krueger, Amr Alexandari, Khyati Dalal, Robin Fropf, Charles McAnany, Julien Gagneur, Anshul Kundaje, and Julia Zeitlinger. Baso-rosolution models of transcription-facten binding reveal soft motif syntax. Nat Genet, 53(3):354-38) February 2021b.
Timothy L Bailey. STREME: accurate and versatile sequence motif discovery. 37(18):2834 2840, 2021. ISSN 1367-4803, 1460-2059. doi: 2a.109//biiiffommn hatrrmaattissbbbb2a333 btab203. URL https://academic. oup. com/
Timothy L. Bailey, Nadya Williams, Chris Misleh, and Wilfred W. Li. MEME: discovering and analyzing
DNA and protein sequence motifs. 34:W369 W373, 2006. ISSN 0305-1048. doi:
URL
/ttps://www.i ncbi nlm. mih, www.nnii...o.........
J. E. Balmer and R. Blomhoff. robust characterization of retinoic acid response elements based on comparison of sites in three species. The Journal of Steroid Biochemistry and Molecular Biology, 96(5): 347-354, September 2005. ISSN 0960-0760. doi: 10. 10.1016/j.jsbmb..
Leon Bungert and Philipp Wacker. The lion in the attic a resolution of the borel kolmogoro paradox, 2022. URL bttt:://arxiy http//www.vi......aa.ooo/. wwwp:/r////b///////// httpaabb///00..../
Amandine Chatagnon, Philippe Veber, Valérie Morin, Justin Bedo, Gérard Triqueneaux, Marie Sémon, Vincent Laudet, Florence d' Alché Buc, and Gérard Benoit. RAR/RXR binding dynamics distinguish pluripotency from differentiation associated cis-regulatory elements. Nucleic Acids Research, 43(10))4833- 4854, May 2015. ISSN 1362-4962. doi: 10.1093/nar/gkvv//
Dexiong Chen, Laurent Jacob, and Julien Mairal. Biological sequence modeling with convolutional kernel networks. 2017. doi: 10.1101/217257. URL http://biorxiv..o htt.://iwi..rr//oooo/////////... 1101/217277.
Jonas C. Ditz, Bernhard Reuter, and Nico Pfeifer. Convolutional motif kernel networks, 2022. URL http: 02272.
ENCODE Project Consortium. The ENCODE (ENCyclopedia of DNA elements) project. 306(669)):636- 640, 2004. ISSN 1095-9203. doi: 10.1126 htts//ccecc..o00 105136.
Shobhit Gupta, John A. Stamatoyanmopoulos, Timothy L. Bailey, and William Stafford Noble. Quantifying similarity between motifs. 8(2):R24, 2007. ISSN 1474-760X. doi: 10.1186//b-2007--8 2-r24. URL https: / //doi.org/10. 11w.:/1b-/bb-00------
24

---

Published in Transactions on Machine Learning Research (06/2023)
R. Harr, M Häggström, and P Gustafsson. Search algorithm for pattern match analysis of mucleic acid sequences. Nucleic Acids Res, 11 (999229943220) 29477 May 1983.
Trevor Hastie, Robert Tibshirani, and Martin Wainwright. Statistical learning with sparsity, Monographs on stalistics and applied probability, 143:143, 2015.
Magali Jaillard, Leandro Lima, Maud Tournoud, Pierre Mahé, Alex van Belkum, Vincent. Lacroix, and Laurent Jacob. A fast and agnostic method for bacterial geneme-wide association studies: Bridging the gap between k-mers and genetic events. PLOS Genetics, 14(11):1-28, 11 2018. doi: /0.187//journal.pgom 1007758. URL https://doi: wwwpeiiorgg/... journal. 1007758.
Arttu Jolma, Jian Yan, Thomas Whitington, Jarkko Toivonen, Kazuhiro R. Nitta, Pasi Rastas, Ekaterina Morgunova, Martin Enge, Mikko Taipale, Gonghong Wei, Kimmo Palin, Juan M. Vaquerizas, Renand Vincentelli, Nicholas M. Luscombe, Timothy R. Hughos, Patrick Lemaire, Esko Ukkonen, Teemu Kivioja, and Jussi Taipale. DNA-binding specificities of human transcription factors. 152 339- 2013. ISSN 00928674. doi: 00.1016/j.cell.20122 22.00.0.. 12.009. URL https://1inkkngguut wew.eegehhu.beeeeeer.eeeeee..eee elsevier http:///rr...crrteeee///////// 50092867412014961.
David R Kelley, Yakir Reshef, Maxwell Bileschi, David Belanger, Cory Y McLean, and Jasper Snoek. Soquential regulatory activity prediction across chromosomes with convolutional neural networks. Genome Res, 28(5):739 750, March 2018.
Peter K. Koo and Sean R. Eddy. Representation learning of genomic sequence motifs with convolutional neural networks. 15(12):e 2019. ISSN 1553-7358. doi: p0. 1007560 URL 1007560. Publisher:
attps://journals. .po... www.ppll.........pp.pp.p Public Library of Science.
Alex Krizhevsky, Πyya Sutskever, and Geoffrey Hinton. Imagenct classification with deep convolutional neural networks. In F. Pereira, C.J. Burges, L. Bottou, and K.Q. Weinberger (eds.), Advances in Newral Information Processing Systems, volume 25. Curran Associates, Inc., 2012. URI. https://proceeding pdf.
Yann LeCun and Yoshua Bengio. Convolutional Networks for Images, Speech, and Time Series, pp. 255-258.
MIT Press, Cambridge, MA, USA, 1998. ISBN 0262511029.
Jason D. Lee, Dennis L. Sun, Yuekai Sun, and Jonathan E. Taylor. Exact post-selection inference, with application to the lasso. 44(3):907-927, 2016. ISSN 0090-5364. doi: 10.1214/15-A0S1371 URL http: //arxiv.. wwwxiv..rg//ba...... harg//bb//// 11/111111111 6238.
M.A. Lifshits. On the absolute continuity of distributions of functionals of random processes. Theory of Probability & Its Applications, 27(3):600-66 -07, 1983.
Joshua R. Loftus and Jonathan F. Taylor. Selective inference regression models with groups of variables. arXiv e-prints, art. 1rX1vv1511.01478, November 2015.
Christoph Molnar. Interpretable Machine Learning. 2 edition, 2022. URL https://christophi ..n..ggii github. Ao/interpreeablemm1--2 mol.rretttalemmmbook..
Gherman Novakovsky, Nick Dexter, Maxwell W. Libbrecht, Wyeth W. Wasserman, and Sara Mostafavi. Obtaining genetics insights from deep learning via explainable artificial intelligence. Nature Reviews Genetics, Oct 2022. ISSN 1471-0064. doi: (011///1157--22-000555 URL https://doi.org/// 0098/841576-022-0000 00532-2.
Gabrielle Ras, Ning Xie, Marcel Gerven, and Derek Doran. Explainable deep learning: A field guide for the uninitiated. J. Artif. Int. Res., 73, may 2022. ISSN 1076-9757. doi: 10.1613/jair.1.3300 URL /ttps:://doii.rg/1....
13200.
Stephen Reid and Robert Tibshirani. Sparse regression and marginal testing using cluster prototypes. 2013.
URL http:://rxiv http:://aardvv..//aa//////////// org/abs///0/ 00334.
25

---

Published in Transactions on Machine Learning Research (06/2023)
Stephen Reid, Jonathan Taylor, and Robert Tibshirani. A general framework for estimation and inference from clusters of features. 2015. URL http://arxiv.. abs/1511.111...//11/11//11111111
Stephen Reid, Jonathan Taylor, and Robert. Tibshirani. A general framework for estimation and inference from clusters of features. Journal of the American Statistical Association, 113(521):280 293, 2018. doi: 10.1080/001224459.200 1246668.
URL https: attps:///oi... http://///oooo////////o wwg/oo..//..o.......o 108000162199999 01000/00000001010000000011114111 2016. 1246368.
Guillaume Rizk, Dominique Lavenier, and Rayan Chikhi. DSK: k-mer counting with very low memory usage.
29(5):652-655 2013. ISSN 1367-4811. doi: 10. http://innoomaaiisssttt/tto
Hector Roux de Bézieux, Leandro Lima, Fanny Perraudeau, Arnaud Mary, Sandrine Dudoit, and Laurent Jacob, CALDERA: finding all significant de Bruijn subgraphs for bacterial GWAS. Bioinformatics, 38: 136-i44, 06 2022. ISSN 1367-4803. doi: 10.1093/bion hita///oooormmmttee/raa22...... URL https://doi. org/
1098/bioinffrrmtiit/bbaaaa
Shuxiang Ruan and Gary D. Stormo. Inherent limitations of probabilistic models for protein-DNY binding specificity. 13(():e1005638, 2017. ISSN 5553-7358. doi: 10.
1005638 URL https:
//journals, petooralssppls......... http:///cccpppp///aaaaaeeaaeee.. http:///d.....1///rrrrrraoao
Publisher: Public
Library of Science.
T D Schneider and R M Stephens. Sequence logos: new way display consensus sequences. Nucleic Acids Res, 10000)0007,- October 1990.
Avanti Shrikumar, Katherine Tian, Ziga Avsec, Anna Shcherbina, Abhimanyu Banerjee, Mahfuza Sharmin, Surag Nair, and Anshul Kundaje. Technical note on transcription factor motif discovery from importance scores (tf-modisco) version 0.5.6.5, 2018. URL https://ent //arxiv org/abs/1811 2abbb//11100016.....111111111111 .00166.
Lotfi Slim, Clément Chatelain, Chloe-Agathe Azencott, and Jean- Phin-Philippe Vert. kernelPSI: a post-sclection inference framework for nonlinear variable selection. In Kamalika Chaudhuri and Ruslan Salakhutdinov (eds.), Proceedings of the 36th International Conference Machine Learning, volume 97 of Proceedings of Muchine Learning Research, pp. press/v97/sll http://////mmmm...///// html.
PMLR, 09 15 Jun 2019. URL httpss// mlr.
Mahdi Soltanolkotabi, Adel Javanmard, and Jason D. Lee. Theoretical insights into the optimization land- scape of overpaaramcterized shallow neural networks. IEEE Transactions on Information Theory, 65))): 742 769, 2019. doi: 10.1109//TI 202.2....225555
Le Song, Alex Smola, Arthur Gretton, Justin Bedo, and Karsten Borgwardt. Feature selection via dependence maximization. Journal of Machine Learning Research, 13(47):1393- 1434, 2012. URL http://jmlz: org/ papers/v13/song122 html.
Jonathan Taylor and Robert J. Tibshirani. Statistical learning and selective inference. Proceedings of the Nalional Academy of Sciences, 112(25):7629 7634, 2015. doi: 10.1073 1t773//mas.1005531111 URL https: //www. pnas. org/ /o////o///a adoo//ab///i a0..////.... 1073 //aas////asaa./ 1507583112.
Jonathan Taylor, Richard Lockhart, Robert Tibshirani, and Ryan J Tibshirani. Exact post-selection inference for forward stepwise and least angle regression. pp. 32, 2014.
V.S. Tsircison. The density of the distribution of the maximum of gaussian process. Theory of Probability & Its Applications, 20(4):847 856, 1976.
Fahad Ullah and Asa Ben-Hur. A self-attention model for inferring cooperalivity between regulatory features. 1971 2021. ISSN 0008-1048. doi: 200093/aar/gabb// URL https://doi.org gkab349.
Ashish Vaswani, Noam Shazoer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaisor, and Illia Polosukhin. Attention is all you need. Advances in neural information processing systems, 30, 2017.
26

---

Published in Transactions on Machine Learning Research (06/2023)
P M Visscher, N R Wray, Q Zhang, P Sklar, M McCarthy, M A Brown, and J Yang. 10 years of GWAS discovery: Biology, function, and translation. Am J Hum Genet, 101(1):5-22 22, July 2017. doi: :...11//).jjgg.211.0600000 URL https://www.nctt nlm alm.nih. www..ii..ooo...
Larry Wasserman and Kathryn Roeder. High-dimensional variable selection. 37(5):2178-22 2000)22211222 2009. ISSN 0090-5364, 2168-8966.
doi: 10. 10.1214/08-AOS644 Publisher: Institute of Mathematical Statistics.
Ronald L. Wasserstein and Nicole A. Lazar. The ASA statement on p-values: Context, process, and purpose. The American Statistician, 70(2):129-133, 2016. doi: 00000/00003335220150 1154111..188...1 URL https:///d8 2tg/1...10800/0033355.22222 2016. 1154108.
Makoto Yamada, Yuta Umezu, Kenji Fukumizu, and Ichiro Takeuchi. Post selection inference with kernels. 2018.
Jian Zhou and Olga G Troyanskaya. Predicting effects of noncoding variants with deep learning-based sequence model. Nat Methods, 12:931 4, 2015 Oct 2015. ISSN 1548-7105. doi: 10.1038/ 10.1038/nmeth.B35 m54th
27

---

Published in Transactions on Machine Learning Research (06/2023)
Supplemental Materials of "Neural Networks beyond explainability: Selective inference for sequence motifs"
A Tuning the activation bandwidth hyperparameter
The data representation u0 x depends hyperparameter controlling the bandwidth of the gaussian non-linearity (Equation 3): exp 2 - Assuming that the positions are independant, we know that
the expected value of the distance between a motif and k-mer with length proportional to
In order to get an activation that does not depend on the length of the motifs, we simply set to be proportional to Vk. From empirical tests, we set w = v0.9-k to achieve good selection results by choosing 2
the motif that maximizes the association score among a set of possible lengths.
B Disintegration of the selection event given by sequence motifs
In this section we consider the selection event:
Ecomt. (Z) (Z). y/ E E Vi fl.......) /) arg magaax s(z, P.y') 21).
(S1)
#E2
given by the sequence of selected motifs Z = , zn). a Gaussian distribution on E.
We denote by the law of as given by Eq. (9),
A first remark on the uniqueness of the selection
Consider the mapping 3 + 2a given by (3y) Z where Z = (21,..., 2g) is the sequence of motifs such that y' E Eeont. (Z). It is not clear that π is well defined as a same y' may lead to the selection of at least two different motifs sequences Z and Z'. As a first remark, we can see that the set of problematic y is exactly
: U Ecoat. (Z) Econt. (Z)).. Z
When one assumes that Z = (zi...., 2g) is unique, one implicitly assumes that μ(P) 0. For sufficiently regular scores, this is however the case. For sake of readability, we will not comprehensively study this issue but we will present an argument for the scores HSIC and s ridge. In this case, we can circumvent this difficulty considering the Gaussian random field
H (xp= . for (HSIC) and 2 ((1p==x11 1p An)
1/2 112) x=xxxxx,, . for (Ridge)
indexed by where y is distributed with respect to a multivariate Gaussian distribution Eq. (9). Its autocovariance function is given by (2,z") + 02(gz. xx from Eq. (9) (one has to multiply by (Ilp=:X11 An) -1/2 -n)-1/2 (Ilp=''x112 + Xn) for the Ridge). The score is just the largest norm of this Gaussian random field. It is well established in theory of Gaussian random fields that the law of this maximum is regular and the argument maximum is unique. The interested reader may consult the pioncering work of Tsirelson (Tsirelson, 1976) and Lifshits (Lifshits, 1983). In Tsirelson's theorem, the parameter set is count- able. This says that the same result holds true for separable bounded Gaussian processes, since in this case, the distribution of the supremum coincides a.s. with the one of the supremum on some countable nonrandom set. To avoid a cumbersome presentation, we will assume that almost surely the selected sequence motifs Z (Z1,..., 2g) is uniquely defined, hence is well defined.

---

Published in Transactions on Machine Learning Research (06/2023)
The disintegration steps
To sample conditionally on (S1), need to consider the conditional law with respect this event. We will denote this law by (zz, it depends only μ4, 2 and æ. This law is described by the theorem of disintegration, see for instance (Ambrosio et al., 2003, Theorem 5.3.1). Denote v the pushforward mcasure of μ by denoted by v = и#,, a probability measure on the set 2a of Z. By the disintegration theorem, there exists a v-almost everywhere uniquely determined Borel family of probability measures pzz (the though-after conditional distributions) such that
Supported by Econt. (Z); uz (0 =0 for v-almost every Z;
Expectation of the conditional expectation is the expectation: It holds that, for every Borel test. map f - [0, too),
| fdu fdu= (/-112) fdpz fduz)du(z) du(Z).
(S2)
where one can remark that π T Econt. (Z) by definition of x. Let us comment on this result regarding our purposes. First, we have mentioned that we known that the support Econt. (Z) is included in some subspace, say S, defined by the first order conditions. Second, although strategy the subspace S to draw points on the support Econt. (Z) (viewed subset of the same Hausdorff dimension as the subspace S), it is not clear at all what should be the density of uzz Indeed, the family of probability measures Pz is the unique family that satisfies Eq. (S2). It implies that a. measure μzz depends on the others measures jAz' and this dependency is geometrically given by the (picce-wise) topological sub- manifold given by the function z + x from to E.
rejection sampling
From a practical view point, we tried various law for pz such as the uniform, or rejection sampling based the Gaussian distribution (9), but none of them matched the condition (S2). In the next subsection, we recall toy example: the disintegration of the uniform measure on the sphere is not the uniform measure. Even in this simple geometrical example, the calculus of the conditional law might be seen as tedious. We believe that the calculus of jzz is somchow out reach for our purposes and our analysis with selection events defined by meshes more suited.
A toy example on the sphere
Let $ be the 2-sphere embedded in the B-Euclidean space. Let μ be the uniform measure on the sphere S. Let (So A E [0, ])} be family of sub-spaces of co-dimension 1 (hyper-plancs) sharing Spanf(0. O, 1)} (say the north pole) as revolution axis parameterized by 0. The parameter 0 can be interpreted as the longitude.
So
Figure S1: For 0= So is the light red plan, the conditional measure doo(e) is depicted with red area and is proportional to (ii( which is not the uniform measure.
2

---

Published in Transactions on Machine Learning Research (06/2023)
Let & be the function that maps point to its longitude modulo л. By spherical symmetries, the pushforward
measure # is the uniform measure on [), so that dv(0) (1 (1/x)] //00. Condition (S2) (the Ihs of the equality below) is given by the coordinate integration system (the rhs) in:
[ (// sdwm) -
2T ((0, S(,000-(((A SS²
sin Φ dØ d0, )d)
4x
where & is the latitude. Note that (0) = SnS, and it is in bijection with [0, 2T) using the mapping that to point maps its latitude. Using this representation, it is not hard to see that the uniform mcasure on (=) is given by (0.2x)) while the above equality shows that the conditional measure μg of the uniform measure on the sphere has density (1/4) sin p01 see Figure S1. It proves that the disintegration of the uniform mcasure on the sphere not the uniform mcasure, but rather distribution that will put few mass around the poles and large mass around the equator.
C Proof of Proposition 4.4
Consider the orthogonal decompositon
R
where R is the span of 4p x T is the span of μ (orthogonal to R by Proposition 4.1), and S such that the equality holds. Consider y E & and its othogonal decomposition y = = + le where 14/l|al22 unit norm vector that spans T. Let T > and note that is enough prove that
s Z Y) Y
8 (Z,
PA
$T} Pof
S
2 ST).
where Y is a random variable with the same distribution le + ae (resp. σε) on the probability space defined by Pa (resp. Po). Not that the event decomposed
y s(Z,y) 1
{T ((),,)) s(2,r) S (( + ||r|| + [s|]?]]
By othogonality, note that (u(r, 8) Lo(r, 8) and this law is centered Gaussian multivariate law. We deduce that the aforementioned probabilities are of the form
{ s( sz,Y)
ST] = Lo" (o)())p(()d) 21
PA
where
w/ (1) Po Pofs(Z,r) S T(( 12 + [|s|
(t) = exp (-(t - -ee) + exp (-(t + /=)²/2)
with Me (e. |||| Note that wa (0) + (0, 1) is an increasing continuous function. is an increasing homeomorphism and the Fubini's equality yields
p 3 (Z.Y) 1
[ wo(D)pp())d
S
du (u
1 1,m, (w)s du up dt
Ju
(t)dddu
No ((a)
3

---

Published in Transactions on Machine Learning Research (06/2023)
By Anderson's theorem, the measure of the interval [- (u), wy (a)) for the centered Gaussian density is greater than the one for a non-centered Gaussian density with the same variance. As result, we deduce that
2 S po(t)dt,
Px(1)dt
'(u)
w (w)
which achieves the proof.


Reviews:
[{"review_id": "op5m5oLuWe", "reviewer": "Reviewer_XYhN", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "This paper focuses on the formalism of gradient-based interpretability for identifying DNA sequence motifs. Given this topic I don't think I am the right reviewer for analysing the strength of this paper.\n\nAs far as I can tell, the paper does not claim to significantly advance the performance over current state of the art, by any metric.\n\nInstead, the argument is that the formalism is defined and would be useful in extending the method further. It is difficult for me to evaluate this analysis. The paper does not articulate explicitly what is possible with the new formalism or show actual results on extending it.\n\nAs a result, the reader is left to evaluate the utility of the formalism on its own. However, I'm afraid this is beyond my area of expertise to review. I am not an expert in neural network interpretability or on methods for identifying DNA sequence motifs.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper formalizes an approach to interpretability for sequence motif identification for neural networks. The paper does not seem to add any significant methodological advances over prior methods.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Actually showing an extension that the proposed formalism enables would go a long way to showing a broader audience the utility of the method.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Ums5QhKOW0", "reviewer": "Reviewer_rRMK", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\nThe work proposes a holistic framework for extraction and association of features with learning targets.\n\nThe proposed approach for post-selection inference can deal with continuous feature space\n\nThe authors provide theoretically guided parameterisation of null distribution, which enables sampling for post-selection inference procedure. They show that under the assumed null hypotheses and the Gaussian outcome model, the parameterisation of the distribution is invariant to scaling and translation.\n\nThe proposed method achieves competitive performance on motif discovery tasks. Experiments on synthetic data validate the proposed inference procedure  \n\nThe framework may outperform data-split ]earning and inference strategy in small data regime\n\nWeaknesses:\n\nThe major weakness of the proposed approach is that it is not scalable. Parameters such as the size of the input alphabet, size and number of features, granularity of the mesh structure each cause super-linear/exponential increase in the computational cost of the inference procedure. Even the experiments in the paper are rather small scale. It is therefore not obvious if the proposed framework is broadly applicable on a real scale compared to what the experiments demonstrate in the paper.\n\nAlthough post-selection inference is the main contribution of this paper, the work falls short of demonstrating its utility in associating/inferring the right number of features, which bears significance for any real application\n\nThe framework assumes Gaussian noise model for targets, but it is not clear if it in general is a valid assumption for real application scenarios of the framework and what could be the implications when the assumption does not hold. Furthermore, the Gaussian activation units used in the convolutional neural network do not make an obvious choice for modelling one-hot encoded categorical input.\n\nThe framework is proposed to be generic and flexible, but it is only applied to one type of data. In regards to flexibility no results or variations in the components of the framework have been reported.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "\nThis work proposes SEISM, a framework for extracting features/motifs from biological sequences and (statistically) associating them with target phenotypes. For feature extraction, the framework employs a single-layer (convolutional) neural network with radial basis activation and max pooling. The framework furthermore proposes a hypothesis testing procedure to establish (the statistical significance of) associations between the extracted features (z_j) and prediction targets/phenotypes y. \n\nThe authors define an association score as their optimisation objective, which is a function of feature and prediction targets, and it is derived from the solution of linear ridge regression. They propose an iterative greedy approach for learning features such that in every iteration, a new feature is learnt to model the residuals targets from previously learnt features. The size (k) and number of features or motifs (q) are hyper-parameters of the the feature extraction module. \n\nOnce the features are learnt, the authors propose a post-selection inference step to estimate p-values for establishing the statistical significance of association between each feature and target vector y. To that end, authors propose hypothesis testing which assumes a null hypothesis that is defined to have zero association score between a motif and the vector of target outcomes, implying orthogonality between the two entities. To estimate p-values for hypothesis testing, the authors propose to empirically compute the null distribution of association scores. They partition the feature space into a mesh structure so that the feature space spanned by the mesh partitions of extracted motifs can be sampled for the estimation of the null distribution of association scores. To simulate alternative target outcomes y_prime which lead to the selection of other motifs from the mesh partitions of learnt motifs, the authors propose a hit-and-run MCMC sampler. For computing the null distribution and corresponding test statistic, the authors put forth two definitions of null hypothesis based on two varying representations of mesh partitions. They also present theoretical analysis on the approximation of the null distribution through sampling from a global distribution.\n\nIn numerical analysis, SEISM is benchmarked against three other methodologies on a motif discovery task. The authors use a publicly available set of datasets to compare the accuracies/alignment of motifs discovered by different methodologies with the known ground-truth. They compute so called Tomtom score as a measure of accuracy. They apply thresholds to the best Tomtom scores achieved for each dataset by each of the compared methods to compute their success rates across all datasets. The authors find that for low threshold levels, their proposed method is comparable or better than the other approaches, whereas for higher levels of threshold, two of the compared methodologies tend to achieve better success rates in general.\n  \nNext the authors present calibration experiments using synthetic data, where they show that SEISM is on-par with data-split strategy for learning and inference in terms of empirically recovering the underlying null distribution of p-values. In further experiments on small synthetic datasets where the outcomes were generated under alternative hypothesis, SEISM is shown to be more salient in deviating from the null distribution than data-split strategy. \n\nLast set of experiments demonstrate non-linear/exponential scaling behaviour of SEISM with respect to various hyper-parameters.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Critical:\n\nI would like to understand why to define a selection region for pose-selection inference, the authors propose to partition the feature space into a rigid mesh structure instead of for instance defining the selection regions to be within a certain distance from extracted features. Wouldn’t that lead to more uniform sampling from the space surrounding the extracted features?\n\nIt seems that in the greedy optimisation proposed by the authors, the order in which features are learnt would influence their relevance to the outcome. Can the authors show if there is in general any correlation between the order of extracted features and their associated p-values? In case that the order of features is found to have an influence on the association between features and outcome, would it not undermine the relevance of (computationally costly) post-selection inference procedure for assessing the associations between individual features and outcome?\n\nTo substantiate the claim that the proposed framework is generic and flexible, the authors should perform numerical experiments on different types of data. They should also demonstrate the scalability of their framework beyond the small scale experiments (2 or 5 features with coarse meshing) in the paper.  \n\nAs post-selection inference is the main contribution of the paper, it feels imperative that the method is evaluated for its performance on extraction/selection of the “right” number of latent features from observed data. This can be demonstrated through (generated) data with known ground-truth.\n\nNice to have:\n\nThrough experiments, authors can also show how post-selection inference works for features learnt by different methods (for different input and targets). They can perform experiments with varying CNNs\n\nIt would be nice to see experiments that show how the bin size affects the performance of post-selection inference, especially when the input cardinality, size and/or number of features increase.\n\nFor feature selection/importance assessment, it would be nice to contrast the proposed method with existing approaches e.g., the Shapley Value.\n", "broader_impact_concerns": "No. concerns in this regard", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "kUNO5lh4vd", "reviewer": "Reviewer_FHuH", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\nA code implementation is provided\n\nProposed method performs comparably to other baselines in the motif discovery benchmark\n\nOther questions/comments:\n\n“For the sake of simplicity, we choose to restrict this presentation to simple one-layer CNNs and sequence motifs. The procedure that we introduce, however, is by no means limited to this framework, and could be applied to any of the more expressive features proposed in the explainable machine learning literature.” -- But the authors mention the one-layer CNNs are a special straightforward case for providing interpretable results since the trained filters can be seen as position weight matrices. What about multilayer CNNs where interpretability is more challenging?\n\n“The user must also specify the number of motifs to find, as well as a parameter controlling the meshing of the motif space, that is the precision with which the found motifs will be tested.” -- How sensitive are the results to the choice of these hyperparameters? How should the number of motifs be determined a priori?\n\n“To that end, we enumerates the k-mers contained in X using the DSK software (Rizk et al., 2013) and compute their scores s(·, ·)” -- How is k set?\n\n“In order to assess the statistical validity and of the SEISM procedure with the different strategies, we simulate datasets under the null hypothesis. To that end, we draw one sequence motif z˜ with length k = 8 for each simulated dataset using a uniform distribution on Z restricted to motifs with an information level fixed at 10 bits” -- Why only analyze a sequence with length k = 8?\n\nIt is unclear to me what are the advantage/disadvantage of the “de novo” motif discovery methods (STREME, MEME) compared to the proposed SEISM method (their performances in Figure 4 are very similar) \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper describes an interpretable method called SEISM for detecting and testing the associations between sequence motif features learned by a convolutional neural network and the outcome phenotypes predicted by the CNN. The method uses a one-layer CNN as an interpretable sequence motif discovery tool, and leverages a post-selection inference framework to select over a continuous set of features. The motif discovery performance of SEISM is benchmarked against de novo motif discovery baselines (STREME, MEME, CKN-seq), and shows comparable performance. ", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Some minor comments:\nFigure 4 right, the error bars for the different methods are overlapping and hard to read\n", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper has gone through a thorough review, discussion and the authors have already updated the paper accordingly. The paper has some issues about general applicability and scalability but is valid and therefore meets the TMLR evaluation criteria. 

The authors are encouraged to add a more realistic test case for the final version.

Minor comment: 

"Over the past decade, neural networks have been successful at making predictions from biological sequences, especially in the context of regulatory genomics." Arguably, there has been more activity in the protein space. So suggest to remove ", especially in the context of regulatory genomics". 
Decision: Accept (poster)