Martin Drance

Fleur Mougin

Akka Zemmari

Gayo Diallo

Martin Drancé

Pre-Trained Embeddings for Enhancing Multi-Hop Reasoning

Knowledge graphs are an efficient way to represent heterogeneous data from multiple sources or disciplines by utilizing nodes and their relations. Nevertheless, they are frequently incomplete in terms of the subject they represent. Link prediction methods are used to discover additional links (or even to create new ones) between entities present in the Knowledge Graph (KG). In order to achieve this, multi-hop reasoning models have demonstrated good predictive performance and the ability to generate interpretable decisions, thereby enabling their application in high-stakes domains such as finance and public health. A multi-hop reasoning model usually has two tasks: 1) construct an accurate representation of the entities and relationships of the KG; 2) use these representations to explore the reasoning paths in the KG that support the newly predicted links. In this paper, we investigate how the performance of a multi-hop reasoning model changes when using pre-trained embeddings for the KG's nodes and relations. The experiments conducted on three benchmark datasets, respectively WN18RR, NELL-995 and FB15K-237, suggest that using pre-trained embeddings improves: (i) the predictive performance of multi-hop reasoning models for all three datasets, (ii) the number of newly predicted links, and (iii) the quality of paths used as explanations.

Introduction

A knowledge graph (KG) is a large data structure that use a graph database to describe real-world entities and their relationships. KGs represent entities and their relationships in the form of a schema, allowing entities to be linked together, and can aggregate information from various data sources and domains [START_REF] Paulheim | Knowledge graph refinement: A survey of approaches and evaluation methods[END_REF]. The facts stored in a KG are in the form of triples, i.e. two nodes connected by a direct and typed link, such as (Tim Berners-Lee →wonAward →Turing Award). As KGs are highly incomplete [START_REF] Min | Distant supervision for relation extraction with an incomplete knowledge base[END_REF], link prediction methods can be used to infer missing links between pairs of nodes, creating and updating the representation of entities and relations as embeddings.

Following the TransE [START_REF] Bordes | Translating embeddings for modeling multi-relational data[END_REF] effort, knowledge graph embeddings (KGE) models capable of learning good representations (or embeddings) of entities and relations in a KG have been proposed. These embeddings can then be used to infer new links by evaluating the plausibility of new triples. Although they have achieved state-of-the-art results for link prediction, these methods are not transparent in nature and their use is problematic when a clear understanding of a prediction is required.

Recently, multi-hop reasoning models (referred as MHR or walk-based methods) [START_REF] Das | Go for a walk and arrive at the answer: Reasoning over paths in knowledge bases using reinforcement learning[END_REF][START_REF] Liu | Neural multi-hop reasoning with logical rules on biomedical knowledge graphs[END_REF]Lin et al., 2018;[START_REF] Lei | Learning collaborative agents with rule guidance for knowledge graph reasoning[END_REF][START_REF] Safavi | Evaluating the Calibration of Knowledge Graph Embeddings for Trustworthy Link Prediction[END_REF] have proposed to address this issue by searching for reasoning paths in KGs using reinforcement learning (RL) to learn a policy, which is then used to predict an answer to a query (entity 1 , relation 1 , ?). As introduced in [START_REF] Das | Go for a walk and arrive at the answer: Reasoning over paths in knowledge bases using reinforcement learning[END_REF], MHR models have the strong advantage of producing directly interpretable predictions. Indeed, the paths used during inference are used to understand what are the existing links between the query and the predicted entity. Interpretability in machine learning is a growing challenge for high-stakes decisions [START_REF] Rudin | Interpretable machine learning: Fundamental principles and 10 grand challenges[END_REF] and in the context of link prediction, transparency is a desirable attribute as KGs are particularly good at storing highly symbolic concepts and relationships.

Despite the encouraging results of MHR models recently, they still struggle with the sparsity and large size of KGs. Since there is no notion of a correct path to reach a correct answer, link prediction can be performed by the agent using an incorrect path but still reach a valid entity in the graph. Also, due to the incompleteness of KGs, true predictions that are missing in the training data won't be rewarded. In both cases, these predictions will result in the alteration of entities and relations representations based on the use of false information.

In this paper, we investigate and report how using pretrained embeddings for nodes and relations can improve the performances of a multi-hop reasoning model and discuss our findings. This improvement is analyzed in terms of link prediction scores as well as the quality of the paths generated as explanations. The use of high-quality pre-trained embeddings is intended to facilitate the training process of MHR models by bringing prior knowledge in the KG representations. In contrast to previous work and state-of-the-art models, we exclude from the MHR model the requirement to learn entity and relationship representations from scratch. These representations are first learned by the KGE method ConvE [START_REF] Dettmers | Convolutional 2d knowledge graph embeddings[END_REF] and then used as is or finetuned by the MHR model. We conduct our experiments on the three benchmark datasets WN18RR [START_REF] Dettmers | Convolutional 2d knowledge graph embeddings[END_REF], NELL-995 [START_REF] Das | Go for a walk and arrive at the answer: Reasoning over paths in knowledge bases using reinforcement learning[END_REF] and FB15K-237 [START_REF] Toutanova | Representing text for joint embedding of text and knowledge bases[END_REF]. The main contributions of our work are the following: (i) we show experimentally that pre-trained embeddings for entities and relations improve the performance of MHR model for three benchmark datasets WN18RR, (ii) we show that the use of pre-trained embeddings allows the model to find at least one reasoning path for more queries, and (iii) pre-trained embeddings increases the overall quality of the reasoning paths used as explanations.

The next section outlines related work in the area of MHR models for link prediction and KG embedding methods. Then, in section 3, we describe our proposed approach to improving multi-hop reasoning 1 . In section 4, we report its performance evaluation before discussing and interpreting it in section 5.

Background

In this section, we first define the properties of KGs and justify their use. Then, we outline how the main KGE models work. Finally, we outline related work in the area of multihop reasoning.

Knowledge Graphs

Given a set of entities E and a set of binary relations R, a KG G ⊂ E × R × E represents the collection of facts expressed as triples (h, r, t) with h, t ∈ E and r ∈ R. All relations r are directed and define the role of entities in each triple, with the left node h being the head (or source) and the right node t being the tail (or target) of the triple. All entities e and relations r belong to specific categories or types, clearly defining what type of triples can be found in the KG. For example, if the KG defines that the relation worksF or should link a node of type person to a node of type organization, it can then be inferred from (T im, worksF or, Oxf ordU niversity) that T im is a person and Oxf ordU niversity is an organization.

KGs are highly flexible and can be used to represent general knowledge, such as that in Wikidata [van Veen, 2019], or more domain-specific data, such as Hetionet [START_REF] Scott Himmelstein | Systematic integration of biomedical knowledge prioritizes drugs for repurposing[END_REF] which stores biomedical data or FinKG [START_REF] Cheng | Knowledge graph-based event embedding framework for financial quantitative investments[END_REF] which stores financial data. In 2013, Min et al. pointed out that KGs are mostly incomplete, with some of the entities in the graph lacking crucial information [START_REF] Min | Distant supervision for relation extraction with an incomplete knowledge base[END_REF]. For example, in the Freebase KG [START_REF] Bollacker | Freebase: A collaboratively created graph database for structuring human knowledge[END_REF], 93.8% of the P erson nodes are not linked to a node of type P laceOf Birth and 78.5% of them are not linked to a node of type N ationality. From this observation, two main tasks emerged to discover the missing knowledge in a 1 Code to reproduce the experiments is available at https://gitub.u-bordeaux.fr/mdrance/ pre-trained-embeddings-for-enhancing-multi-hop-reasoning KG: (i) KG completion, which aims to retrieve the missing actual facts in the KG, and (ii) link prediction, which aims to discover the unknown links between the entities. These tasks have been used for specific uses cases, such as drug repurposing [START_REF] Drancé | Neuro-symbolic XAI for computational drug repurposing[END_REF][START_REF] Edwards | Explainable biomedical recommendations via reinforcement learning reasoning on knowledge graphs[END_REF] or in recommendation systems for social networks [START_REF] Wang | Link prediction in social networks: the state-of-the-art[END_REF].

Knowledge Graph Embedding

Knowledge graph embedding methods aim to construct latent vector representations for each entity and relationship in the KG. Given these embeddings, each KGE model defines its own scoring function f to evaluate the correctness of a given triple (h, r, t). KGE models are trained to distinguish between true and false triples. Given that, every triple belonging to the KG is considered true and false triples are constructed by corrupting either the head or the tail of an existing triple (h, r, t) → (h ′ , r, t ′ ). For each triple in the graph, the model learns to give a better score to the true triple than to its corresponding negatives. KGE models can be divided into two different categories [Ali et al., 2021a], defined by how the model uses embeddings to score each triple.

Translational Distance Interaction Models. These models use the notion of distance between embeddings to compute the plausibility of a triple. For example, TransE [START_REF] Bordes | Translating embeddings for modeling multi-relational data[END_REF] uses the relation embedding as a translation from head to tail embeddings. This translation e h + e r ≈ e t is defined by the scoring function as f (h, r, t) = -∥ e h + e r -e t ∥ p , with p ∈ {1, 2} the l p norm applied to the scoring function. The simplicity of TransE makes it an efficient KGE model to work with on large KGs, but its scoring function makes it impossible to model 1:N, N:1 and N:N relations. Following TransE, more robust models were designed, such as RotatE [START_REF] Sun | RotatE: Knowledge graph embedding by relational rotation in complex space[END_REF], aiming at addressing the limitations of TransE. They allow in particular to model symmetry, anti-symmetry, inversion and composition using embeddings laying in complex space C, modeling relations as rotations from the head to the tail: e t = e h ⊙ e r with e h , e r , e t ∈ C d and ⊙ being the Hadamard (or element-wise) product. Thus, the scoring function is defined as f (h, r, t) = -∥ e h ⊙ e r -e t ∥.

Semantic Similarity Matching Models. These models exploit the similarity of the latent features to compute the plausibility of a triple. RESCAL [START_REF] Nickel | A three-way model for collective learning on multi-relational data[END_REF] and Dist-Mult [START_REF] Yang | Embedding entities and relations for learning and inference in knowledge bases[END_REF] use vectors to model entities and matrices to model relations. The goal of the relation matrices W r ∈ R d×d is to learn the weights w i,j that quantify the interaction between head entities h ∈ R d and tail entities t ∈ R d . Hence, the RESCAL scoring function is given by f

(h, r, t) = h T W r t = d i=1 d j=1 w (r) ij h i t j .
The sole difference between RESCAL and DistMult is that the latter reduces the scoring function by utilizing only diagonal relation matrices:

f (h, r, t) = h T W r t = d i=1 h i • diag(W r ) i • t i .
As with RotatE, ComplEx [START_REF] Trouillon | Complex embeddings for simple link prediction[END_REF] improves DistMult by learning representations for entities and relations in C. The scoring function is defined as f (h, r, t) = Re(e h ⊙ e r ⊙ e t ) where Re(e t ) is the real component of the complex valued vectors representing t. Finally, ConvE [START_REF] Dettmers | Convolutional 2d knowledge graph embeddings[END_REF] makes use of a convolutional neural network layer to learn interactions between h and r. ConvE first takes the head and relation of each triple in the batch and concatenates them to create a matrix B ∈ R m×n where the first half m/2 rows represent all the h and the second half m/2 rows represent their corresponding r. Convolutional filters are then applied on B to capture interactions between h and r. These interactions are reshaped in order to obtain a feature vector v, that is then mapped to the entity space using a linear transformation W , finally creating conjoint representations of each pair of (h, r): e h,r = v T W . This representation of the first half of the triple is then scored against all the potential tails using the scoring function f (h, r, t) = e h,r • e t .

To the best of our knowledge, KGE models are the most efficient methods for the link prediction task. State-of-the-art results are obtained by using these models or improving existing KGE models [START_REF] Lu | DensE: An enhanced non-commutative representation for knowledge graph embedding with adaptive semantic hierarchy[END_REF][START_REF] Zhang | Quaternion knowledge graph embeddings[END_REF][START_REF] Chen | Relation prediction as an auxiliary training objective for improving multi-relational graph representations[END_REF][START_REF] Pan | Hyperbolic hierarchy-aware knowledge graph embedding for link prediction[END_REF]. The main issue is that embedding-based methods are not interpretable, although they can handle uncertainty and noise in the data. This is because each embedding represents a combination of latent factors, encoding the meaning of the entity/relation [START_REF] Bianchi | Knowledge graph embeddings and explainable AI[END_REF].

Related Work

MHR models are neuro-symbolic approaches in that they combine learning the vector representation of entities and relations with reasoning properties and interpretabilty provided to the user via the reasoning paths exploited for each prediction. MHR aims to create a new direct link between two entities using reasoning paths supporting the new direct link. To do so, a policy network (typically a fully-connected neural network) is trained to find the next best action for the agent given its current state. The sequence of successive actions then forms the reasoning path explaining the prediction. DeepPath [START_REF] Xiong | DeepPath: A reinforcement learning method for knowledge graph reasoning[END_REF] is the first attempt to model the path finding problem in a KG using the Monte-Carlo Policy Gradient algorithm, REINFORCE [START_REF] Williams | Simple statistical gradient-following algorithms for connectionist reinforcement learning[END_REF]. The particularity of DeepPath is that the source and target entities must be known, i.e. DeepPath learns to find correct reasoning paths between the two entities (e source , e target ), but not to predict a new link between unconnected nodes in the graph. MINERVA [START_REF] Das | Go for a walk and arrive at the answer: Reasoning over paths in knowledge bases using reinforcement learning[END_REF] is the first MHR model to actually address the link prediction problem using reinforcement learning. Indeed, instead of focusing on finding true paths between two entities, MINERVA trains its policy to reach a correct answer node given a query (e source , r query , ?), by traversing the best sequence of relations and entities supporting the choice of the predicted node. PoLo [START_REF] Liu | Neural multi-hop reasoning with logical rules on biomedical knowledge graphs[END_REF] is the first approach that attempts to decrease the effect of a noisy reward signal received by the agent when spurious As REINFORCE is a policy-based learning algorithm, the agent may be inclined to use spurious but rewarding paths encountered earlier in the training phase. Action dropout randomly disconnects some outgoing edges at each step, forcing the exploration of more diverse paths. RuleGuider [START_REF] Lei | Learning collaborative agents with rule guidance for knowledge graph reasoning[END_REF] proposes another strategy to improve MultiHopKG's reward shaping mechanism using high-quality logical rules, also dissociating the walk-based agent between a relation agent and an entity agent. As with PoLo, a symbolic method first mines the logical rules, which are then used to modify the agent's reward based on the confidence of the type of logical rule used to make the prediction. In addition, the relation agent will now select the type of outgoing relation to choose for the next step, then the entity agent will select the best node to move to, based on the starting node and the relation selected by the relation agent. This method significantly prunes the search space when selecting the next action.

Approach of Exploiting Pre-trained Embeddings

The current work proposes to analyze how pre-trained embeddings for entities and relations can be useful to enhance the performance of the link prediction task with MHR models. We tested this new approach on the state of the art model MultiHopKG, introduced in subsection 2.3, which we call hereafter PT-MultiHopKG. We describe here how MHR methods, and in particular MultiHopKG, model the problem of link prediction on a KG, followed by the details of our contribution based on the use of pre-trained embeddings.

Problem Definition

Environment and States. The KG G represents the environment, such that G ⊂ E × R × E where E represents the set of all entities in G and R represents the set of all relations in G. The link prediction task consists in finding the set of all possible answers e o ∈ E o given a query (e s , r q , ?), where e s is the source node and r q is the query relation, such that each (e s , r q , e o ) is a missing triple in G. These answer nodes are selected after the agent moves successively from one node to another until it reaches a plausible target. The state s t of the agent at step t should encode the source entity, the query relation and the location of the agent e t at step t. Thus, the current state s t is defined by s t = (e s , r q , e t ) ∈ S.

Action Space. The action space at step t comprises all outgoing edges from the current location e t : A t = {(r t ′ , e t ′ )|e s , r q , e t } with r t ′ ∈ R and e t ′ ∈ E. At each step, the action selection is made by selecting an outgoing edge, knowing the type of r t ′ and the next node e t ′ .

For each starting node e s , the search is limited to a certain number of steps T. To allow the agent to stay at its current position if it reaches a plausible answer at step t < T , an additional action "NO OP" is added, which corresponds to a self-loop on the current node e t .

Rewards. The base reward function defines a reward of 1 if the agent reaches a correct target entity and 0 otherwise: R = 1{(e s , r q , e o ) ∈ G}. For MultiHopKG, the reward is shaped using the scoring function f (e s , r q , e o ) of a KGE model. The idea is to give a reward of 1 if the triple (e s , r q , e o ) ∈ G, otherwise the reward is only defined by the KGE scoring function: R ′ = R+(1-R)f (e s , r q , e o ).

Policy Network. As defined in [START_REF] Das | Go for a walk and arrive at the answer: Reasoning over paths in knowledge bases using reinforcement learning[END_REF], the policy network makes use of three pieces of information to choose the appropriate action: the agent's current position e t , the query relation r q and the history of all previous actions of the agent. Each entity and relation in G is respectively assigned an embedding e ∈ R d and r ∈ R d , all actions A t at step t are represented as a t = [r t ′ ; e t ′ ] where [; ] is the vector concatenation of the chosen relation and its corresponding destination node. The history h t ∈ R 2d , encoded using an LSTM [START_REF] Hochreiter | Long Short-Term Memory[END_REF], represents the sequence of past observations and actions performed up to step t and is defined as follows:

h 0 = LST M (0, [r 0 ; e s ]) (1)
h t = LST M (h t-1 , a t-1 ) (2)
where equation ( 1) is used at step t 0 to encode the "starting action" in the history with r 0 a special "starting relation", and equation ( 2) is used otherwise. Based on this history h t , the current state s t and the query relation r q , the policy network π is defined as a two-layers feed-forward network that outputs the probability distributions of all possible actions A t at step t:

π θ (a t |s t ) = σ(A t × W 2 ReLU (W 1 [h t ; e t ; r q ])) (3)
with σ being the softmax function. Equations ( 2) and ( 3) are repeated for each transition step, until the maximum number of steps T is reached. The learnable parameters are the LSTM parameters, the feed-forward network parameters W 1 , W 2 and the entity and relation embeddings e, r.

Training

To train the policy network and find the best parameters θ, the REINFORCE algorithm is used to maximize the expected reward E as follows:

J (θ) = E (es,rq,eo)∈G [E a1,a2,...,aT ∼πθ [R(S T |e s , r q )]] (4) using the following stochastic gradient:

∇ θ J (θ) ≈ ∇ θ t R(s T |e s , r q )logπ θ (a t |s t ) (5)
During training, MultiHopKG proposes to add an action dropout mechanism to randomly hide some outgoing edges in the sampling step of REINFORCE. The goal is to encourage the agent to search for more diverse paths, as exploration in equation ( 5) can be biased towards spurious paths leading to correct answers.

Pre-trained Embeddings

The purpose of MHR models can be divided into two distinct parts:

• The learning of the entity and relation embeddings e and r as latent vector representations, describing the information carried by each node and each relation in the KG. This task is also performed in the KGE models. • The learning of the LSTM parameters in equation ( 2) and the feed-forward parameters W 1 and W 2 in equation (3), given the information of the agent history h t , the current node e t and the query relation r q . This task corresponds specifically to the walk-based method. During the training process, node and relation representations e and r are updated using equations ( 4) and ( 5), which are dependent on the rewards obtained at steps T . This training procedure implies that, in some cases, the representations will be updated using false information. First, when a false path is used to reach a good target node, the nodes and relations embeddings are updated using a positive reward based on a path that is completely illogical. Secondly, for correct predictions not present in the training data due to the incompleteness of the KGs, the representations will be updated using a wrong reward signal caused by a false negative.

In this work, we propose to use pre-trained embeddings for Method / Dataset WN18RR NELL-995 FB15K-237 H@1 H@10 MRR H@1 H@10 MRR H@1 H@10 MRR MultiHopKG ( ConvE entities and relations in the KG, generated by the KGE model ConvE. Consequently, the e, r embeddings in equations ( 1), ( 3), ( 4) and ( 5) are substituted with the corresponding pre-trained representations e + , r + corresponding to the pre-trained learnable parameters, or e -, r -corresponding to the pre-trained non-learnable parameters. The goal of this approach is to mitigate the impact of false information on entities and relations embeddings during the training process.

Experiments

Experimental Setup

Datasets. To compare our approach to other MHR methods, we evaluated its performance on the following three benchmark datasets: WN18RR [START_REF] Dettmers | Convolutional 2d knowledge graph embeddings[END_REF], Nell-995 [START_REF] Das | Go for a walk and arrive at the answer: Reasoning over paths in knowledge bases using reinforcement learning[END_REF] and FB15K-237 [START_REF] Toutanova | Representing text for joint embedding of text and knowledge bases[END_REF]. Their characteristics are given in Table 1. For each existing triple (h, r, t) in G, the inverse triple (t, r -1 , h) is added to allow bidirectional movement of the agent. At each step, the maximum number of outgoing edges is limited to a number n, to avoid GPU memory overflow. The top-n neighbors are selected using their PageRank scores [START_REF] Page | The PageRank Citation Ranking: Bringing order to the web[END_REF].

Hyperparameters. As previously indicated, we used the pre-trained embeddings generated by ConvE as defined in [Lin et al., 2018]. We keep the size of all embeddings at 200. We performed a grid search on embedding dropout rate [0, 0.5], feed-forward layers [0, 0.5], action dropout rate [0.1, 0.9] and learning rate [0.001, 0.003], respectively. We conducted the experiment using pre-trained embeddings as learnable parameters e + or having them frozen e -.

KGE model. Experiments were conducted using the embeddings generated by the KGE model used for reward shaping, resulting in no computational overhead. Tests carried out with the KGE model implemented using the Pykeen librairy [Ali et al., 2021b] showed exactly the same results as with the KGE model used for reward shaping. The process of training and testing the KGE model is made on exactly the same data splits used for training and testing the MHR model.

Path Analysis

A key feature of MHR models is their ability to provide reasoning paths to explain each prediction. MHR methods are expected to be more reliable than KGE methods, as the reason for each triple prediction can be easily understood. However, it has been shown that these reasoning paths are often unreasonable, i.e. paths that do not make sense but lead to a correct answer, or incomplete [START_REF] Lv | Is multi-hop reasoning really explainable? Towards benchmarking reasoning interpretability[END_REF]. We chose to compare how the addition of pre-trained embeddings on MultiHopKG changes the reasoning paths supporting each prediction. We measured: 1) the number of unique paths used as explanations, 2) the diversity of these paths, and 3) the number of test triples for which the model found an answer.

Unique Paths. MHR models often give more than one reasoning path to support each prediction, which leads to many explanations that decrease the interpretability of the results. Given two different queries, the same reasoning paths are often found for both predictions because, even though the entities present in the paths are different, the relations used are the same. The goal here is to abstract the reasoning paths into logical rules using only the relations, and measure how many unique logical rules are used for the explanations. For each path p = (h, r, t) ← (h, r 1 , e 1 ) ∧ (e 1 , r 2 , e 2 ) ∧ (e 2 , r 3 , t), the corresponding logical rule is l = (r 1 ∧ r 2 ∧ r 3 ). We then calculate the number of unique logical rules used on the set of predictions, which represents the model's ability to generalize the paths encountered during training onto the entities in the test set.

Rule Diversity. The number of unique logical rules provides a measure of the variety of explanations, but does not take into account the variety of each reasoning path. them are highly redundant because they use the same type of relation for more than one step. Diversity in the number of relation types employed by each logical rule results in more informative explanations, as it corresponds to rules with a greater variety of node types and semantically distinct relations. We quantify diversity as follows:

d = l U nique(l) |l| (6)
where U nique represents the number of unique relations r ∈ R found in the logical rules l.

Path Recall. We use the path recall score as defined in [START_REF] Lv | Is multi-hop reasoning really explainable? Towards benchmarking reasoning interpretability[END_REF] to quantify the number of triples in the test set that can be retrieved by the model. A significant path recall indicates that the model accurately predict and explain a greater proportion of the test triples. The path recall is defined as follows:

P R = (h,r,t)∈T test Cnt(h, r, t) |T test | (7)
where Cnt(h, r, t) = 1 if the model finds at least one path from h to t, and 0 otherwise.

Results and Discussion

Model Comparison

Table 2 reports the performance comparison between Multi-HopKG and PT-MultiHopKG. For all three datasets, the addition of pre-trained embeddings increases the predictive performance of the model. For NELL-995 and WN18RR, the use of frozen embeddings provides the best results, indicating that the ConvE-constructed representations are most efficient than using embeddings as model parameters. For NELL-995, using pre-trained embeddings increases the performances on both setup. For FB15K-237, the best performance is obtained by fine-tuning the embeddings. FB15K-237 is a more complex KG than both WN18RR and NELL-995 due to the number of triples it contains. Overall, although the two other datasets have more unique nodes, FB15K-237 contains more facts and each entity is much more connected to its neighbors (median degree of 14), making it a more complex KG to work with.

Table 3 shows the results of path analysis produced for MultiHopKG with and without pre-trained embeddings. Except for WN18RR, the addition of pre-trained embeddings increases the number of unique paths that serve as explanations for the model. This suggests that the improvement in performance is not solely attributable to the model's ability to predict more facts using the same rules, but rather to the model's use of different rules. The pre-trained embeddings improve path diversity and recall across all datasets. First, this indicates that the relationships between the logical rules employed are generally more diverse, i.e. the model is able to reuse a larger number of different logical paths seen during training when predicting new links. Second, the path recall values indicate that the model is able to identify at least one path for a greater number of query triples in the test dataset. In the case of FB15K-237, for instance, our model was able to provide a prediction and a reasoning path for 409 additional test queries.

Ablation Study and Model Variations

Using ComplEx. We conducted the same tests than previously, with the exception that we replaced ConvE with Com-plEx to generate the pre-trained embeddings. We used the best configuration for each dataset: frozen pre-trained embeddings for WN18RR and NELL-995 and learnable pre-trained embeddings for FB15K-237. Our results indicate that using ConvE is the best option and yields the best results.

Reward Shaping. The purpose of using pre-trained embeddings is to have prior information about each node and relation in the form of its embedding. Reward shaping mechanisms are also a way to use prior or external knowledge to assist the MHR model during training. As pointed out in [Lin et al., 2018], reward shaping improves performance on FB15K-237 and NELL-995, but decreases it for WN18RR. We compared the effectiveness of the two approaches by excluding the reward-shaping mechanism, as their goals are identical. Table 4 shows the results obtained for MultiHopKG with and without reward shaping mechanism using or not using pretrained embeddings. First, for all three datasets, the use of pre-trained embeddings increases the performance of the base model, highlighting that these embeddings fulfil their role in providing prior knowledge during the MHR training process. Secondly, reward shaping has a stronger impact on FB15K-237 but a negative impact on WN18RR, where the use of trainable pre-trained embeddings still has a positive impact on the MRR.

A summary of all results and a comparison with other symbolic and MHR methods can be found in Table 5.

Method / Dataset

WN18RR

NELL-995 FB15K-237 H@1 H@10 MRR H@1 H@10 MRR H@1 H@10 MRR AnyBURL † 42.9

Discussion

Pre-training has been widely studied in the past few years as a way to improve the predictive performance of a model. From image classification, where model pre-training is used to enhance feature detection [START_REF] Russakovsky | ImageNet Large Scale Visual Recognition Challenge[END_REF], to natural language processing, where pre-trained word embeddings are fine-tuned to match a specific domain [START_REF] Lee | BioBERT: a pre-trained biomedical language representation model for biomedical text mining[END_REF], pre-training allows for prior knowledge to be used or for the model to be tailored to a specific task. We believe that the experimental evaluation results presented in subsection 5.1 suggest that pre-trained embeddings have a positive impact on MultiHopKG. As in other areas of machine learning, model parameters need to be fine-tuned when applied to more complex tasks. Consequently, while embeddings can be used as is for the WN18RR and NELL-995 datasets, they must be included in the model parameters for the FB15k-237 dataset, which is significantly more complex. Since the value of MHR models lies in their interpretability, the analysis of reasoning paths has shown the interest of using pre-trained embeddings regarding explanations. Possessing more query triples for which a prediction can be made is a significant characteristic, but having more interpretable reasoning routes is even more crucial. Each prediction can be supported by hundreds of alternative reasoning paths, however these paths are frequently redundant in the type of relation they employ, traversing via various nodes while employing the exact same type of relation. Diverse reasoning paths enable more sophisticated and exhaustive explanations. Explanations can be more complex and exhaustive when there are a variety of reasoning paths. Compared to reward shaping, subsection 5.2 showed that pre-trained embeddings consistently improve the model performance, while reward shaping hurt performance for the WN18RR dataset. For FB15K-237, reward shaping has a greater impact than pre-trained embeddings, confirming the complexity of the link prediction task on this dataset, as it benefits the most from the addition of prior knowledge, but the best results were obtained using both methods simultaneously. Finally, when testing with the KGE model Com-plEx to build the embeddings, we observed inferior results to those obtained with ConvE. It is important to note that, when it comes to link prediction, ConvE always produces better results than ComplEx, as shown in [Ali et al., 2021a].

Conclusion

The current study investigates the impact of pre-trained embeddings on the state-of-the-art model MultiHopKG for the link prediction task. We demonstrated that adding pre-trained embeddings improves the performance of such model on the three benchmark datasets WN18RR, NELL-995 and FB15K-237. In addition to improving the results of the link prediction task, the results demonstrate that the pre-trained embeddings enable the model to predict new links for a larger number of queries using a larger number of unique and more diverse paths. This suggests that the method is ideally suited for use cases in which explanations are required, such as biomedical KGs. An ablation study shows that pre-trained embeddings are a valid method for providing prior knowledge during the MHR model training process. Furthermore, these results suggest that KGE methods are appropriate for building these pretrained embeddings. In future work, we plan to investigate the use of other methods to construct these pre-trained embeddings and validate the lessons learned from this study.

Ethical Statement

This work does not rely on sensitive personal data that could raise any ethical concern. However, the authors are fully aware that a great attention should be paid to ensure positive ethical and societal consequences of machine learning-based technologies.

Table 1 :

 1 Characteristics of the KGs used in our experiments. Degree corresponds to the median degree of nodes in the KG.paths are used but lead to a correct prediction. Since there exist many paths to connect two nodes, some of them may not be valid to support the final prediction. However, in MIN-ERVA the reward is only binary {0, 1}, which allows the agent to use false paths for valid predictions, i.e. paths that are meaningless regarding the new predicted link. Based on the implementation of MINERVA and to help guide the agent and remove noise from the reward signal, PoLo uses a set of known and efficient logical rules as a reward shaping mechanism: if a prediction uses one of these correct rules, the reward is increased.MultiHopKG [Lin et al., 2018] takes it further and proposes two modeling advances for MHR methods. First, the authors tackle the problem of missing true facts in the KG by adopting a new reward shaping mechanism. Be-

	Dataset	#Nodes #Relations	#Facts	Degree
	WN18RR	40,945	11	86,835	2
	NELL-995	75,492	200	154,213	1
	FB15K-237 14,505	237	272,115	14

cause KGs are incomplete

[START_REF] Min | Distant supervision for relation extraction with an incomplete knowledge base[END_REF]

, the agent may arrive at a correct answer that is not present in the training data, and thus not receive a reward for that prediction. Instead of using a binary {0, 1} reward, a trained KGE model is used to estimate a soft reward for predictions not existing in the training KG. For each new triple (h, r, t) predicted by the MHR model, if the triple is not present in the training KG, the KGE model scoring function f (h, r, t) is used to evaluate the plausibility of the prediction and shape the reward accordingly. Second, they enhance the agent's exploration ability by adding an action dropout for each training step.

Table 2 :

 2 Performance comparison on MultiHopKG with and without pre-trained embeddings. PT-MultiHopKG + corresponds to trainable pre-trained embeddings, PT-MultiHopKG -corresponds to frozen pre-trained embeddings. Best results in bold.

		)	41.4	51.7	44.8	65.6	84.4	72.7	32.7	56.4	40.7
	MultiHopKG (ComplEx) 42.5	52.6	46.1	64.4	81.6	71.2	32.9	54.4	39.3
	PT-MultiHopKG +	43.0	52.4	46.0	67.3	84.6	74.0	32.7	58.1	41.2
		PT-MultiHopKG -	44.1	52.8	46.8	68.1	85.6	74.9	31.1	57.2	39.9
		Unique Paths			Diversity			Path Recall
	Datasets	MultiHopKG PT-MultiHopKG MultiHopKG PT-MultiHopKG MultiHopKG PT-MultiHopKG
	WN18RR	1918	1790		2.76		2.78		0.62		0.63
	NELL-995	29,091	29,396		2.94		2.96		0.66		0.67
	FB15K-237	82,426	82,992		2.94		2.95		0.74		0.76

Table 3 :

 3 Path analysis. The number of unique paths used as explanations, the diversity of these paths and the path recall on the test set are reported.

Table 4 :

 4 Performance comparison (MRR) of the reward shaping and the impact of pre-trained embeddings on the model. MultiHopKG -RS corresponds to the model without reward shaping, MultiHopKG +RS corresponds to the model with ConvE reward shaping, Ours PT-MultiHopKG + -RS and Ours PT-MultiHopKG --RS correspond respectively to trainable pre-trained embeddings and frozen pre-trained embeddings without reward shaping. The percentage of improvement for each method is provided in parentheses.

Table 5 :

 5 Summary of our results and comparison with other symbolic and MHR models. Best scores are in bold. † Results taken from[START_REF] Lei | Learning collaborative agents with rule guidance for knowledge graph reasoning[END_REF]

			53.7	-	44.0	57.0	-	26.9	52.0	-
	MINERVA †	41.3	51.3	44.8	66.3	83.1	72.5	21.7	45.6	29.3
	MultiHopKG (ConvE)	41.4	51.7	44.8	65.6	84.4	72.7	32.7	56.4	40.7
	MultiHopKG (ComplEx) 42.5	52.6	46.1	64.4	81.6	71.2	32.9	54.4	39.3
	RuleGuider (ConvE) †	42.2	53.6	46.0	66.0	85.1	73.1	31.6	57.4	40.8
	RuleGuider (ComplEx) †	44.3	55.5	48.0	66.4	85.9	73.6	31.3	56.4	39.5
	PT-MultiHopKG +	43.0	52.4	46.0	67.3	84.6	74.0	32.7	58.1	41.2
	PT-MultiHopKG -	44.1	52.8	46.8	68.1	85.6	74.9	31.1	57.2	39.9
	PT-MultiHopKG + -RS	44.8	54.4	48.1	66.8	83.5	73.3	28.5	51.6	36.2
	PT-MultiHopKG --RS	45.5	55.8	49.0	63.9	84.6	71.4	27.2	51.2	35.2


Marie-Sophie C Ogloblinsky
email: marie-sophie.ogloblinsky@inserm.fr

Ozvan Bocher

Chaker Aloui

Anne-Louise Leutenegger

Ozan Ozisik

Anaïs Baudot

Elisabeth Tournier-Lasserve

Helen Castillo-Madeen

Daniel Lewinsohn

Donald F Conrad

Emmanuelle Génin

Gaëlle Marenne
email: gaelle.marenne@inserm.fr

PSAP-genomic-regions: a method leveraging population data to prioritize coding and non-coding variants in whole genome sequencing for rare disease diagnosis

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Each rare disease affects, by definition, a small number of individuals. However, as a whole, rare diseases affect about 350 million people world-wide [START_REF] Sequeira | The economic and health impact of rare diseases: A meta-analysis[END_REF]. Approximately 80% of rare diseases have a genetic origin that mostly follows a Mendelian mode of inheritance [START_REF] Amberger | McKusick's Online Mendelian Inheritance in Man (OMIM®)[END_REF][START_REF] Amberger | OMIM.org: Online Mendelian Inheritance in Man (OMIM®), an online catalog of human genes and genetic disorders[END_REF][START_REF] Ehrhart | A resource to explore the discovery of rare diseases and their causative genes[END_REF]. The advent of Next Generation Sequencing (NGS) and the development of variant pathogenicity prediction tools have allowed, in recent years, the identification of many genes involved in rare Mendelian diseases.

Nonetheless, despite extensive efforts, the molecular diagnosis is still unknown for more than 50% of rare diseases cases [START_REF] Wright | Paediatric genomics: diagnosing rare disease in children[END_REF][START_REF] Boycott | International Cooperation to Enable the Diagnosis of All Rare Genetic Diseases[END_REF][START_REF] Chong | The Genetic Basis of Mendelian Phenotypes: Discoveries, Challenges, and Opportunities[END_REF]. This can mainly be explained by the fact that many rare diseases are characterized by an extreme genetic heterogeneity, which results in only one individual carrying a specific pathogenic causal variant. This issue is referred to as the "n-of-one" problem [START_REF] Wilfert | Genomewide significance testing of variation from single case exomes[END_REF].

With the advent of high throughput sequencing technologies in clinics, molecular diagnosis is now often sought through whole exome or whole genome sequencing (WES and WGS respectively).

However, due to the large number of rare variants in each individual genome, causal variants are sought among very rare and highly pathogenic variants in genes relevant to the current known disease mechanism. The limited knowledge about gene functions and disease mechanisms can make this strategy unfruitful. To address the issue of variant prioritization at the level of an individual, the Population Sampling Method (PSAP) [START_REF] Wilfert | Genomewide significance testing of variation from single case exomes[END_REF] was developed. PSAP computes, for each gene, a null distribution, which is the probability to observe in the general population a genotype with a CADD pathogenicity score [START_REF] Kircher | A general framework for estimating the relative pathogenicity of human genetic variants[END_REF] greater than or equal to the highest one to the highest one observed in the patient for this gene. This initial version of the PSAP method, which we will refer to as PSAP-genes, has been successfully applied to identify variants of interest in diverse phenotypes, including male infertility [START_REF] Wyrwoll | Bi-allelic Mutations in M1AP Are a Frequent Cause of Meiotic Arrest and Severely Impaired Spermatogenesis Leading to Male Infertility[END_REF][START_REF] Kasak | Bi-allelic Recessive Lossof-Function Variants in FANCM Cause Non-obstructive Azoospermia[END_REF][START_REF] Salas-Huetos | Disruption of human meiotic telomere complex genes TERB1, TERB2 and MAJIN in men with non-obstructive azoospermia[END_REF], recurrent pregnancy loss [START_REF] Kasak | Recurrent Pregnancy Loss and Concealed Long-QT Syndrome[END_REF] and ciliary diskynesia [START_REF] Bustamante-Marin | Mutation of CFAP57, a protein required for the asymmetric targeting of a subset of inner dynein arms in Chlamydomonas, causes primary ciliary dyskinesia[END_REF].

A current hindrance to the application and generalization of PSAP-genes as a tool for diagnosis is its restriction to the coding parts of the genome. Indeed, the majority of variants reside in non-coding parts of the genome [START_REF] Hindorff | Potential etiologic and functional implications of genome-wide association loci for human diseases and traits[END_REF]. Non-coding variants may contribute to explain part of the etiology of rare diseases [START_REF] Posey | Genome sequencing and implications for rare disorders[END_REF], as suggested by the large number of GWAS hits located in non-coding regions of the genome [START_REF] Buniello | The NHGRI-EBI GWAS Catalog of published genome-wide association studies, targeted arrays and summary statistics 2019[END_REF]. The involvement of non-coding pathogenic variants in rare diseases is further corroborated by the fact that non-coding regions are heavily involved in the regulation of gene expression. Several prediction tools have been developed to this end [START_REF] Gussow | Orion: Detecting regions of the human non-coding genome that are intolerant to variation using population genetics[END_REF][START_REF] Huang | Fast, scalable prediction of deleterious noncoding variants from functional and population genomic data[END_REF][START_REF] Vitsios | Prioritizing non-coding regions based on human genomic constraint and sequence context with deep learning[END_REF], but most of them lack a variant-based score for both coding and non-coding regions. In addition, to be performant, they often require multiple annotations like Human Phenotype Ontology (HPO) terms [START_REF] Robinson | The Human Phenotype Ontology: A Tool for Annotating and Analyzing Human Hereditary Disease[END_REF] to characterize the symptoms or disease of a patient . Thus, they rely on previous knowledge and rarely go beyond candidate genes.

To move beyond the gene as a natural unit of testing for the PSAP method, we need to use predetermined regions across the whole genome. These regions also need to be defined using functional information to be used as a cohesive unit for the construction of PSAP null distributions. This challenge of defining regions along the whole genome has been tackled by Bocher et al. in the context of rare-variant association testing [START_REF] Bocher | Testing for association with rare variants in the coding and non-coding genome: RAVA-FIRST, a new approach based on CADD deleteriousness score[END_REF]: they describe CADD regions, which are characterized by a lack of observed variants with high functionally-Adjusted CADD Scores (ACS) in the gnomAD database [START_REF] Karczewski | The mutational constraint spectrum quantified from variation in 141,456 humans[END_REF].

CADD regions are expected to reflect functional constraints. CADD regions present the key advantage of providing pre-defined and functionally-informed regions which can be used to construct PSAP null distributions.

We have made available a new implementation of the PSAP method using Snakemake [START_REF] Köster | Snakemake-a scalable bioinformatics workflow engine[END_REF] workflows, called Easy-PSAP (https://github.com/msogloblinsky/Easy-PSAP), which features null distributions constructed with up-to-date allele frequency data and pathogenicity scores. Here, we introduce PSAP-genomic-regions, an extension of the PSAP method to the non-coding genome by using the pre-defined CADD regions as testing unit instead of genes. This is an innovative strategy to prioritize variants at the scale of an individual genome. PSAP-genomic-regions is now available in Easy-PSAP. We devised an evaluation protocol using artificially-generated disease exomes and genomes, obtained by inserting coding and non-coding ClinVar [START_REF] Landrum | ClinVar: improving access to variant interpretations and supporting evidence[END_REF] variants in general population whole genomes from the 1000 Genomes Project [START_REF] Auton | A global reference for human genetic variation[END_REF] and exomes from the FrEnch EXome (FREX) project [START_REF] Génin | The French Exome (FREX) Project: A Population-based panel of exomes to help filter out common local variants[END_REF]. We show the consistent improvement in prioritization by using PSAP-genomic-regions over pathogenicity scores alone for non-coding and then coding variants. For coding variants, we also demonstrate the good performance of PSAP-genomic-regions compared to PSAP-genes. On real-life data, we illustrate the power of PSAP-genomic-regions on WES data from six resolved cases of Cerebral Small Vessel Disease (CSVD) and WGS data from three families affected by male infertility. These two diseases are particularly relevant to test our method, monogenic forms of CSVD [START_REF] Rannikmäe | Beyond the Brain[END_REF] and male infertility (29) being extremely heterogeneous.

Results

Construction of PSAP null distribution in coding and non-coding regions

The idea behind the original PSAP method, referred to as PSAP-genes, relies on the calculation of gene-specific null distributions of CADD pathogenicity scores. More precisely, for an individual exome or genome and in a given gene, PSAP-genes considers the genotype with the highest CADD score and evaluates the probability to observe such a high CADD score in this gene in the general population (see S1 File for a detailed explanation of the calculation of PSAP null distributions). PSAP-genes deals separately with heterozygote and homozygote variants in the autosomal dominant (AD) and the autosomal recessive (AR) models respectively. As a result, PSAP-genes gives a p-value to the genotype with the highest CADD score in the gene for each gene, model, and individual. This p-value allows the ranking of the genes for an individual exome or genome. The PSAP principle can be generalized to any genomic unit.

Here, with PSAP-genomic-regions, we extended the PSAP method to analyze whole-genome data using predefined CADD regions as testing units instead of genes (Fig 1). The same principle as before is employed, with the difference being that the genotype with the highest CADD score in the region can be coding or non-coding. We thus constructed PSAP-genomic-regions null distributions with two pathogenicity scores : the initial CADD score (PHRED scaled across the whole genome), or the ACS [START_REF] Bocher | Testing for association with rare variants in the coding and non-coding genome: RAVA-FIRST, a new approach based on CADD deleteriousness score[END_REF] (PHRED scaled CADD scores by "coding", "regulatory" and "intergenic" regions) to mitigate the higher CADD scores of coding variants. Our two novel strategies will be referred to as PSAP-genomic-regions-CADD and PSAP-genomic-regions-ACS. They were compared to the initial PSAP-genes strategy, also referred to as PSAP-genes-CADD. We calculated PSAP null distributions for SNVs in genes and CADD regions, in the hg19 and hg38 assemblies of the human genome. In hg19, PSAP null distributions were obtained for 19,283 genes and 119,695 CADD regions. In hg38 PSAP null distributions were obtained for 18,395 genes and 123,991 CADD regions. PSAP null distributions and their parameters (unit of testing, allele frequencies and pathogenicity score) can be found in S1 Table. Using the ACS scores further improved the performance to detect non-coding-variants: 56.6% and 83.3% of variants reached the top 10 with PSAP-genomic-regions-ACS for the AD and AR models, respectively. Nonetheless, we can note the pattern is different for the top 1 for the AR model: 51% with PSAP-genomic-regions-CADD to 5.5% with PSAP-genomic-regions-ACS. Indeed, switching from CADD score to ACS score has lowered the PSAP p-value of non-coding variants shared by more than 10% of NFE genomes. This led to a defect of the top rank reached by the ClinVar variants, as we considered the lowest rank reached in at least 90% of individuals. For instance, a variant in the CADD region R109138 shared by 70 of the NFE genomes went from a CADD score of 18.1 and a PSAP-genomic-regions-CADD pvalue of 0.1 to an ACS of 22.2 and a PSAP-genomic-regions-ACS p-value of 5.18x10 -10 . Thus, the ClinVar variants inserted in these individuals having a higher p-value than 5.18x10 -10 do not rank first.

Evaluating the performance of PSAP-genomic-regions on artificiallygenerated disease exomes and genomes using ClinVar variants

Prioritization of non-coding pathogenic variants

We further explored PSAP results for splicing ClinVar variants versus other type of non-coding ClinVar variants. Indeed, we observed that splicing variants are the major type of non-coding ClinVar The full results of ranking by PSAP-genomic-regions-ACS for the non-coding non-splicing pathogenic ClinVar variants can be found in S3 File. With PSAP-genomic-regions-ACS, around half of the non-coding non-splicing variants are ranked in the top 100 of variants for more than 90% of NFE genomes (46 out of 73 variants for the AD model and 19 out of 31 variants for the AR model). The other half of variants present a less significant PSAP-genomic-regions-ACS p-value and a poorer ranking. To confirm this pattern of ranking for non-coding non-splicing pathogenic variants on another set of variants, we evaluated with our artificially generated disease genomes protocol 320 non-coding SNVs used to train Genomiser [START_REF] Smedley | A Whole-Genome Analysis Framework for Effective Identification of Pathogenic Regulatory Variants in Mendelian Disease[END_REF]. These variants were not associated with a mode of inheritance. Hence, we inserted them in the NFE genomes and scored them with both AD and AR PSAP-genomic-regions-ACS null distributions. Among the 320 non-coding variants, 169 reached the top 100 in at least 90% of NFE genomes, with either the AD or AR model (S4 File). This can be explained by the distributions of CADD scores compared to ACS scores for the ClinVar variants: the non-coding variants that do not reach the top 100 have a significantly lower CADD and ACS scores compared to all the other types of variants (S4 PSAP-genomic-region is also relevant for the analysis of exome data. Indeed, exome sequencing captures variants outside of the bounds of coding regions [START_REF] Guo | Exome sequencing generates high quality data in nontarget regions[END_REF], such as intronic variants. We explored the prioritization of non-coding ClinVar variants located within the WES-targeted regions of the FREX individuals using our artificially-generated disease exomes protocol (N=48 variants for the AD model and N=64 variants for the AR model, Panel A in S5 Fig) . For both PSAP-genomic-regions-CADD and PSAPgenomic-regions-ACS, there was a large increase in prioritization performance compared to using only the pathogenicity scores. Because there are fewer variants in an exome background than in a genome background, the rankings of these non-coding ClinVar variants were better in FREX than in NFE genomes. The best ranking was achieved using PSAP-genomic-regions-ACS, with 82% and 90.3% of variants reaching the top 10 for the AD and AR models, respectively. Most of these non-coding pathogenic variants were splicing variants (40 out of 73 variants for the AD model and 56 out of 64 variants for the AR model), and half of them were considered as having a functional "HIGH IMPACT" ( 26variants for the AD model and 22 variants for the AR model). Hence, prioritizing variants with PSAP on CADD regions allows identifying more variants even in exome data, that are in addition functionallyrelevant. We also compared the number of coding ClinVar variants reaching the tops in NFE genomes between PSAP-genomic-regions-CADD strategy and the initial PSAP-genes-CADD strategy (Fig 3). More differences were observed across the two PSAP strategies for the AD than for the AR model (Fig 3A).

Prioritization of coding pathogenic variants

There were 362 variants ranked first and 1,017 variants ranked [START_REF] Amberger | McKusick's Online Mendelian Inheritance in Man (OMIM®)[END_REF][START_REF] Amberger | OMIM.org: Online Mendelian Inheritance in Man (OMIM®), an online catalog of human genes and genetic disorders[END_REF][START_REF] Ehrhart | A resource to explore the discovery of rare diseases and their causative genes[END_REF][START_REF] Wright | Paediatric genomics: diagnosing rare disease in children[END_REF][START_REF] Boycott | International Cooperation to Enable the Diagnosis of All Rare Genetic Diseases[END_REF][START_REF] Chong | The Genetic Basis of Mendelian Phenotypes: Discoveries, Challenges, and Opportunities[END_REF][START_REF] Wilfert | Genomewide significance testing of variation from single case exomes[END_REF][START_REF] Kircher | A general framework for estimating the relative pathogenicity of human genetic variants[END_REF][START_REF] Wyrwoll | Bi-allelic Mutations in M1AP Are a Frequent Cause of Meiotic Arrest and Severely Impaired Spermatogenesis Leading to Male Infertility[END_REF] in common between the two strategies. However, 908 variants that were ranked [START_REF] Amberger | McKusick's Online Mendelian Inheritance in Man (OMIM®)[END_REF][START_REF] Amberger | OMIM.org: Online Mendelian Inheritance in Man (OMIM®), an online catalog of human genes and genetic disorders[END_REF][START_REF] Ehrhart | A resource to explore the discovery of rare diseases and their causative genes[END_REF][START_REF] Wright | Paediatric genomics: diagnosing rare disease in children[END_REF][START_REF] Boycott | International Cooperation to Enable the Diagnosis of All Rare Genetic Diseases[END_REF][START_REF] Chong | The Genetic Basis of Mendelian Phenotypes: Discoveries, Challenges, and Opportunities[END_REF][START_REF] Wilfert | Genomewide significance testing of variation from single case exomes[END_REF][START_REF] Kircher | A general framework for estimating the relative pathogenicity of human genetic variants[END_REF][START_REF] Wyrwoll | Bi-allelic Mutations in M1AP Are a Frequent Cause of Meiotic Arrest and Severely Impaired Spermatogenesis Leading to Male Infertility[END_REF] with PSAP-genes-CADD were with PSAP-genomic-regions-CADD, and 395 variants that were ranked [2-10] with PSAP-genes-CADD were ranked first with PSAP-genomic-regions-CADD. Regarding variants that are ranked more than a 100 with PSAP-genomic-regions-CADD, 278 of them are ranked

Application of PSAP-genomic-regions to real data with different modes of inheritance

To illustrate our method in real-life settings, we analyzed two datasets (S4 Table ), one with an AD mode of inheritance and the other with an AR mode of inheritance. The first dataset consisted of WES data for six individuals affected by monogenic forms of CSVD [START_REF] Aloui | End-Truncated LAMB1 Causes a Hippocampal Memory Defect and a Leukoencephalopathy[END_REF]. Using PSAP-genomic-regions-CADD, all of the causal variants were ranked at least in the top 100 in each patient (Fig 4). The contribution of CADD regions as a unit of testing was especially visible for the variant in COL4A2 and one variant in HTRA1 which were not well-ranked using genes as testing unit (rank 110 and 193 respectively with genes, and rank 3 and 69 with CADD regions). Using their maximal CADD score by gene or CADD region alone, these variants would not have been prioritized in the top 100 for five out of six individuals.

Fig. 4. Prioritization of 6 known CSVD mutations and 3 male infertility candidate variants with PSAPgenomic-regions-CADD, PSAP-genes-CADD and the maximal CADD score on genes or CADD regions.

The second dataset consisted of WGS data for 9 individuals from three families with clinically diagnosed male infertility [START_REF] Khan | Genome sequencing of Pakistani families with male infertility identifies deleterious genotypes in SPAG6, CCDC9, TKTL1, TUBA3C, and M1AP[END_REF]. All causal variants fell within the top 20 of variants with prioritization by PSAP-genes-CADD, and within the top 50 for at least one case per family with PSAP-genomic-regions-CADD (within top 100 for all cases, Fig 4). PSAP-genomic-regions-CADD did not improve the ranking of these coding variants, which was expected considering the large number of variants in a WGS analysis (see S4 Table for the total number of variants in each analysis). The prioritization from PSAP-genomic-regions-CADD was still interesting to narrow the set of candidates for causal variants. In clinics when the CADD score alone is used, these variants would not have been prioritized (CADD score < 25, and rank > 100 with the maximal CADD score strategy). PSAP-genomic-regions-CADD thus allow a relevant prioritization of coding pathogenic variants in WGS sequencing and an unbiased exploratory analysis at the scale of the whole genome.

Using PSAP-genomic-regions-ACS or the ACS score alone, almost all of the CSVD and male infertility coding pathogenic variants had a rank greatly exceeding the top 100 (S4 Table ). The only exception is one variant in HTRA1 (10:124266885 G/A) that was ranked 3 by PSAP-genomic-regions-ACS and 10 by the maximal ACS score alone. This HTRA1 variant was a splicing variant, which confirms the good performance of the PSAP-genomic-regions-ACS strategy on this type of variant.

Discussion

Variant prioritization, especially in the case of very heterogeneous rare diseases, is a clinicallyrelevant methodological challenge for both clinicians and researchers. Mounting evidence suggests that current methods of analysis and their restriction to the coding genome are a hindrance to the discovery of new genetic variants implicated in rare diseases [START_REF] Posey | Genome sequencing and implications for rare disorders[END_REF]. We have developed PSAP-genomic-regions, an extension of the PSAP method to the whole genome using functionally-relevant genomic regions. PSAPgenomic-regions broadens the scope of variants evaluated by PSAP and addresses the issue of variant prioritization at an individual whole-genome scale.

PSAP-genomic-regions has been thoroughly tested and validated by using simulations emulating real-life scenarios of causal variant prioritization. PSAP-genomic-regions achieves a prioritization of coding pathogenic SNVs in the top 100 variants of an exome or genome which is a relevant number of variants to analyze for clinicians. Without use of prior knowledge on the disease, PSAP-genomic-regions achieves relevant variant prioritization within millions of variants to analyze, which is illustrated by the ranking of 6 variants involved in CSVD and 3 variants involved in familial cases of male infertility in the top 100 of WES and WGS data respectively. PSAP-genomic-regions thus helps with the diagnosis of such heterogeneous diseases in conjunction with other relevant information like the mode of transmission, prevalence or type of variant involved. PSAP-genomic-regions also allows the scoring of variants otherwise discarded from the analysis, like splicing variants with a high predicted functional impact, and other non-coding variants of proven clinical significance. The only scenario for which PSAP-genomic-regions is not advantageous compared to the PSAP-genes strategy is for prioritizing coding variants in WGS data. In that case, using coding CADD regions, i.e. the coding parts of CADD regions for the analysis still yields better results compared to PSAP-genes (S7 Fig) . Our simulations using known pathogenic variants have shown which PSAP strategy performs the best depending on the type of data and variant expected to be involved in the disease mechanism (S8 Fig) . To effectively prioritize non-coding variants in WES and WGS, we advise the use of PSAP-genomic-regions-ACS. For coding variants, PSAP-genomic-regions-CADD gives the best results in WES, and PSAP-coding-genomic-regions-CADD performs best in WGS data. A two-step approach can also be carried out if there is no expected type of variant: first, the PSAP-genomic-regions-CADD or PSAP-coding-genomic-regions-CADD strategy is applied depending on the type of data, and if no coding variant of interest for the disease is found within the top results, PSAP-genomic-regions-ACS can be applied to look for non-coding variants of interest.

To the best of our knowledge, there is no other score of predicted pathogenicity for all possible SNVs comparable to CADD. Other methods have been developed to distinguish between coding pathogenic and neutral variants [START_REF] Niroula | PON-P2: Prediction Method for Fast and Reliable Identification of Harmful Variants[END_REF][START_REF] Alirezaie | ClinPred: Prediction Tool to Identify Disease-Relevant Nonsynonymous Single-Nucleotide Variants[END_REF][START_REF] Choi | PROVEAN web server: a tool to predict the functional effect of amino acid substitutions and indels[END_REF][START_REF] Ng | SIFT: predicting amino acid changes that affect protein function[END_REF][START_REF] Rogers | FATHMM-XF: accurate prediction of pathogenic point mutations via extended features[END_REF][START_REF] Li | CAPICE: a computational method for Consequence-Agnostic Pathogenicity Interpretation of Clinical Exome variations[END_REF], but often restrict to non-synonymous variants. These methods were shown to perform better or have advantages compared to CADD for the limited set of variants they explore [START_REF] Niroula | PON-P2: Prediction Method for Fast and Reliable Identification of Harmful Variants[END_REF][START_REF] Alirezaie | ClinPred: Prediction Tool to Identify Disease-Relevant Nonsynonymous Single-Nucleotide Variants[END_REF][START_REF] Choi | PROVEAN web server: a tool to predict the functional effect of amino acid substitutions and indels[END_REF][START_REF] Ng | SIFT: predicting amino acid changes that affect protein function[END_REF][START_REF] Rogers | FATHMM-XF: accurate prediction of pathogenic point mutations via extended features[END_REF][START_REF] Li | CAPICE: a computational method for Consequence-Agnostic Pathogenicity Interpretation of Clinical Exome variations[END_REF]. Similar types of methods aim at prioritizing more constrained regions in the noncoding genome [START_REF] Gussow | Orion: Detecting regions of the human non-coding genome that are intolerant to variation using population genetics[END_REF][START_REF] Vitsios | Prioritizing non-coding regions based on human genomic constraint and sequence context with deep learning[END_REF] or distinguishing deleterious non-coding variants from neutral ones [START_REF] Gussow | Orion: Detecting regions of the human non-coding genome that are intolerant to variation using population genetics[END_REF][START_REF] Caron | NCBoost classifies pathogenic non-coding variants in Mendelian diseases through supervised learning on purifying selection signals in humans[END_REF].

Other well-known methods for identification of pathogenic variants in exome and genome data rely on the use of HPO terms to make a prediction, like Exomiser [START_REF] Smedley | Next-generation diagnostics and disease-gene discovery with the Exomiser[END_REF] or Genomiser [START_REF] Smedley | A Whole-Genome Analysis Framework for Effective Identification of Pathogenic Regulatory Variants in Mendelian Disease[END_REF], making in comparison PSAP an unmatched prioritization tool. As any other bioinformatics variant prioritization method, it has to be used in conjunction with other lines of evidence to ultimately lead to any genetic diagnosis of a patient. PSAP-genomic-regions does not make assumption on the type of variants and does explore the whole genome. The ranking by p-values coming from the application of PSAP-genomic-regions to an individual's variants is a useful way to narrow-down the list of variants to further investigate for both researchers and clinicians in different scenarios.

The method most comparable to the strategy followed by PSAP-genomic-regions is the recentlydeveloped machine-learning algorithm FINSURF [START_REF] Moyon | Classification of non-coding variants with high pathogenic impact[END_REF]. FINSURF aims to predict the functional impact of non-coding variants in regulatory regions and has been applied to known pathogenic variants inserted in WGS data like we did. Nonetheless it has been difficult to compare properly the two methods considering FINSURF only scores non-coding variants in predefined regulatory regions, and the set of variants used to train the method is not available.

The main limitation of PSAP-genomic-regions comes from the score used to calibrate null distributions, namely the CADD score. We have observed that known pathogenic non-coding ClinVar variants that were not well-ranked by PSAP-genomic-regions had significantly lower CADD and ACS scores compared to splicing and better-ranked non-coding variants. Because such CADD score is likely to be seen in the general population, PSAP-genomic-regions will not be able to prioritize such a variant with at a low rank. We also observed that some CADD regions were badly-calibrated and resulted in the assignment of very low PSAP-genomic-regions p-values to putatively neutral variants in the 1000 Genomes Project. As allele frequencies from larger databases and more accurate pathogenicity scores become available, this will lead to an improvement of the PSAP method as well. The most recent release of the CADD score v1.7 [START_REF] Schubach | CADD v1.7: using protein language models, regulatory CNNs and other nucleotide-level scores to improve genome-wide variant predictions[END_REF] notably integrates regulatory annotations and may further improve the prioritization of non-coding pathogenic variants when integrated in PSAP-genomic-regions.

Many avenues of further development and improvement are open for PSAP-genomic-regions, including the inclusion and scoring of InDel variations and structural variants. Exploring the combination of the PSAP-genomic-regions p-values with other metrics or information coming from omics analysis could also improve prediction. Finally, the flexibility of the PSAP method makes it potentially adaptable to other more complex models like digenic and oligogenic models of inheritance, considering the increasing availability of information coming from gene networks and biological pathways. population cohorts. These different PSAP strategies (see Table 1) were compared in terms of their performances to prioritize the known pathogenic variants.

The pathogenic ClinVar (25) SNVs with coordinates in hg19 and hg38 were downloaded from the NCBI website (https://www.ncbi.nlm.nih.gov/clinvar/, accessed on the 3rd of June 2022). Some of these ClinVar variants had an annotated mode of inheritance ("moi autosomal recessive" and "moi autosomal dominant"). From ClinVar, there were 12,776 variants annotated as AD and 12,776 variants annotated as AR. Variants were filtered out to keep only autosomal pathogenic SNVs having as review status either "reviewed by expert panel" or "criteria provided, multiple submitters, no conflicts", which are the two best review status in ClinVar. There were 1,518 AD and 1,118 AR variants meeting these criteria.

For variants which did not have an annotated mode of inheritance, we used a curated version of the database OMIM, hOMIM [START_REF] Blekhman | Natural selection on genes that underlie human disease susceptibility[END_REF] to retrieve a mode of inheritance, and kept variants that were always associated with an AD or AR mode of inheritance in hOMIM. The same filtering was applied, which left 3,641 additional variants for the AD and 1,706 for the AR model. In total, we had a set of 5,159 variants for the AD model and 2,824 variants for the AR model. Among these ClinVar variants, 4,965 and 2,680 variants were coding SNVs respectively for the AD and AR models. Similarly, 175 and 96 variants were non-coding variants for the AD model and AR models, among which 48 variants for the AD model and 64 for AR model fell within the boundaries covered by FREX exomes. The list of pathogenic ClinVar variants and their mode of inheritance can be found in S2 File.

We inserted each variant from our curated list of pathogenic ClinVar variants successively in each of the 533 high coverage genomes of Non-Finnish Europeans (NFE) from the 1000 Genomes Project phase 4 (NFE genomes) and each of the 574 exomes from the FREX project. An individual-focused QC was applied on both datasets using the RAVAQ R package (50): we performed a genotype and variant QC with default parameters corresponding to standard GATK hard filtering criteria, mean allele balance The rank of the known CSVD variants among other heterozygote variants in the patient's exome according to its PSAP p-value for the 2 strategies was then retrieved.

The PSAP strategies were also applied to WGS data of three families with clinically diagnosed forms of male infertility [START_REF] Khan | Genome sequencing of Pakistani families with male infertility identifies deleterious genotypes in SPAG6, CCDC9, TKTL1, TUBA3C, and M1AP[END_REF] and for which a pathogenic recessive variant was prioritized using a computational pipeline featuring the initial PSAP-genes implementation. Three affected individuals were analyzed for each family. The description of the whole dataset and candidate variant filtering process can be found in [START_REF] Khan | Genome sequencing of Pakistani families with male infertility identifies deleterious genotypes in SPAG6, CCDC9, TKTL1, TUBA3C, and M1AP[END_REF]] (33), except for the QC that was performed in the same way as for the CSVD data. Two other families were resolved from the same dataset, but considering that the causal variants were deletions we did not include them in the current analysis. The prioritized pathogenic variants were in the genes: SPAG6 (chr10:22389235 C/T) for family 3, TUBA3C (chr13:19177247 C/T) for family 7 and CCDC9 (chr19:47260609 C/T) for family 4. We applied PSAPgenes-CADD and PSAP-genomic-regions-CADD in hg38 to the 9 cases and retrieved the rank of the known male infertility variants among other homozygote variants in the patient's genomes according to its PSAP p-value for the 2 strategies.

Fig 1 .

 1 Fig 1. Description of the PSAP-genomic-regions strategy.

First, to evaluate

  how PSAP-genomic-regions performed to prioritize non-coding pathogenic variants, we used artificially-generated disease genomes created by inserting non-coding ClinVar variants in the NFE genomes from 1000G project (see Material & Methods and S2 File for the list of variants). Because the 1000 Genomes project is population-based, we expect that some individuals might carry one or a few pathogenic variants in their genome. These pathogenic variants are characterized by a high CADD score and a low PSAP p-value. Indeed, there is large variation in the maximal CADD score or lowest PSAP p-value, whereas the rest of the distribution is extremely similar between individuals (S1 Fig). Thus, in order to summarize the rank of a ClinVar variant in an evaluation setting, we considered the best rank reached by the variant in at least 90% of the individuals. Most of the NFE genomes carried a variant with a higher pathogenicity score or a lower PSAP pvalue than most of the ClinVar variants (S2 Fig). We thus compared the percentage of the non-coding pathogenic variants ranked among the top N (N = 1, 10, 50 and 100) in at least 90% of the NFE genomes. The ranking at the individual level was done among all heterozygous variants for the ClinVar variants under the AD model, and across homozygous variants for the ClinVar variants under the AR model. (Fig 2A). With both CADD and ACS pathogenicity scores, PSAP-genomic-regions performed systematically better than using the pathogenicity scores alone. The improvement was especially large for the top 10 ranking: 24.6% and 79.2% of ClinVar variants reached the top 10 with PSAP-genomic-regions-CADD for the AD and AR models, respectively, while no ClinVar variant reached the top 10 with CADD scores alone.

Fig 2 .

 2 Fig 2. Comparison of the PSAP-genomic-regions strategy versus a pathogenicity score alone for in

  variants. These splicing variants often had a very good ranking, especially with PSAP-genomic-regions-ACS (n=115 splicing variants among 175 non-coding AD variants and n=72 splicing variants among 96 non-coding AR variants; S3 Table; Panel A in S3 Fig). Splicing ClinVar variants have a much higher ACS than CADD scores (Panel B in S3 Fig) which results in better ranking than for other types of non-coding ClinVar variants using PSAP-genomic-regions-ACS p-values (Panel C in S3 Fig). As a consequence, the percentage of splicing ClinVar variants ranked in the top 10 was largely improved when using PSAPgenomic-regions-ACS, for the AD model especially which was less powerful with PSAP-genomic-regions-CADD to begin with (Panel D in S3 Fig).

Fig

  Fig). Overall, PSAP-genomic-regions-ACS prioritizes around half of non-coding ClinVar and Genomiser

  Similar evaluations were performed for ClinVar coding variants inserted in either WGS from 1000G NFE individuals or WES from FREX. As observed for non-coding pathogenic variants, PSAPgenomic-regions outperformed the pathogenicity scores alone (Fig 2B, Panel B in S5 Fig). However, in the context of coding pathogenic ClinVar variants, we observed that the strategy of PSAP-genomicregions-CADD provided better prioritization compared with the PSAP-genomic-regions-ACS strategy. We observed that 18.2% and 74.6% of the coding variants reached the top 1 in at least 90% of genomes backgrounds with the PSAP-genomic-regions-CADD for the AD and AR model respectively, against no variants with the CADD score alone, and against 5.3% and 2.5% reaching the top 1 with PSAP-genomicregions-ACS. In the exome background and with PSAP-genomic-regions-CADD, 38.7% and 89.8% of AD variants reached the top 1 and top 50, respectively; 80.3% and 97.9% of AR variants reached the top 1 and the top 50, respectively.

Fig. 3 .

 3 Fig.3. Comparison of PSAP-genomic-regions-CADD and PSAP-genes-CADD strategies in artificially-

  and are ranked [51-100] by PSAP-genes-CADD. Regarding the AR model (Fig 3B), PSAP-genomic-regions-CADD performed similarly to PSAPgenes-CADD, and the majority of variants were ranked first with both strategies (1,550 variants). Even more promising results can be found when looking at the same comparison of ranks within the FREX exomes (S6 Fig). For instance, in the AD model, 592 variants that were ranked [2-10] with PSAP-genes-CADD are ranked first with PSAP-genomic-regions-CADD, against 115 variants ranked [2-10] with PSAPgenomic-regions-CADD that become first with PSAP-genes-CADD.

( 19 :

 19 15303053 G/A and 19:15303260 G/A), one individual in the gene COL4A2 (13:111132702 G/T) and three individuals in the gene HTRA1 (10:124266285 T/G, 10:124266281 C/A and 10:124266885 G/A).

Materials and Methods

Construction of PSAP null distributions

The first parameter is the units in which to construct the PSAP null distribution. Here we considered two unit strategies: the genes and the CADD regions (S1 Table ). For the genes, the coding regions of genes were defined based on the biomaRt R package: the gene coding sequences were retrieved from Ensembl (44) by requesting the "genomic_coding_start" and "genomic_coding_end", on both the hg19 and hg38 builds. To account for splicing regions, the coding regions were extended by two bases on both sides of the gene coding regions. In total, 19,780 genes were retrieved in hg19 and 23,163 in the hg38 build. For the CADD regions, their coordinates were downloaded from https://lysine.univ-brest.fr/RAVA-FIRST/ for the hg19 build and were lifted over to hg38 using the Ensembl Assembly Converter. CADD regions coordinates in hg38 are available on Easy-PSAP GitHub (https://github.com/msogloblinsky/Easy-PSAP). There were 135,224 CADD regions in hg19 and 131,970 in hg38. For the coding CADD regions, i.e. the coding parts of CADD regions, we considered the intersection of the CADD regions and the gene coding regions for each build, which yielded 37,978 coding CADD regions in hg19 and 52,340 in hg38.

The second parameter is the allele frequencies database. Here we considered the global allele frequencies from the gnomAD database to calibrate the PSAP null distributions: gnomAD genome r2.0.1 for hg19 and gnomAD V3 (45) for hg38. For our purpose, we considered only single nucleotide variants (SNVs) annotated as PASS by the Variant Quality Score Recalibration (VQSR) of GATK [START_REF] Mckenna | The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data[END_REF] and located in well-covered regions. Well-covered regions in gnomAD genome were defined as regions for which 90% of individuals have coverage at depth 10. Variants not seen in gnomAD genome, not annotated as PASS or not located in well-covered regions (gnomAD genome version according to the build) have a frequency of 0 and thus did not contribute to the construction of the null distributions.

To ensure reliability of PSAP null distribution, it is crucial that the units are well covered in the database from which the allele frequencies are taken. Thus, we only considered units for which at least half of the unit was well-covered (as defined previously) in gnomAD genome (version according to the build). Coding regions of genes and well-covered regions in gnomAD genome were intersected to get the percentage of each gene's coding regions that were well-covered in the database. The same steps were carried out with CADD regions as genomic units for PSAP, for hg19 and hg38 builds. PSAP null distributions were thus constructed for 19,283 and 18,395 genes in hg19 and hg38 respectively, 119,695 and 123,991 CADD regions, and 34,397 and 35,226 coding CADD regions in hg19 and hg38 respectively.

The third parameter is the pathogenicity score. Here, for the evaluation of PSAP on coding variants, we used the version 1.6 of CADD [START_REF] Rentzsch | CADD-Splice-improving genome-wide variant effect prediction using deep learning-derived splice scores[END_REF] for each build, accessible on the CADD website (https://cadd.gs.washington.edu/). For the evaluation on non-coding variants, which tend to have lower CADD scores than coding variants (48), we followed the strategy described in Bocher et al. [START_REF] Bocher | Testing for association with rare variants in the coding and non-coding genome: RAVA-FIRST, a new approach based on CADD deleteriousness score[END_REF] to adjust the RAW CADD score v1.6 of all possible SNVs on a PHRED scale stratifying by type of genomic regions: "coding", "regulatory" and "intergenic", resulting in "adjusted CADD scores", referred to as "ACS".

Easy-PSAP (https://github.com/msogloblinsky/Easy-PSAP) was used to generate null distributions according to the previously described input files and parameters. This resulted in 4 sets of null distributions for the AD and AR models for both hg19 and hg38 assemblies (S1 Table ).

Evaluating the performance of PSAP-genomic-regions using artificially-generated disease exomes and genomes

To evaluate the ability of PSAP-genomic-regions to prioritize known pathogenic variants in an individual, we leveraged artificially-generated disease exomes and genomes using available general computed across heterozygous genotypes and call rates, except for MAX_AB_GENO_DEV = 0.25, MAX_ABHET_DEV, MIN_CALLRATE and MIN_FISHER_CALLRATE "disabled".

We conducted the artificially-generated disease genome and exome evaluation with PSAP null distributions in hg19 and hg38 respectively, to match with the build of the data. We then applied the 3 PSAP strategies mentioned previously (PSAP-genes-CADD, PSAP-genomic-regions-CADD and PSAPgenomic-regions-ACS). For each strategy, we kept the maximal pathogenicity score (CADD or ACS) for each unit (gene or CADD regions) and then ranked the units according to their PSAP p-value or to their pathogenicity score alone within each genome or exome. We compared the PSAP-genes-CADD and PSAP-genomic-regions-CADD strategies to using the maximal CADD score alone by gene or CADD regions, respectively; and the PSAP-genomic-regions-ACS strategy to using the maximal ACS score by CADD region. For each ClinVar variant, we retrieved its rank within each genome or exome. Coding ClinVar variants were evaluated with the 3 PSAP strategies whereas non-coding ClinVar variants were evaluated with the novel PSAP-genomic-regions-CADD and PSAP-genomic-regions-ACS strategies (see S2

Patient data analysis

The PSAP strategies were applied to real WES data from six unrelated patients affected by a CSVD for which the causal variant is known, which allowed a comparison of performance between the different strategies. The full description of the dataset can be found in [START_REF] Aloui | End-Truncated LAMB1 Causes a Hippocampal Memory Defect and a Leukoencephalopathy[END_REF]] (32), with the exception of the QC process. For this analysis, the same QC as for the FREX and 1000 Genomes Project datasets was performed. We applied PSAP-genes-CADD and PSAP-genomic-regions-CADD in hg19 to the six resolved CSVD patients' exome data. The other PSAP parameters were the ones by default as

S1