Ania De La Nuez Veulens

Yoanna Ginarte

Rolando Fernandez

Fabrice Leclerc

Luis A Montero Cabrera

Yoanna M Álvarez Ginarte

Rolando E Rodríguez Fernandez

Montero Cabrera

Prediction of molecular interactions and physicochemical properties relevant for vasopressin V2 receptor antagonism

Keywords: G Protein Coupled Receptor, PKD: Polycystic Kidney Disease, ADPKD: Autosomal Dominant Polycystic Kidney Disease, V2R, AVP, Vasopressin Antagonist, GPCR, QSAR, Docking, Molecular Dynamics Simulations

INTRODUCTION

Autosomal Dominant Polycystic Kidney Disease (ADPKD) is a genetic condition with an incidence of 1:400, to 1:1000 in the world population. Patients develop multiple fluid-filled cyst in both kidneys, increasing the total kidney volume and leading to Chronic Kidney Disease. In contrast to the normal renal cells, the PKD cystic cells have an increased cAMP induced proliferation by activating the Ras/B-Raf/MEK/ERK pathway, mediated by the decrease in intracellular calcium levels [START_REF] Yamaguchi | Calcium restriction allows cAMP activation of the B-Raf/ERK pathway, switching cells to a cAMP-dependent growthstimulated phenotype[END_REF]. The absence of renal cysts in PCK rats lacking arginine vasopressin (AVP) might indicate that the receptors on the collecting ducts activating other Adenylate Cyclases, do not have a significant role on the generation of cysts [START_REF] Torres | Effective treatment of an orthologous model of autosomal dominant polycystic kidney disease[END_REF][START_REF] Grantham | Rationale for early treatment of polycystic kidney disease[END_REF].

AVP is a nonapeptide hormone consisting of a six amino acid ring closed by a disulphide bridge between cysteines 1 and 6, followed by a tripeptide tail. This hormone is synthesized in the hypothalamus, principally produced by neurons with the cell body within the supraoptic and the paraventricular nuclei, and their axon terminations in the neural lobe of the posterior pituitary gland in which AVP is released into the circulation [START_REF] Decaux | Non-peptide arginine-vasopressin antagonists: the vaptans[END_REF]. The primary function of AVP is to maintain body fluid balance by keeping plasma osmolality within narrow limits [START_REF] Bankir | Vasopressin: physiology, assessment and osmosensation[END_REF]. Increase in plasma osmolality or decrease in plasma volume trigger its release to induce expression of water transport proteins in the late distal tubule and collecting ducts of the kidneys, to increase water reabsorption [START_REF] Bankir | Vasopressin: physiology, assessment and osmosensation[END_REF]. Due to its role in the regulation of osmolarity by increasing the ability of the kidney to reabsorb water reducing the urinary volume, it is also known as Antidiuretic Hormone (ADH).

AVP physiological roles are mediated by three receptor subtypes V1a, V1b (also called V3), and V2 all belong to vasopressin/oxytocin receptor family and they are class-A G-protein coupled receptor (GPCR). The V1a receptors are mainly distributed on vascular smooth muscle, but also present in myocardium, platelets, and hepatocytes. V1a stimulation is associated with vasoconstriction and cardiac hypertrophy, together with platelet aggregation, and glycogenolysis [START_REF] Decaux | Non-peptide arginine-vasopressin antagonists: the vaptans[END_REF][START_REF] Thibonnier | Molecular pharmacology and modeling of vasopressin receptors[END_REF][START_REF] Facciorusso | Role of vaptans in the management of hydroelectrolytic imbalance in liver cirrhosis[END_REF]. The V1b receptors have little selective distribution and their activation is part of the adaptive reaction to stress, leading to stimulation of adrenocorticotropic hormone and endorphin release [START_REF] Decaux | Non-peptide arginine-vasopressin antagonists: the vaptans[END_REF][START_REF] Thibonnier | Molecular pharmacology and modeling of vasopressin receptors[END_REF]. The activity of each receptor is mediated by G proteins which activate a phosphatidyl-inositol-calcium second messenger system. The V2 receptor (V2R) is expressed predominantly in the principal cells of the renal collecting duct system, in which its activation leads to increased resorption of free water [START_REF] Decaux | Non-peptide arginine-vasopressin antagonists: the vaptans[END_REF][START_REF] Bankir | Vasopressin: physiology, assessment and osmosensation[END_REF][START_REF] Thibonnier | Molecular pharmacology and modeling of vasopressin receptors[END_REF]. V2R is the major activator of adenylyl cyclase signaling pathway in principals cell of collecting ducts in kidney. The increase of cAMP intracellular concentration by the activation of V2R promote proliferation in PKD cystic cells, suggesting that V2R antagonists can be used as treatment for PKD to retard development and growth of the cysts [START_REF] Aihara | Tolvaptan delays the onset of end-stage renal disease in a polycystic kidney disease model by suppressing the increases in kidney volume and renal injury[END_REF][START_REF] Boertien | Short-term Effects of Tolvaptan in Individuals With Autosomal Dominant Polycystic Kidney Disease at Various Levels of Kidney Function[END_REF][START_REF] Kelsey | Polycystic kidney disease: Tolvaptan in ADPKD-TEMPO 3:4 trial results[END_REF].

Selective peptide antagonist of V2R were developed [START_REF] Manning | Peptide and non-peptide agonists and antagonists for the vasopressin and oxytocin V1a, V1b, V2 and OT receptors: research tools and potential therapeutic agents[END_REF][START_REF] Manning | Oxytocin and Vasopressin Agonists and Antagonists as Research Tools and Potential Therapeutics[END_REF] but these efforts have encountered many obstacles due to the residual agonist activity, heterogeneity in species response and very low oral bioavailability limiting their clinical use [START_REF] Rondon-Berrios | Vasopressin receptor antagonists: Characteristics and clinical role[END_REF]. These limitations make the development of new non-peptide V2R antagonists more attractive.

Orally and intravenously active non-peptide vasopressin receptor antagonists are called vaptans.

The first success in this field was mozavaptan (OPC-31260), a benzazepine derivative and a potent, selective, competitive and orally active vasopressin V2 receptor antagonist [START_REF] Yamamura | Characterization of a novel aquaretic agent, OPC-31260, as an orally effective, nonpeptide vasopressin V2 receptor antagonist[END_REF], soon followed by its use in humans [START_REF] Ohnishi | Potent aquaretic agent. A novel nonpeptide selective vasopressin 2 antagonist (OPC-31260) in men[END_REF] and the first to gain approval for clinical use in Japan since 2006 for the treatment of tumor-associated to Syndrome of inappropriate Antidiuretic Hormone secretion (SIADH).

Among the non-peptide V2R antagonists developed and experimentally tested, only two compounds of this class have been approved in the United States, Canada and the European Union [START_REF] Rondon-Berrios | Vasopressin receptor antagonists: Characteristics and clinical role[END_REF]. The U.S. Food and Drug Administration (FDA) approved conivaptan and tolvaptan for euvolemic and hypervolemic hyponatremia [START_REF]Drugs@FDA: FDA-Approved Drugs[END_REF]. Tolvaptan is also approved to slow kidney function decline in adults at risk of rapidly progressing Autosomal Dominant Polycystic Kidney Disease (ADPKD), the only drug approved to treat this condition so far [START_REF]Drugs@FDA: FDA-Approved Drugs[END_REF].

Besides the retardation of progressive renal failure in ADPKD and the treatment for euvolemic or hypervolemic hyponatremia, experiments show that V2R antagonists can be used for rescue treatment in Congenital Nephrogenic Diabetes Insipidus [START_REF] Ranieri | Vasopressin-aquaporin-2 pathway: recent advances in understanding water balance disorders[END_REF], treatment of diabetic nephropathy [18], congestive heart failure [START_REF] Izumi | Therapeutic potential of vasopressin-receptor antagonists in heart failure[END_REF], and also in the prevention of ascites formation in cirrhosis [START_REF] Gassanov | Arginine vasopressin (AVP) and treatment with arginine vasopressin receptor antagonists (vaptans) in congestive heart failure, liver cirrhosis and syndrome of inappropriate antidiuretic hormone secretion (SIADH)[END_REF].

Other indications for treatment with vasopressin-receptor antagonists will probably emerge.

In this work we have used a computational modeling approach to study a family of non-peptide V2R antagonists, with their IC 50 being determined experimentally, identifying physicochemical properties relevant to the biological activity of these compounds and relevant interactions with V2R.

METHODS

Data set of V2R Antagonist

To obtain a reliable QSAR model we used chemical information from the assays AID- binding site and the inhibition of intracellular cAMP accumulation. The IC 50 value is the concentration of compound which inhibits [3H]-AVP binding by 50%. In our study, the negative logarithm of the biological activity, pIC 50 , was used as the dependent variable to determine QSAR correlation equations. The 3D structure of each antagonist were generated from its SMILES in Pubchem database.

Estimation of molecular properties

The specific action of drugs depends on many intrinsic features such as hydrophobic, electronic, and steric properties. In a QSAR model, the biological activity is expressed as a function of molecular descriptors. A molecular descriptor encodes as a number, the result of a mathematical and logical procedure using the information of specific properties of molecules. In this study we calculated as a hydrophobic descriptor, the logarithm of the octanol/water partition coefficient (LogP) and hydration energy; as steric descriptors: approximate surface area (ASA), grid surface area (GSA), molar volume (MV) [START_REF] Bodor | A new method for the estimation of partition coefficient[END_REF][START_REF] Hasel | A rapid approximation to the solvent accessible surface areas of atoms[END_REF], and molar refractivity (MR); as for electronic descriptors: polarizability (P) [START_REF] Miller | Additivity methods in molecular polarizability[END_REF], dipole moment (μ), total energy (TE), highest occupied molecular orbital eigenvalue (eHOMO), lowest unoccupied molecular orbital eigenvalue (eLUMO), partial atomic charges of the pharmacophore atoms (q1 to q11), electrophilicity index (ω), chemical hardness (η), chemical softness (s) [START_REF] Parthasarathi | Electrophilicity index as a possible descriptor of biological activity[END_REF]. Electronic descriptors were calculated by Kohn-Sham's DFT B3LYP/6-31G method as included in Gaussian 09 program routines [START_REF] Frisch | HyperChem: a software package for computational chemistry and molecular modeling[END_REF]. The other descriptors were calculated with QSAR properties available in Hyperchem v8 software [26].

Cluster analysis

Cluster analysis is used in QSAR models to build the training and test sets as well as to determine the structural diversity of the dataset. In cluster analysis, the antagonists were classified in groups, called clusters, with a relative homogeneity. The structural diversity or similarity between the compounds is determined by calculating the Euclidean distance between each couple of objects: the smaller the distance, the more of the objects are considered similar to each other [START_REF] Mcfarland | Multivariate Data Analysis of Chemical and Biological Data. Cluster Significance Analysis[END_REF]. To check the structural diversity of the dataset and to define the number of possibles clusters, a hierarchical cluster analysis of these molecules was performed using k-NNCA algorithm to construct the dendrogram. The complete linkage distance (Euclidean metric) was used as the connection function to merge the objects into clusters. The complete linkage measures the proximity between two groups, calculating the distance among the farthest objects, or the similarity among the objects with lesser similarities. The Euclidean distance is the square root of the sum of the squared differences among the values of two objects for each variable [START_REF] Alvarez-Ginarte | Integration of ligand and structure-based virtual screening for identification of leading anabolic steroids[END_REF][START_REF] Yu | Optimized Data Fusion for Kernel k-Means Clustering[END_REF].

To select the training and test sets we used the k-mean cluster algorithm (k-MCA). Such algorithms use a switching method to divide N data points into k groups (clusters) to minimize the sum of distances/dissimilarities among the objects within the same cluster. The k-mean approach requires that k (the number of clusters) must be known before clustering [START_REF] Yu | Optimized Data Fusion for Kernel k-Means Clustering[END_REF]. The k values were set taking into account the dendrogram obtained for the first cluster analysis. Both hierarchical and partitional (non-hierarchical) cluster analyses were implemented using the STATISTICA 8 software [START_REF] Weiß | STATISTICA, Version 8[END_REF]. After the cluster analysis, the compounds were separated in two sets: 80% of compounds in each cluster were selected for the training set and 20% of each cluster for the test set. The training set was used to develop the QSAR model and the test set was used for external cross-validation of the model.

QSAR model

A correlation matrix was performed to determine among the calculated molecular descriptors the ones that do not correlate to each other. A Genetic Algorithm (GA) was used as a metaheuristic method for the molecular descriptors selection and optimization of the functions [START_REF] Alvarez-Ginarte | A novel in-silico approach for QSAR Studies of Anabolic and Androgenic Activities in the 17β-hydroxy-5αandrostane Steroid Family[END_REF][START_REF] Liu | Current mathematical methods used in QSAR/QSPR studies[END_REF][START_REF] Pourbasheer | QSAR study of IKKβ inhibitors by the genetic algorithm: multiple linear regressions[END_REF][START_REF] Pourbasheer | Quantitative structure activity relationship study of p38α MAP kinase inhibitors[END_REF]. The length of the equation was set for three terms and a constant, and the GA was used for input selection to establish which of the descriptors will have the best multiple linear regression (MLR).

Several statistical parameters were employed to validate the model. A good QSAR model should have the highest squared correlation coefficient, R² and Fisher-test, with the lowest standard deviation (S). The P-value is another important parameter used for modeling validation and it should be lower than 0.01. The predictive power of the model was then determined by examining the leave-one-out (LOO) cross-validation (q²). The q² is known as the predictive variance, with a value higher than 0.5. To validate the QSAR model, an external prediction test set of compounds (in the model range) was used, as the predictive ability of a QSAR model shall only be estimated using an external test set [START_REF] Tropsha | Best Practices for QSAR Model Development, Validation, and Exploitation[END_REF][START_REF] Gramatica | Principles of QSAR Modeling: Comments and Suggestions From Personal Experience[END_REF]. All the procedures used to build the QSAR model were performed with BuildQSAR software [START_REF] Oliveira Db De | BuildQSAR: A New Computer Program for QSAR Analysis[END_REF] and validated with STATISTICA 8 software [START_REF] Weiß | STATISTICA, Version 8[END_REF].

The applicability domain

The applicability domain is the theoretical region of the chemical space, defined by the model descriptors and the modeled response, therefore by the nature of the compounds in the training and test sets, as represented in each model by specific molecular descriptors. The applicability domain of a QSAR model is "the range within which it tolerates a new molecule" [START_REF] Alvarez-Ginarte | Integration of ligand and structure-based virtual screening for identification of leading anabolic steroids[END_REF], a QSAR model is only valid within the same domain for which it was developed. Even if the models are developed for the same chemical structures, the applicability domain for new structures can differ from model to model, depending on specific descriptors.

In multiple predictor models, performing simple single-variable range checks is not sufficient to verify the applicability domain. For MLR, one of the most used approaches with normally distributed data for a multiple predictor problem is a distance-based measure like the leverage (h).

As the leverage of a compound measures its influence on the model, it becomes possible to verify whether a new chemical will fit within the structural model domain. The leverage used as a quantitative measure of model applicability domain is also suitable for evaluating the degree of extrapolation, which represents a sort of compound distance from the model experimental space.

The warning leverage (h*) is a critical value or cut-off. Predictions should be considered unreliable for compounds with high leverage (h>h*, being the critical value h*=3p'/n, where p' is the number of model variables plus one, and n is the number of the compounds) [START_REF] Alvarez-Ginarte | Integration of ligand and structure-based virtual screening for identification of leading anabolic steroids[END_REF].

V2R modeling

GPCRdb template tool (https://gpcrdb.org/structure/template_selection) was used to identify the possible templates for V2R and the human OX2 orexin receptor (PDB ID 4S0V) was selected for V2R modeling. The alignment between the template sequence and V2R was made with a Clustal X v2.1 profile [START_REF] Larkin | Clustal W and Clustal X version 2.0[END_REF] taking into account the alignment already obtained with the GPCRdb template tool. The model was built by homology with YASARA program v12.8.26 [START_REF] Krieger | Increasing the precision of comparative models with YASARA NOVA--a self-parameterizing force field[END_REF]. YASARA uses knowledge-based energies to validate the receptor model normalizing them to remove the dependencies on the size and shape of the protein, and also on its amino acid composition, obtaining estimates for the expected average energy and its standard deviation from gold-standard reference structures. Then it calculates how many standard deviations it is away from the average, thereby obtaining a Z-score to evaluate the quality of the models.

To minimize and equilibrate the receptor model, it was inserted in a 1-palmitoyl-2-oleoyl-snglycero-3-phosphocholine (POPC) membrane patch, generated with the VMD [START_REF] Humphrey | VMD: visual molecular dynamics[END_REF] Build Membrane plugin, which mimics its natural environment. POPC was selected because its abundance in biological membranes and because it does not introduce any curvature in the structure. The receptor model was oriented on the membrane, according to the orientation of the orexin receptor (ID PDB 4S0V) on the membrane in the MemProtMD database [START_REF] Stansfeld | MemProtMD: Automated Insertion of Membrane Protein Structures into Explicit Lipid Membranes[END_REF]. The membrane with the receptor inserted was oriented in the XY plane. The system was solvated in a 93 × 92 × 113 Å periodic box of TIP3P water and NaCl was added at physiological concentration, neutralizing the system. The energy of the system was minimized with 1000 steps of conjugate gradient and a further equilibration for 10ns (0.5ns to 300K lipids tails; 0.5ns to 300K membrane; 9ns to 310K, on the whole system) was performed with NAMD2 [START_REF] Phillips | Scalable molecular dynamics with NAMD[END_REF] with an integration time step of 2 fs. A Langevin thermostat and barostat were used to maintain an NPT system, the cut-off for non-bonding interactions was 12 Å, a smooth switching function at 10 Å was used for van der Waals interactions, and Particle Mesh Ewald (PME) for electrostatic interactions. The membrane parameters were checked with MEMBPLUGIN [START_REF] Guixà-González | MEMBPLUGIN: studying membrane complexity in VMD[END_REF] of VMD.

Molecular docking

The receptor model obtained in the previous simulation was used to perform docking studies to predict the binding modes of three antagonist of the studied family. The receptor model and the antagonists were prepared using AutoDockTools [START_REF] Morris | AutoDock4 and AutoDockTools4: Automated Docking with Selective Receptor Flexibility[END_REF] to perform the molecular docking with

AutoDock Vina [START_REF] Trott | AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading[END_REF], where the antagonist and the sidechains of the receptor residues (Q96, W99, F105, K116, and F307) were flexible. The search space was restricted to a 28 x 20 x 14 Å box. The default parameters for configuration files were used in Autodock Vina, running it 5 times and saving 10 conformations of each compound for each run to generate a total of 50 conformations for each compound. The docking results were visually analyzed using UCSF Chimera [START_REF] Pettersen | UCSF Chimera--a visualization system for exploratory research and analysis[END_REF]. All docking results were clustered using a tolerance value of 2.0 Å RMSD and three representative orientations in the binding site were identified to select one conformation per compound. The interactions between the selected antagonist conformation and the V2R were analyzed using BINANA [START_REF] Durrant | BINANA: A Novel Algorithm for Ligand-Binding Characterization[END_REF].

Molecular dynamic simulations of complexes in POPC

The topologies and parameters of the selected antagonists for CHARMM 36 force field [START_REF] Best | Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone φ, ψ and side-chain χ1 and χ2 dihedral angles[END_REF] were generated using CHARMM-GUI [START_REF] Jo | CHARMM-GUI: A web-based graphical user interface for CHARMM[END_REF][START_REF] Lee | CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM Simulations Using the CHARMM36 Additive Force Field[END_REF]. The complex in the POPC membrane was minimized (10000 steps) and equilibrated for 100 ps at 310K with NAMD2. A molecular dynamic simulation was performed for 50 ns (310K, NTP and constant area) with NAMD2, saving frames and calculating energy every 5000 steps.

For the equilibration and the production simulations, a Langevin thermostat and barostat were used to maintain an NPT system, the cut-off for non-bonding interactions was 12 Å, a smooth switching function at 10 Å was used for van der Waals interactions and Particle Mesh Ewald (PME) for electrostatics interactions. The integration time step was 2 fs. VMD [START_REF] Humphrey | VMD: visual molecular dynamics[END_REF] was used for the analysis and visualization of the molecular dynamic simulations.

Complex free energy calculations using linear interaction energy methods

The linear interaction energy (LIE) method was used to estimate the free energy of antagonistreceptor binding. For this purpose, in addition to the previous simulation of the complex in the membrane in a cubic water box, a second simulation of the antagonist only in the water box is needed, which was carried out using the same parameters as the simulation of the antagonistreceptor complexes. Eq.1 shows the improved LIE formula suggested by Almlöf et al. [START_REF] Almlöf | Improving the Accuracy of the Linear Interaction Energy Method for Solvation Free Energies[END_REF][START_REF] Miranda | Improving the LIE Method for Binding Free Energy Calculations of Protein-Ligand Complexes[END_REF] and takes into account the intra-ligand electrostatic interactions.

(

Where 〈V el l-S 〉and 〈V vdw l-S 〉are MD-generated interaction energy averages from the nonbonded electrostatic and van der Waals interactions of the ligand with its surrounding environment (s). 〈V el l-l 〉 is the electrostatic intramolecular ligand-ligand average energy. The ∆'s denote the change in average values when transferring the ligand from solution (free state) into the binding site of the solvated receptor (bound state). Coefficients α and β are scaling factors for the energy terms, while γ is an empirical constant. In this study, α was considered as 0.18 which is considered to be a robust value from previous works [START_REF] Almlöf | Improving the Accuracy of the Linear Interaction Energy Method for Solvation Free Energies[END_REF][START_REF] Miranda | Improving the LIE Method for Binding Free Energy Calculations of Protein-Ligand Complexes[END_REF][START_REF] Hansson | Ligand binding affinity prediction by linear interaction energy methods[END_REF]. The β specific values for each antagonist was calculated using the parameterization model E proposed by Almlöf et al. [START_REF] Almlöf | Improving the Accuracy of the Linear Interaction Energy Method for Solvation Free Energies[END_REF][START_REF] Miranda | Improving the LIE Method for Binding Free Energy Calculations of Protein-Ligand Complexes[END_REF] (Eq.2).

(

) 2
Where w i , β 0 , and Δβ i , were calculated from explicit solvent FEP calculations of single chemical group (w i = 1 if group is neutral and 11 if it is an anion or a cation), β 0 = 0.43 and Δβ i was obtained by the model proposed by Almlöf et al [START_REF] Almlöf | Improving the Accuracy of the Linear Interaction Energy Method for Solvation Free Energies[END_REF].

The balance (difference) between the electrostatic (polar) and the van der Waals (nonpolar) contributions to the free energy binding in the LIE method was defined as the parameter D (Eq. 3)

(

LIE-D is an approach based on the linear correlation between the γ coefficient and the D parameter that accounts for the balance (difference) between the polar and nonpolar binding free energy contribution. The relationship between the γ coefficient and D parameter takes the form:

The values of f and g were estimated by Miranda et al. [START_REF] Miranda | Improving the LIE Method for Binding Free Energy Calculations of Protein-Ligand Complexes[END_REF] as -0.95 and -2.06 respectively.

RESULTS AND DISCUSSION

Construction of training and test sets using Cluster Analysis

We selected a series of 53 compound antagonists of V2R to construct the training and test sets. All the selected molecules have the same core substructure

γ = f x D + g [kcal/mol]
compounds A03 and A07. We display in Table 1 the compounds we used in this study with PubChem ID and the experimental biological activities of V2R antagonists' activity expressed as IC 50 and pIC 50 .

To classify the molecules of the datasets, depending on their structural variability we performed a hierarchical cluster analysis, the resulting dendrogram was constructed using the Euclidean distance (x-axis) and the complete linkage (y-axis), illustrating the results of the k-NNCA developed in this dataset. The dendrogram shows 6 different subsets demonstrating the molecular variability among the compounds of this dataset (Figure 2). To evaluate the output dendrogram and to split the whole dataset into training and test sets, we performed a k-mean cluster analysis (k-MCA) [START_REF] Yu | Optimized Data Fusion for Kernel k-Means Clustering[END_REF].

The selection of the training and test sets was carried out by randomly taking molecules belonging to each cluster. From the initial 53 compounds, 42 (80% of the dataset) were chosen to form the training set and the remaining 11 compounds, (20% of the dataset) were used as a test set for the external cross-validation of the model.

Development and validation of the QSAR model

GA combined with MLR is widely used for QSAR and QSPR studies [START_REF] Alvarez-Ginarte | A novel in-silico approach for QSAR Studies of Anabolic and Androgenic Activities in the 17β-hydroxy-5αandrostane Steroid Family[END_REF][START_REF] Liu | Current mathematical methods used in QSAR/QSPR studies[END_REF][START_REF] Pourbasheer | QSAR study of IKKβ inhibitors by the genetic algorithm: multiple linear regressions[END_REF][START_REF] Pourbasheer | Quantitative structure activity relationship study of p38α MAP kinase inhibitors[END_REF]. In this method, a GA is performed to search the feature space and select the major descriptors relevant to the activities or properties of the compounds. This method can deal efficiently with a large search space, and it has fewer chances to only find a local optimal solution than other algorithms. GA is a well-estimated method for parameter selection and to overcome the shortages of MLR in variable selection. After a GA, the MLR is employed to correlate the selected descriptors with the activity values using a classic regression method to yield the explicit equations. The R 2 (R-square statistic or coefficient of determination) indicates that the model could explain 80 % of the variance for the experimental values of pIC 50 . The model shows a q 2 of 0.75. This value of more than 0.5 could be considered as proof of the high predictive ability of the model, along with the good prediction of the test set (R 2 = 0.74). The good R 2 and q 2 values obtained in equation 5 for both training and test set can be explained with the experimental values for all the compounds of the series. The calculated values for pIC 50 , are highly similar to the experimental, sustaining the reliability of the QSAR model (Figure 3, Table 2).

In the correlation study with the calculated descriptors, a low correlation was observed between the variables, indicating the reliable information content on each term in the equation (Table 4). The selected variables by the genetic algorithm were P (polarizability) E H2O (hydration energy), and q6

(partial charge of nitrogen in the common substructure 4-formamidobenzamide (Figure 1)). For each one of the variables, the coefficients were significant (Table 5), indicating their relative contribution to the combined prediction of the biological activity as the dependent variable.

Calculating the value of the coefficients on the regression analysis we ensure a good prediction starting from the group of the independent variables (q6, E H2O and P), facilitating the interpretation of the independent influence of each variable on the final equation.

The variable q6 represents the partial charge of the N7 in the 4-formamidobenzamide common substructure, involved on an amide bond associated with the variable zone of the compounds. The partial charge q6 is the most negative, with the highest module value for all the studied charges, having a fully negative value range and a negative coefficient on equation 5, could indicate the favorable tendency of an increased antagonist activity with more negative values of q6. It could be explained by the fact that N7 is involved in a hydrogen bond, or because it just reflects the variation of partial charge depending on the nature of substituent R1. The partial charge q6, calculated only for the common substructure is -0.712 and the substituent on R1 is making it more negative, except for compounds A53, A24, and A12. In the case of A12, the N7 is substituted by a methyl, which has a strong inductive effect (+i) over the nitrogen. Compounds A53 and A24 both have a benzene ring with a nitro substituting in the ortho position as R1. The common substructure in most compounds is formed entirely by a conjugated system with a substituted benzene ring as R1, which might contribute to a whole conjugated system. If we compare the compounds by the position of the substituent on the benzene ring at R1, we observe a lower partial charge on q6 associated with an ortho substituent, calling for a combination of steric and electronic factors, where the ortho substituent could break the conjugation planarity.

Hydration energy (E H2O ) is the amount of energy released when one mole of a compound is hydrated, and represents the measure of the water molecules affinity for the compound. More negative hydration energy values could be associated with more polar groups in the compound, and less negative hydration energy could be attributed to the presence of a higher number of nonpolar groups [START_REF] Schauperl | Enthalpic and Entropic Contributions to Hydrophobicity[END_REF]. The range of values for this variable in the data set is negative and it has a positive coefficient on equation 5, thus indicating that more negative values of hydration energy are unfavorable for the antagonist activity, suggesting a binding site with possible hydrophobic interactions, as more hydrophilic compounds are shown unfavorable for the activity.

Polarizability refers to the tendency of any compound to acquire an electric dipole moment in proportion to an applied electric field, on our model the Polarizability component is having a correlation coefficient of 0.7 with the pIC 50 , making it the descriptor with the highest correlation to the activity, being the other variables the fine adjustments necessary to improve the model in general. The range of polarizability values in the data set is positive and has a positive coefficient in equation 5, so it has a favorable contribution to antagonist activity.

As mentioned above, the compounds in the studied dataset have conjugated systems in their structure, and systems with delocalized π electrons exhibit high polarizabilities. The aromatic systems' planarity with their high polarizability and multipole moment, are all factors of key importance for the 3D architecture of aromatic complexes [START_REF] Riley | On the Importance and Origin of Aromatic Interactions in Chemistry and Biodisciplines[END_REF]. Soft interactions like dispersion, are predominant in stacking and can be estimated from the polarizability [START_REF] Mignon | Influence of the π-π interaction on the hydrogen bonding capacity of stacked DNA/RNA bases[END_REF]. In general, the studied compounds have three aromatic rings in their structure that could be involved in pi-interactions with the residues in the receptor-binding site.

In general, the obtained QSAR model provides indications that the binding mode of V2R antagonists might fundamentally be involving hydrophobic and electron density interactions.

The applicability domain (AD) of the QSAR model

A QSAR model needs to show not only a good accuracy, but also some reliability for predictions of new compounds. These models cannot be universal and should be constrained to a defined chemical space, commonly known as the applicability domain (AD). The AD can be described as the physicochemical, structural or biological spatial information based on which the model training set is developed. The QSAR model applies to make predictions for new compounds within the specific domain [START_REF] Mignon | Influence of the π-π interaction on the hydrogen bonding capacity of stacked DNA/RNA bases[END_REF]; in summary, the AD is the degree to which a QSAR model tolerates (reliably) new compounds.

A crucial problem in chemometrics and QSAR studies is the definition of the AD with a regression model. We will define it here as a squared area within ± 2 bands for standardized residuals and a leverage threshold of h=0.23 for inhibitory activity (Eq. 5). Thus, compounds with standardized residuals greater than 2 standard deviations will be considered unreliable. For the graphical visualization of outliers for the response (standardized residuals > 2) or for the structure (leverage > 0.23) in the regression model, the Williams plot for Eq. 5 is shown in Figure 4. Of the 53 compounds in the dataset, only two compounds (A02 and A12) have a leverage higher than the critical value.

A02 (conivaptan) has the highest value of polarizability (57.51) of the dataset, while the other compounds are between 36.87 to 55.95. A02 shows a diphenyl moiety as substituent of the amide in the common substructure, while the other compounds exhibit only a single aromatic ring or an aliphatic substituent, A02 also have a condensed 3 ring system of 3,4,5,6-tetrahydroimidazol [4,5d][1]benzazepine, while the other compounds have only a 2 ring system. The presence of extra rings on A02, might account for the increase on the polarizability for this compound with a different electronic structure than the rest.

A12 shows the highest value of q6 (-0.515) of the dataset, while the other compounds are between -0.705 to -0.844, and it also has the highest value of hydration energy (E H2O ), with -3.44 kcal/mol, while the other compounds are between -5.05 to -14.07 kcal/mol. Compound A12 show a minor difference in the common substructure with all the antagonists of this family having a methyl group as substituent for the N amide, directly altering the partial charge (q6) of N7 on the 4formamidobenzamide and increasing the hydration energy value. The structure of A12 and A01 differ only in the aforementioned methyl group, A12 has a difference of 0.203 on q6, and of 2.36 kcal/mol in the value of hydration energy compared to A01, evidencing the influence that this single methyl group can have.

V2R modeling

We selected the human OX2 orexin receptor (PDB ID 4S0V) as the template for V2R modeling, as suggested by the GPCRdb template tool. The selected template has 27% of identity and 46% of similarity with V2R. In the corresponding alignment, the fragments corresponding to transmembrane helices and the conserved motifs are preserved (Figure 5). In both the OX2 receptor and V2R, the natural ligand is a peptide, and the selected structure has an antagonist bound being on an inactive conformation, suitable to study the binding modes of antagonists to the V2 receptor.

To relax the obtained model in a more natural environment, it was minimized and equilibrated for 10ns in a POPC membrane, solvated and with NaCl added at physiological concentration (0.15M).

At the end of the simulation the membrane parameters, like thicknesses and per lipid area, were calculated to check the correct packing. The membrane thickness (distances between phosphates of each monolayer) was 38.91 Å, and the per lipid area was 64.54 Å 2 . These parameters are reasonable for a POPC membrane at 310K, according to experimental parameters obtained at different temperatures [START_REF] Kučerka | Fluid phase lipid areas and bilayer thicknesses of commonly used phosphatidylcholines as a function of temperature[END_REF].

During the current work, another X-ray structure from the class A of GPCR was released in the Protein Data Bank (PDB ID 6TPK): the oxytocin receptor (OXTR) which is also a member of the vasopressin receptor family. Although OXTR exhibits slightly better sequence identity and similarity: 41% and 56% respectively, its lower resolution (3.20Å versus 2.50Å), a missing region (loop and helix 8), and a shorter loop (…) make it a less relevant template. However, in an effort for further validation, the V2R model was compared to OXTR; the superimposition between the V2R model and OXTR is shown in Figure 6, where we considered the two conformations of V2R before and after membrane relaxation.

The RMSD values, using OXTR as reference and comparing the model before and after the membrane relaxation were 1.02Å and 1.15Å respectively. The main differences between OXTR and the models were that OXTR lacks the ICL3, a long intracellular loop involved in the interaction with the G-protein usually missing in GPCR solved structures and the helix 8, parallel to the membrane and useful for orienting the receptor in the membrane. Comparing the bundled helices, the main difference relative to the binding site is that the TM2 of relaxed model in the membrane is in the same position as that of the OXTR, while the model before relaxing has the TM2 slightly tilted towards the interior of the cavity decreasing its volume. The difference in the orientation of TM2 in the model before relaxation could be caused by the difference in the proline position in this transmembrane section between V2R and the template (Figure 5). Proline in the middle of alpha helices cause a kink, by being unable to complete the H-bonding chain of the helix and because of steric and/or rotameric effects keeping it out from the preferred helical geometry [START_REF] Von Heijne | Proline kinks in transmembrane α-helices[END_REF]. Proline one position earlier on the sequence of the template TM2 with respect to that of the model, can affect the orientation of the kink and result on a different orientation in the model. This odd orientation could later be corrected during the relaxation in the membrane showing this process as very favorable.

The comparison between the V2R model after relaxation and OXTR brings more confidence in the model's quality and the protocol used for relaxation.

V2R-antagonist complexes

Visual inspection of the binding site revealed that the side chains of residues W99 and F307 are occluding its entrance, therefore these two residues were considered as flexible for the molecular docking. To improve docking results other residues of the binding site (Q96, F105 and K116) were also considered as flexible.

The antagonists selected as ligands for the molecular docking were: mozavaptan (A01), conivaptan (A02), and tolvaptan. Tolvaptan also shares the common substructure of the studied compound series for the QSAR model and it is the only drug approved to treat Polycystic Kidney Disease. All the rotatable bonds of the ligands were flexible.

The 50 complexes for each antagonist obtained by molecular docking were clustered for analysis.

Three orientations of the antagonists in the binding site were identified and shown in Figure 7. The first with the condensed ring of the antagonist toward TM2 and TM7 (CR-27), the second with the condensed ring towards TM5 and TM6 (CR-56) and the last one with the condensed ring towards the entrance of the cavity (CR-UP). The binding energies of the complexes obtained for the three identified conformations differ by less than 1 kcal/mol in each of the antagonists studied. This difference in the energy value is lower than the standard error of the Autodock Vina scoring function [START_REF] Trott | AutoDock Vina: Improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading[END_REF], so we might need more studies to select the best conformation.

We expect a good antagonist to bind with high affinity to the receptor binding site but failing to activate it, blocking the access of any agonist to the binding site. A study with meta-dynamics enhanced sampling revealed the existence of three binding sub-sites for V2R, proposed to respond to the vasopressin entry pathway [START_REF] Saleh | A Three-Site Mechanism for Agonist/Antagonist Selective Binding to Vasopressin Receptors[END_REF]. The compounds that bound in both the vestibule and the intermediate sites block the access to the orthosteric site so that an agonist will never be able to bind, if there is an antagonist already bound to any of the non-activating sites. Two of the antagonists studied by Saleh et al. [START_REF] Saleh | A Three-Site Mechanism for Agonist/Antagonist Selective Binding to Vasopressin Receptors[END_REF], with high structural similarity to those in this study, were predicted to bind to vestibule site for V2R and intermediate site for V1aR, so it is to be expected that the antagonists in our study are located in one of these sites. Therefore, we eliminated from our subsequent analysis the CR-UP conformation, showing the antagonist penetrating deeper into the cavity with part of it located in the orthosteric site.

To study the stability of the compounds in the binding site and make a better estimation of the binding energy, a molecular dynamic simulation of the best complex of the two remaining orientation was performed. Figure 8 In order to predict which binding conformation could be the best for each antagonist, we estimated the binding free energy variation using LIE-D method. This method is flexible enough to consider different interaction patterns, even though the ligands share some common chemical scaffolds and are bound to the same protein receptor [START_REF] Miranda | Improving the LIE Method for Binding Free Energy Calculations of Protein-Ligand Complexes[END_REF].

For tolvaptan, the conformation with the best binding free energy is tolvaptan CR-27 (-14.34 kcal/mol) with overn 2 kcal/mol of difference with CR-56, that was thus discarded from further analysis. The conformations of conivaptan have similar energies between them, so it is not possible to select which of the two conformations is the best using this criterion. The inhibition constant (K i )

was calculated from the estimated binding free energy for these conformations and compared with experimental values (Table 6). The values of the calculated K i follow the same trend than the experimental inhibition constants [START_REF] Crombie | Synthesis and evaluation of azabicyclo[3.2.1]octane derivatives as potent mixed vasopressin antagonists[END_REF][START_REF] Yamamura | OPC-41061, a highly potent human vasopressin V2-receptor antagonist: pharmacological profile and aquaretic effect by single and multiple oral dosing in rats[END_REF], except for the tolvaptan, with the K i for conformation CR-27 one order lower than the experimental value and the K i conformation CR-56 one order higher, suggesting as expected that CR-27 was the best conformation for tolvaptan.

The binding free energy and K i values obtained for the different conformations of the antagonists bound to V2R allowed us to select the CR-27 conformations (or the CR-56 conformation for conivaptan only) as the possible binding modes of these antagonists to V2R. From these conformations, the analysis of the main antagonist-receptor complex interactions can be carried out.

Interactions analysis of the best antagonist-V2R complexes

The antagonists must block the access or interact with those residues favoring the union of any agonist for the receptor activation, and/or also sterically blocking the residues involved in triggering the activation mechanism. The most relevant contacts between the antagonist and the receptor are summarized in Table 7. The interactions observed for the mozavaptan and tolvaptan complexes are very similar, while conivaptan interacts with a greater number of residues, since it is a compound with a greater volume than those mentioned above.

In the analysis of the obtained complexes, there are some common interaction involving the hydrophobic (C192, A194, L310, M311), aromatic (W99, F105, F178, F307) and polar (Q96, K116, Q119, Q291) residues. The residues Q96, Q119, Q291, K116 are highly conserved in all the AVP and OXT receptors' family, and are known to have a key role in agonist binding [START_REF] Mouillac | The binding site of neuropeptide vasopressin V1a receptor. Evidence for a major localization within transmembrane regions[END_REF]. Previous studies with V1aR suggest that the residues Q96, Q119, Q291 are K116 are specifically involved with the ligand binding process but do not intrinsically modulate the efficacy of the functional response [START_REF] Mouillac | The binding site of neuropeptide vasopressin V1a receptor. Evidence for a major localization within transmembrane regions[END_REF]. An analysis of the presence of H-bonds along the molecular dynamic simulation was performed to study the nature of the interactions between the antagonists and these polar residues of the receptor. The percentage of frames of the trajectory with a determined number of H-bonds is shown in the Table 8. In all the studied conformations the occurrence of H-bonds is low, approximately between 13 and 38%, although interactions with some polar residues are observed in greater percentages of the trajectory (Table 7), which suggests other kind of interactions.

The antagonists interact with the aromatic residues W99, F105, F178 and F307. These residues are not directly involved in the receptor activation and they are near or at the entrance of the binding site, suggesting that the antagonist binding site might not be located deep into the cavity. The fact that the antagonists interact with residues near the entrance of the cavity agrees with the binding sites predicted by Saleh et al. [START_REF] Saleh | A Three-Site Mechanism for Agonist/Antagonist Selective Binding to Vasopressin Receptors[END_REF].

Some of these residues have been shown to be involved in vasopressin binding. W99 plays a fundamental role in stabilizing the vasopressin/receptor interactions responsible for the high-affinity binding of agonists to the V2 receptor and receptor selectivity. A mutation of W99 (W99R) greatly impaired the binding properties of the receptor and had a minor effect on its intracellular routing [START_REF] Albertazzi | Nephrogenic Diabetes Insipidus: Functional Analysis of New AVPR2 Mutations Identified in Italian Families[END_REF]. Other important residue for AVP binding is F105, the mutation F105V was reported to show cell surface expression and a maximal AVP-induced cAMP formation (Vmax) comparable to the wild type, but with a reduced ligand binding ability [START_REF] Pasel | Functional characterization of the molecular defects causing nephrogenic diabetes insipidus in eight families[END_REF][START_REF] Makita | V2 vasopressin receptor mutations[END_REF].

An interesting interaction for the V2R antagonists is with F307, a non-conserved residue in vasopressin/oxytoxin family since V1aR has a threonine in this position. The relevance of this interaction is because some antagonists could bind to both V2R and V1aR due to the similarity of its binding site, but the interaction with F307 would be unique for V2R making it attractive for the design of antagonists having less selectivity for V1aR.

Other antagonist-receptor interactions found were with residues C192, A194 and M311. While C192 and A194 are conserved among the entire family, M311 is not seeming to cooperate in the selective binding of some antagonists [START_REF] Cotte | Conserved aromatic residues in the transmembrane region VI of the V1a vasopressin receptor differentiate agonist vs. antagonist ligand binding[END_REF]. The M311V mutation in the TM7 of V2R has impaired the ligand capacity and binding [START_REF] Neocleous | Identification and characterization of a novel X-linked AVPR2 mutation causing partial nephrogenic diabetes insipidus: A case report and review of the literature[END_REF] suggesting that in V2R the residue M311 could take part in the binding of peptide agonists [START_REF] Makita | V2 vasopressin receptor mutations[END_REF][START_REF] Sahakitrungruang | Functional characterization of vasopressin receptor 2 mutations causing partial and complete congenital nephrogenic diabetes insipidus in Thai families[END_REF] Taking into account the interactions and the estimated binding free energy, we considered CR-27 conformation (and CR-56 for conivaptan) as the best for the three antagonists studied by molecular dynamics. Figure 10 shows the CR-27 conformation of tolvaptan in complex with V2R.

In general, the main interactions observed here are involved in the binding of ligands to V2R, but are not involved on the receptor activation, which suggest that the studied conformations of the antagonists can block the binding of agonists and unable to activate the receptor. The presence of few H-bonds with polar residues and the other interactions observed are with aromatic residues and non-polar residues suggest that the main antagonist-receptor interactions are mostly hydrophobic in nature and could involve pi-clouds.

CONCLUSION

In summary, two computational approaches, ligand and receptor based were developed to study the physicochemical properties relevant to the biological activity of V2R antagonists and to predict their binding mode to V2R. The proposed QSAR model allows us to clarify the contribution of three molecular descriptors to the biological activity. Our model described the antagonist activity in correlation with polarizability, hydration energy and partial charge on atom N7, explaining the molecular properties contributing to the antagonist-receptor interaction and relevant to the antagonist activity, which is also in agreement with the binding modes for the complexes obtained by molecular docking and molecular dynamics simulation.

A good quality model based on the structure of OX2 orexin receptor was obtained and used to estimate the antagonist orientations in the binding site of V2R. The conformations of studied antagonist were analyzed by molecular dynamics. In general, the CR-27 conformation is considered as the best conformation for the antagonist binding (through interaction analysis and binding free energy estimation). Most of the relevant interactions observed along the molecular dynamics simulation involve the electronic density by the interaction of the antagonist rings mainly with the aromatic residues (W99, F105, F178 and F307), which is in correspondence with what is expected according to the polarizability variable of our QSAR model. Other relevant interactions are hydrophobic in nature (A194 and M311) which agrees with the expected effect of the hydration energy to the antagonist activity in the QSAR model.

The results obtained by both developed approaches are in fair agreement and contribute to a better understanding of V2R antagonism. These results represent a step forward for the efficient search of potential new V2R antagonist molecules.

Figure 1 ,

 1 but show structure variability due to substituent structural diversity. All compounds have a common substructure (4-formamidobenzamide) remarked in a rectangle in the top panel of Figure 1. The nitrogen of the benzamide in the common substructure (represented the partial charge q1) is part of a benzazepine, benzenepiperidine, or benzoxazine condensed ring. The R1 substituent is generally a ring except for ΔG bind = β (Δ〈V el l-S 〉+ Δ〈V el l-l 〉)+ α Δ〈V vdw l-S 〉+ γ D = β (Δ〈V el l-S 〉+ Δ〈V el l-l 〉) -α Δ〈V vdw l-S 〉 [kcal/mol]

  The variables selected by the genetic algorithm as the best model of V2R antagonist activity are shown in equation 5. To further validate the variables thus obtained, we performed an MLR analysis of the 43 compounds on the initial training set, with the 11 compound test set for the external crossvalidation. pIC 50 = -7.968 (± 3.584) q6 + 0.095 (± 0.059) E H2O + 0.161 (± 0.027) P -5.842 (± 2.894) (5) n = 43; R = 0.89; R 2 = 0.80; s = 0.40; F = 53,61; p< 0.0001; q 2 = 0.75 Test set: n = 11; R = 0.86; R 2 = 0.74; s = 0.41; F = 25.04; p< 0.0007; q 2 = 0.56

  shows the RMSD for the two different conformations of the three studied antagonists. The conformations of the studied antagonists tend to stabilize along the molecular dynamics simulations, being CR-27 the conformation with lower RMSD for each of the antagonist. For mozavaptan, the conformation CR-56 is the conformation with more fluctuations along the trajectory and it has the highest RMSD value among all antagonist conformations. The change in the antagonist's conformation for the six complexes in the molecular dynamics simulation is shown in Figure9. The two representative conformations for each antagonist are represented from left to right (CR-27 and CR-56 respectively). Tolvaptan is represented in green, conivaptan in blue and mozavaptan in pink. The starting conformation at 0 ns is represented by the light colored ligand and the final conformation at 50 ns by the dark-colored ligand. Mozavaptan CR-56 conformation showed a significant change in the orientation of the antagonist with respect to the starting conformation, also reflected in the high RMSD value observed for this conformation. For this reason, this conformation was eliminated from the subsequent analysis.

Acknowledgment

The authors thank the assistance from the Department of Applied Physical Chemistry of Autonomous University of Madrid, Spain, and the support of the Faculties of Chemistry and Biology from the University of Havana, Cuba. This research was supported through highperformance computing resources provided by the I2BC CNRS-CEA-Université Paris Sud, Paris Saclay, France.

Conflicts of interest

The authors confirm the article content has no conflicts of interest.


Yoshihiro Shibuya

Djamal Belazzougui

Gregory Kucherov

Set-Min sketch: a probabilistic map for power-law distributions with application to k-mer annotation

Keywords: sketching, Set-Min sketch, Max-Min sketch, k-mer counting, k-mer spectrum, power-law distribution

k-mer counts are important features used by many bioinformatics pipelines. Existing k-mer counting methods focus on optimizing either time or memory usage, producing in output very large count tables explicitly representing k-mers together with their counts. Storing k-mers is not needed if the set of k-mers is known, 1 1 INTRODUCTION making it possible to only keep counters and their association to k-mers. Solutions avoiding explicit representation of k-mers include Minimal Perfect Hash Functions (MPHFs) and Count-Min sketches. We introduce Set-Min sketch -a sketching technique for representing associative maps inspired from Count-Min -and apply it to the problem of representing k-mer count tables. Set-Min is provably more accurate than both Count-Min and Max-Min -an improved variant of Count-Min for static datasets that we define here. We show that Set-Min sketch provides a very low error rate, both in terms of the probability and the size of errors, at the expense of a very moderate memory increase. On the other hand, Set-Min sketches are shown to take up to an order of magnitude less space than MPHF-based solutions, for fully assembled genomes and large k. Space-efficiency of Set-Min in this case takes advantage of the power-law distribution of k-mer counts in genomic datasets.

Introduction

Counting k-mer occurrences in genomic sequences is a rather common task in many bioinformatics pipelines. It is often the first step performed before a more complicated analysis, and its main applications go from read trimming [START_REF] Liu | Musket: a multistage k-mer spectrum-based error corrector for Illumina sequence data[END_REF] to alignment-free variant calling [START_REF] Rahman | Association mapping from sequencing reads using k-mers[END_REF][START_REF] Khorsand | Nebula: Ultra-efficient mapping-free structural variant genotyper[END_REF]. In recent years, many k-mer counting algorithms have been proposed, such as Jellyfish [START_REF] Marçais | A fast, lock-free approach for efficient parallel counting of occurrences of k-mers[END_REF], DSK [START_REF] Rizk | DSK: k-mer counting with very low memory usage[END_REF] or KMC [START_REF] Kokot | KMC 3: counting and manipulating k-mer statistics[END_REF]. Each of these methods was conceived with a particular trade-off between speed and main memory consumption. Jellyfish primarily addresses speed, whereas DSK and KMC primarily target memory efficiency. All these tools output a map associating k-mers to their counts. Such a map can require a fairly big amount of disk space, especially when large values of k are used. For example, the KMC output for a human genome with k = 32 weights in at around 28 GB. Memory efficiency remains an important issue even when outputs are compressed. One way to partially solve this issue is to only store counters, without storing k-mers themselves. This idea is supported by the fact that, in many applications, we only deal with k-mers that come from partially assembled reads or succinct representations of k-mer sets such as colored de Bruijn graphs [START_REF] Holley | Bifrost -Highly parallel construction and indexing of colored and compacted de Bruijn graphs[END_REF], or spectrum-preserving string sets [START_REF] Břinda | Simplitigs as an efficient and scalable representation of de Bruijn graphs[END_REF][START_REF] Rahman | Representation of k-mer sets using spectrumpreserving string sets[END_REF], that is, only k-mers present in the original data are queried for their frequencies.

The idea of storing only counter information and not k-mers is also supported by the observation that the number of distinct k-mer counts in genomic data is relatively small. Furthermore, it is known that k-mer counts in fully assembled genomes obey a "heavy-tail" power-law distribution1 with a relatively large absolute value of the exponent [START_REF] Csűrös | Reconsidering the significance of genomic word frequencies[END_REF][START_REF] Chor | Genomic dna k-mer spectra: models and modalities[END_REF]. For such distributions, the number of distinct k-mers makes a linear fraction of the data size, while the number of distinct k-mer counts is relatively small. For example, for a human genome and k = 27, there are about 2.5 billions distinct k-mers but only about 8,000 distinct frequency values. A heavy-tail distribution also imply that the majority of k-mers have a very small count: in the above example, 97% of k-mers are unique and 99% of k-mers have a count of at most 5. Frequent k-mers often tend to have an identical count as well, due to transposable elements: for example, k-mers specific to Alu repeats in primate genomes will likely have the same count.

Our contribution. We propose a new probabilistic data structure that we call Set-Min sketch capable to represent k-mer counts information in a small 1 INTRODUCTION space and with small errors. The sketch guarantees that the expected cumulative error obtained when querying all k-mers of the dataset can be bounded by εN where N is the number of all k-mers (i.e. essentially, the size of the dataset).

We provide a theoretical analysis in order to dimension the sketch according to the desired error bound. Set-Min sketch is a general data structure in that it can be used to efficiently represent a mapping of k-mers to any type of labels.

Nevertheless, it performs best when data is skewed and the number of possible labels is relatively small.

Applications. The problem considered in this paper is the probabilistic compression of k-mer count tables computed from fully assembled genomes, with large k values. A Set-Min sketch is able to provide a more space-efficient map than other representations with low errors when count distributions are very skewed. Set-Min sketch can have different uses.

An application of Set-Min sketches not explored here, is to act as a temporary representation while building more complex structures based on counters.

Consider, for example, the exact computation of weighted pairwise distances between all pairs of genomes in a given set. Examples of such distances are the Bray-Curtis similarity measure, see e.g. [START_REF] Adamic | Multiple comparative metagenomics using multiset k-mer counting[END_REF], or Weighted Jaccard similarity estimation [START_REF] Chum | Near Duplicate Image Detection: min-Hash and tf-idf Weighting[END_REF]. The most naive algorithm is to first produce count tables for each dataset and then compare them pairwisely to produce the desired output. Instead of storing whole tables, one can store multiple Set-Min sketches together with a presence-absence data structure, such as a Bloom filter. By doing so, the weighted comparison computation is reduced to a single pass through the presence-absence data structure with the counters of a given k-mer retrieved on-demand from the Set-Min sketches of the datasets in which the k-mer is found.

Another application is the quick sharing of counter information between different computational units in a distributed setting. Consider a situation where both a server and one of its clients have access to the same genomic representation (say, a set of contigs). Only the server can efficiently perform k-mer counting, while the client has to process incoming data based on the counts, but doen't have the computational resources to perform counting by itself. In this case, the server could send a Set-Min sketch to the other node.

One further feature of Set-Min sketches is their mergeability in case of redundant maps. In this scenario, a large map can be split into m sub-maps without the restriction of having disjoint sets of keys. Even if some maps have redundant information, i.e. share common (key,value) pairs, the Set-Min sketch built by cell-wise union of the m sketches will be equivalent to the sketch built from the whole original map. Note that Count-Min sketches are not mergeable for redundant maps but are mergeable when constituent maps have to be "added up". Set-Min sketches don't have this property. In this respect, Set-Min and Count-Min sketches may have complementary uses.

Outline of the paper. In Section 2, we start by formally introducing the relevant terminology and underlying methods. Section 3 presents the idea and algorithmic foundations of our method. Our results are presented in Section 4 while Section 5 discuss the limitations of our method and possible workarounds.

We conclude in Section 5 with a summary of the work.

Background

k-mer spectrum

A k-mer spectrum is a distribution of k-mer frequencies across all k-mers occurring in the data, showing how many k-mers support each frequency value.

For large values of k, k-mer spectra follow a power-law distribution (Csűrös 2 BACKGROUND et al., 2007;[START_REF] Chor | Genomic dna k-mer spectra: models and modalities[END_REF] characterized by a linear-like dependence when represented in log-log scale. That is, the k-mer frequency distribution fits a dependence f (t) ≈ c • t -a , where a is usually greater than 2. An example of the spectrum for the human reference genome with k equal to 32 is given in Figure 1. According to this distribution, a very large fraction of k-mers have very low frequencies, while a few k-mers have "unexpectedly" large frequencies.

Large value of a implies that there are relatively few distinct frequency values with non-zero support, whose number depends as the 1 a -power of the number of k-mers. Count-Min sketch (Cormode and Muthukrishnan, 2005a) is a method to compactly represent an associative array a of counters in an approximated way.

Count-Min and Max-Min sketches

Count-Min is especially suitable for the streaming framework, when counters associated to keys can be updated dynamically. A Count-Min sketch is an R × B matrix A of counters where each row i is associated with a hash function

h i (•).
To store counters and support their dynamic updates given in a stream, Count-Min works as follows. To process an update which is a (key,value) pair (p, ), we perform A(i, h i (p)) = A(i, h i (p)) + for each row i. The (approximate) current counter associated with a key p is retrieved as â(p) = min i {h i (p)}.

It has been shown that for a Count-Min sketch built on an associative array (vector) a, with R = ln( 1 δ ) and B = e ε for any given ε and δ, the over-estimate error of an individual counter is bounded by ε a 1 with probability at least 1 -δ, where a 1 is the L 1 -norm of a (Cormode and Muthukrishnan, 2005a). If counts follow a Zipfian rank-frequency distribution with parameter b > 1, B can be reduced to O(ε -1/b ) to guarantee the same bounds (Cormode and Muthukrishnan, 2005b). Note that b > 1 corresponds to 1 < a < 2 in the corresponding spectrum distribution [START_REF] Adamic | Multiple comparative metagenomics using multiset k-mer counting[END_REF], while k-mer spectra often fit a distribution with a ≥ 2. Count-Min sketch supports negative updates, i.e. allows < 0 in updates (p, ), provided that the cumulative value for each key stays positive. If updates are only positive, there exists a modification of Count-Min leading to a better accuracy, mentioned in [START_REF] Cormode | Count-Min Sketch[END_REF] (therein attributed to [START_REF] Estan | New directions in traffic measurement and accounting[END_REF]) as conservative update. Under this modification, updates for each row i are made according to A(i, h i (p)) = max{A(i, h i (p)), â(p) + }, where â(p) is the current Count-Min estimate of a(p). It is easily seen that under this scheme, â(p) can still only over-estimate a(p), but cannot be larger than â(p) computed by the original Count-Min. In this paper, we will deal with the static case when the value of any key is given once and never changes after that. In this framework, the conservative update strategy can be further modified by defining updates as A(i, h i (p)) = max{A(i, h i (p)), }. We call this variant of the sketch Max-Min. In the static case, Max-Min sketch improves Count-Min without any computational overhead: it simply replaces addition by max in the update rule. As with the conservative update, one can check that estimates by Max-Min can only be over-estimates which, however, don't exceed estimates by original Count-Min.

Minimal Perfect Hashing

Minimal Perfect Hash Functions (MPHFs) are bijective functions between keys of a set S and integers in the range [0, |S| -1]. By using hash values as indexes for an external array, it is possible to associate any type of information to the k-mers. The construction of MPHFs can be hyper-graph peeling-based [START_REF] Yu | Memory-efficient and Ultra-fast Network Lookup and Forwarding using Othello Hashing[END_REF][START_REF] Esposito | RecSplit: Minimal Perfect Hashing via Recursive Splitting[END_REF] or array-based [START_REF] Müller | Retrieval and Perfect Hashing Using Fingerprinting[END_REF]. The first family of algorithms leads to smaller MPHFs, close to theoretical space lower bound of 1.44 bits per key, while array-based MPHFs are conceptually simpler and have practical implementations for k-mer sets, such as BBHash [START_REF] Limasset | Fast and scalable minimal perfect hashing for massive key sets[END_REF].

Methods

Set-Min sketch in a nutshell

Assume we are given a set K of keys with associated values taken from a set L with |L| |K|. In our case, K is the set of k-mers occurring in the dataset and L includes their frequencies, although our method will hold for any set of labels L. We want to compactly implement the associative map of (key,value) pairs.

A Set-Min sketch is an R × B matrix M where each bucket is treated as a set, initially empty. Similar to the Count-Min sketch, rows in the matrix correspond to hash functions h i , 0 ≤ i ≤ R -1, that we assume pairwise independent.

At construction time, the key of each (key,value) pair (p, ) is hashed by the hash functions to retrieve its buckets and the value is inserted into each set. Formally, we update

M (i, h i (p)) = M (i, h i (p)) ∪ { } for each row i.
To retrieve the value associated with a key p, we compute the intersection of the corresponding sets, that is ∩ 0≤i≤R-1 h i (p). If the intersection is a singleton, the value is returned. If the intersection is empty, p is not present in the map. If the intersection contains more than one value, we have a collision.

Dealing with collisions

Intersect:

> Final answer: In case of collision, the choice is guided by the number of k-mers supporting each label of the intersection: the label with the smallest support is returned

METHODS

(see Figure 2). The rationale for this is that the label with the smallest support has the smallest probability to appear "by chance", as labels with larger support occur in more buckets and are therefore more likely to occur in the intersection by chance. Thus, the algorithm compares spectrum values for all values of the intersection, and returns the label with the smallest value (ties are broken randomly). If the spectrum is monotonically decreasing (as it is usually the case for large k, see Figure 1), then the label returned is simply the largest one among those in the intersection.

In this work, we assume that only k-mers present in the dataset can be queried. In this case, a query can no longer result in an empty intersection. We further optimize by not storing in the sketch the label 1 with the largest support.

For large k, 1 is usually 1 which is the frequency of the largest fraction of k-mers. 1 is retrieved implicitly: when the intersection is empty, 1 is returned.

This optimization allows us to save space and will be further discussed later.

Note that, with these modification, an error may occur even if the resulting intersection is a singleton but the right label is actually 1 .

We now show that with Set-Min sketch, we can bound the total absolute error over all k-mers of the dataset. Consider a sketch S built on a map assigning to each k-mer p ∈ K a value (label) p ∈ L which is the frequency of p in the dataset. We denote by c the number of k-mers with frequency ∈ L (spectrum value).

Consider

D = p∈K | ˆ p -p | (1)
where ˆ p is the label of p returned by the sketch. Our goal is to dimension R and B such that D ≤ ε a 1 , where a 1 is the total number of k-mers in the dataset (roughly, the dataset size) and 0 < ≤ 1.

Querying p returns an incorrect frequency m = p iff m occurs in the intersection and c m < c p . The probability of this event is

1 -1 - 1 B cm R ≈ 1 -e -cm B R (2)
and the expectation of the error when querying p is then

cm<c p m∈L |m -p | 1 -e -cm B R . (3)
Summing up over all k-mers, we obtain

E[D] = ∈L c cm<c m∈L |m -| 1 -e -cm B R . (4)
The total number of k-mers is

a 1 = ∈L c . Given 0 < ε ≤ 1, our goal is to choose B and R in order to ensure ∈L c cm<c m∈L |m -| 1 -e -cm B R < ε ∈L c .
(5)

Assuming that k is sufficiently large and the spectrum is monotonically decreasing, i.e. c m < c iff m > . ( 5) then rewrites to

∈L c m> m∈L (m -) 1 -e -cm B R < ε ≥1 c . ( 6
)
Assume now that the spectrum follows a power-law with large exponent, that is, c = C • -a for some a > 2. Note that under this assumption, the number of unique k-mers is c 1 = C, and the number of all k-mers is

≥1 c = C • ≥1 1 a-1 ≤ C • a -1 a -2 , since ζ(s) = i≥1 1 i s ≤ s s-1 for s > 1.
We then have the following result.

Theorem 1. Given 0 < ε ≤ 1, if B > C and R, B satisfy R • log B C > log 1 ε , (7)
then ( 6) holds.

Proof. Our goal is to estimate

∈L c m> m∈L (m -) 1 -e -cm B R , (8)
where c = C • -a . We assume B > C and approximate 1-e -cm B ≈ cm B = C B m -a . We further lower-approximate (8) by replacing sums by integrals, thus obtaining

∞ 1 C • -a ∞ (m -) C B m -a R dm d . ( 9
)
Routine computation of the integral yields

C C B R 1 (aR -2)(aR -1)(aR + a + 3) . ( 10
)
The inequality of the Theorem becomes

C B R 1 (aR -2)(aR -1)(aR + a + 3) < ε a -1 a -2 . ( 11
)
The Theorem follows.

The theorem allows us to dimension the Set-Min sketch. For example, one can set B = αC for some constant α > 1 and R = log α 1 ε .

Algorithm 1: Heuristic to compute R and B Data:

{c } ∈L , a 1 = ∈L c , ε Result: R and B R ← 1; B ← 1.44 × c max ; T ← ε a 1 ; E ← E(D) (computed by (4)); while E > T do R ← R + 1; E ← E(D); end M ← R × B; while E < T do R ← R -1; B ← M R ; E ← E(D); end R ← R + 1; B ← M R ;

Computing tighter sketch dimensions

Theorem 1 provides a way to dimension a Set-Min sketch, provided that the spectrum follows a power-law distribution with a sufficiently large parameter a. In order to validate these estimates experimentally, and, at the same time, obtain a tool for computing tighter values B and R for arbitrary spectra, we implemented a simple heuristic hill climbing algorithm to compute those values by directly solving equation 5.

Algorithm 1, given below, starts with R = 1 and some initial value of B and then iteratively increments R and recomputes (4) until equation ( 5) holds true. In the implementation, B is initially set to 1.44 × c max , where c max is the largest spectrum value. After such a value of R is found, the algorithm starts decrementing R while incrementing B to maintain the total space RB constant as long as (5) holds. The final R and B are thus the last which satisfied (5) in the decrementing loop. Note that, for both loops, there is a value of R for which

RESULTS

the exit condition is satisfied. The rationale for this step is to have as small R as possible in order to reduce the query time, while maintaining the total space.

Results

We The other two contain unassembled reads:

USAI Sakai strain at 5x coverage from AFproject, SRR low-coverage human data SRR622461 from the 1000 Genomes Project4

Table 1 summarizes the characteristics of each data set for each value of k in our analysis. Observe that, while the number of distinct k-mers is comparable to the total number of k-mers (data size), the number of distinct k-mer counts is small. This is in accordance with the power-law distribution discussed in Section 2.1.

Set-Min vs Count-Min sketch

Table 2 compares Set-Min sketch to Count-Min sketch. Dimensions R and B

were computed using Algorithm 1 to insure bound (5) to hold for ε = 0.01.

Dimensions of Count-Min sketch were set to be the same. Value 1 is the most common count in all reported datasets and it was not inserted into Count-Min sketch in order to make the comparison againts Set-Min as fair as possible. Zero values are thus interpreted as the non-inserted count.

For ease of comparison, column T reports the threshold ε a 1 given to Algorithm 1. Columns E s and E c report the actual total sum of errors for Set-Min and Count-Min, respectively. In all reported cases E s < T , as expected.

The total error of Count-Min, E c is, most of the time, one order of magnitude larger than E s . For SRR with k = 15, it even exceeds the total number of k-mers a 1 in the dataset.

The average error of Set-Min is, in most cases, very close to 1, which suggests that the overwhelming majority of collisions occur between successive counts such as 1 and 2 -the most abundant ones in the spectra considered here. The average error of Count-Min is bigger but of the same order of magnitude, except for small k and unassembled datasets.

RESULTS

Table 1: Data sheet for the data sets used in our study. Columns T k and D k report the total number (in millions) of k-mers and the number of distinct k-mers, respectively. D c reports the number of distinct k-mer counts. C k reports the number (in millions) of distinct k-mers with a count value different from the most common one (which is 1 in all reported cases). On the other hand, the fraction of k-mers producing an error is in striking contrast: in case of Set-Min, about only 1-3% of distinct k-mers produce an error, while for Count-Min, this fraction is much larger. This shows that Count-Min cannot be used when most of k-mer counts are expected to be retrieved precisely, for comparable sketch sizes.

Type Name T k [M] k D k [M] D c C k [M]
Table 2: Set-Min compared to Count-Min. T is the reference upper bound on the sum of errors equal to ε a 1 (right-hand side of ( 5)). E s and E c are the sum of errors for Set-Min and Count-Min respectively. N s and N c are the percentages (rounded to integers) of distinct k-mers producing an error, for Set-Min and Count-Min, respectively. A s and A c are respective average errors, with average taken over the number of distinct k-mers resulting in an error in the respective sketch.

Name k R B T E s E c N s N c A s A c SAI 11

Set-Min vs Max-Min sketch

We compared Set-Min to Max-Min -an optimized version of Count-Min (see Sect. 2.2). Note that in the case when the k-mer spectrum is strictly decreasing for increasing k-mer counts, the maximum count corresponds to the one with the smallest support. In the general case, we use a variant when, instead of directly comparing count values, counts are ordered according to the support size. That is, updates are performed by keeping the label with the minimum number of k-mers in the k-mer spectrum, and a query returns the label with the smallest such number.

Table 3 compares Set-Min with Max-Min. As expected from theoretical considerations, the performance of Max-Min in terms of the average error and the sum of errors is better than for regular Count-Min, but worse than for Set-Min. The same behaviour is observed for the number of k-mers having an erroneously estimated frequency. Therefore, Max-Min falls in-between Set-Min and Count-Min, providing a simple and inexpensive practical method to enhance the latter, without reaching the accuracy of the former.

Altogether, Table 3 shows that the performance of Max-Min is closer to Count-Min than to Set-Min. This is because, by keeping the maximum element in each bucket, we are reducing each set of Set-Min to a single element opening the possibility of increased collisions by potentially sharing a given maximum element between unrelated k-mers. The intersection operation performed by Set-Min during query is thus strictly necessary, to guarantee the desired error bounds.

Set-Min sketch vs KMC output

Not surprisingly, Set-Min achieves better memory consumptions than KMC in all our tests (columns M kmc and M s of Table 4). Values of R and B do not change from Table 2. The compression rate is variable: from a small factor to two orders of magnitude. The best compression is achieved for larger values of k and assembled genomes. The former is primarily explained by the decreasing number of distinct counts, due to the power-law behaviour. As for the difference between assembled genomes and sequencing data, we will discuss it in more details in Section 5.

Set-Min sketch vs MPHFs

Table 4 also reports the space usage for BBHash to obtain the best memoryoptimized hash functions. Column M bbhash is the space (in bytes) required by the hash function only, while M bball is the space required by the hash function plus the external array of frequencies.

As in the previous case, Set-Min sketch is more memory-efficient when k is large, taking about an order of magnitude less memory than a MPHF. For small values of k, BBHash takes slightly less space and, being exact, may therefore be the preferable choice. However, one should keep in mind that MPHF does not support updates, while a Set-Min sketch is updatable to a certain extent with new (k-mer,count) pairs, and also mergeable with another possibly redundant map.

The behaviour of the unassembled datasets is of particular interest. Even for large k's, MPHF appears to be a better choice for this type of data. The causes of this phenomenon and possible solutions are discussed in Section 5. Set-Min average query time performance is 50% slower than Count-Min (and by extention of Max-Min) and comparable to those of BBHash, when data is very skewed. Following the previous trend, Set-Min sketches appear to be the slowest method for small values of k and for unassembled reads.

Time measurements

S A I k = 1 1 S A I k = 1 5 S A I k = 2 1 M N O k = 1 5 M N O k = 2 1 M N O k = 2 7 M N O k = 3 2 R A I k = 1 5 R A I k = 2 1 R A I k = 2 7 R A I k = 3 2 G R C k = 1 5 G R C k = 2 1 G R C k = 2 7 G R C k = 3 2 U S A I k = 1 1 U S A I k = 1 5 U S A I k = 2 1 S R R k = 1 5 S R R k = 2 1 S R R k = 2 7 S R R k =
S A I k = 1 1 S A I k = 1 5 S A I k = 2 1 M N O k = 1 5 M N O k = 2 1 M N O k = 2 7 M N O k = 3 2 R A I k = 1 5 R A I k = 2 1 R A I k = 2 7 R A I k = 3 2 G R C k = 1 5 G R C k = 2 1 G R C k = 2 7 G R C k = 3 2 U S A I k = 1 1 U S A I k = 1 5 U S A I k = 2 1 S R R k = 1 5 S R R k = 2 1 S R R k = 2 7 S R

Discussion

Unassembled datasets

As seen in Table 4, for the unassembled datasets, Set-Min sketch does not seem to have an advantage in memory usage, even for large k's. We found

)). Frequencies of these "erroneous" k-mers do not have the same statistical behaviour as bona fide k-mers, in particular first spectrum values do not decay at the same rate as the rest of the spectrum. Figure 5 shows the spectrum of the unassembled SRR datasets. One can observe a slower decay behaviour for a few first spectrum values. In this situation, additional rows are needed just to make the sketch able to distinguish, with required precision, between small frequency values. Note that in practice, distinguishing between small frequencies is often irrelevant. For example, many read assemblers simply discard low-frequency reads as a way to de-noise the data. In the case of Set-Min, it is possible to collapse together the first m columns of the spectrum by assigning to all k-mers in this subset the same frequency. This would considerably reduce the sketch size. Formally, error guarantee (5) would not hold anymore, but most of newly introduced errors would be small (typically, equal to 1) and would occur for low counts only.

To check the above, we constructed a Set-Min sketch for SRR with dimensions (R, B) = (4, 3310557), merging together the first five columns of the sorted spectrum and assigning count 5 to all merged k-mers. While the final sum of errors was well above the theoretical limit (10 9 against 7 • 10 6 ), the maximum and average error were respectively 55 and 2.8. In many applications this error level could be acceptable.

Presence-absence information

As introduced in Section 3.2, in this work we assumed that only k-mers present in the dataset can be queried. This assumption allowed us to discard the largest value of the spectrum corresponding to unique k-mers, thereby saving space.

Set-Min sketches can seamlessly work without this assumption, but the space required for storing k-mer counters may not be competitive to other solutions.

An alternative could be to build an additional data structure, such as a Bloom filter, representing presence-absence information for the set of k-mers having the largest count. This allows the discrimination between k-mers absent in the dataset from those present but non-represented in the sketch.

Another scenario occurs when working with multiple datasets of very high similarity, such as a large collection of bacterial strains or a collection of RNA-seq data [START_REF] Solomon | Fast Search of Thousands of Short-Read Sequencing Experiments[END_REF][START_REF] Yu | SeqOthello: querying RNA-seq experiments at scale[END_REF]. In this case, it might be beneficial to build a Bloom filter for the k-mers present in the union of the datasets, and maintain multiple Set-Min sketches to represent k-mer counts in each dataset.

Note that Set-Min sketches can also be helpful for long-term storage and transmission of the k-mer composition of a dataset augmented with count information. The k-mers of the dataset can be reassembled into simplitigs [START_REF] Břinda | Simplitigs as an efficient and scalable representation of de Bruijn graphs[END_REF] with a Set-Min sketch storing the (approximated) frequencies. The full count table can be restored from the simplitigs and the sketch.

Conclusions

We presented Set-Min sketch -a novel sketching method inspired by the Count-Min sketch. Its primary use is to associate keys to labels without explicitly storing the former. In this paper, we demonstrated the performance of Set-Min sketch for storing k-mer counts information, where the distribution of labels (k-mer counts) follows a power-law distribution. Under this assumption, we proposed simple bounds for a Set-Min sketch that guarantee the total error sum to be within an ε fraction of the total number of k-mers in the dataset.

We showed that Set-Min sketch allows us to save space compared to the raw output of the popular KMC k-mer counting tool when applied to labels following a skewed distribution, at the price of a very modest error rate. This saving is especially important in the case of whole-genome data and large values of k, where it can achieve a two orders of magnitude reduction in memory usage.

Set-Min has been shown to be more space efficient than the MPHF-based soluton for large values of k. For smaller k's, however, MPHFs provide an implementation with comparable memory consumption. Finally, when compared to Count-Min and Max-Min sketches of comparable dimensions, our sketch achieves better point-query errors thanks to the distribution-aware dimensioning performed on the k-mer spectrum, and the reliance on already computed count tables as input.

Figure 1 :

 1 Figure 1: k-mer spectrum of the human genome for k = 32 in log-log scale.

Figure 2 :

 2 Figure 2: Example of collision resolution in case of multiple items occurring in the intersection. The brown label is returned because it is more rare compared to the blue one.

  implemented Set-Min in a software tool named fress, available at https:// github.com/yhhshb/fress. The fress pipeline compares Set-Min with Count-Min and MPHF implementations. BBHash (Limasset et al., 2017) (https:// github.com/rizkg/BBHash) is the only external library required, as Count-Min is implemented within fress. Rather than storing a set in each bucket of the sketch, fress only stores an index to an array of involved sets. Note that the current version of fress does not include any complex optimisation, such as multi-threading or bit-packing of the final matrix. Sorted spectra and lists of involved sets are explicitly stored in text format. We tested Set-Min on six data sets of different size and complexity. Four of them are fully assembled genomes: SAI Sakai strain of Escherichia Coli taken from (Yi and Jin, 2013) (NCBI accession number B000007), MNO genome of Drosophila melanogaster from FlyBase 2 , RAI genome of Gossypium Raimondii (Hatje and Kollmar, 2012) downloadable from AFproject(Zielezinski et al., 2019), GRC human reference genome assembly GRCh38 3 .

Figure 4 :

 4 Figure 4: Average query time of Set-Min, Count-Min and BBHash.

Figure 5 :

 5 Figure 5: Spectrum in log-log scale of SRR unassembled data sets for k = 32.

  4 1.04 • 10 6 5.50 • 10 4 5.00 • 10 4 1.39 • 10 6 USAI 11 5 4.58 • 10 5 2.57 • 10 5 1.99 • 10 5 7.44 • 10 7 2.6 99 2.46 24.6 USAI 15 4 4.79 • 10 6 2.49 • 10 5 2.45 • 10 5 8.01 • 10 6 1.9 33 1.31 2.53 USAI 21 5 3.99 • 10 6 2.38 • 10 5 2.00 • 10 5 7.91 • 10 6 1.6 34 1.21 2.35 MNO 15 5 1.79 • 10 7 1.43 • 10 6 1.39 • 10 6 8.35 • 10 6 .69 • 10 6 1.43 • 10 6 1.41 • 10 6 2.26 • 10 7 1.1 12 1.03 1.61 MNO 27 5 3.72 • 10 6 1.43 • 10 6 1.11 • 10 6 2.28 • 10 7 0.9 12 1.01 1.48 MNO 32 5 3.81 • 10 6 1.42 • 10 6 1.11 • 10 6 2.10 • 10 7

			1.8 26 1.15 2.24
	SAI	15 5 2.13 • 10 5 5.50 • 10 4 4.66 • 10 4 2.38 • 10 5	0.9	4 1.01	1.1
	SAI	21 5 1.20 • 10 5 5.50 • 10 4 5.29 • 10 4 4.57 • 10 5	0.9	7 1.05 1.18
			1.4	7	1	1.25
	MNO 21 4 40.9 12 1.01 1.44
	RAI	15 4 8.36 • 10 7 7.27 • 10 6 5.65 • 10 6 1.50 • 10 8	2.1 27 1.08 2.25
	RAI	21 4 8.47 • 10 7 7.27 • 10 6 5.73 • 10 6 4.09 • 10 7	1	6 1.01	1.3
	RAI	27 4 7.14 • 10 7 7.27 • 10 6 6.49 • 10 6 3.56 • 10 7	1.1	5 1.01 1.22
	RAI	32 4 6.38 • 10 7 7.27 • 10 6 6.87 • 10 6 3.15 • 10 7	1.1	4 1.01 1.17
	GRC 15 3 2.26 • 10 8 2.93 • 10 7 2.78 • 10 7 1.36 • 10 9	2.9 52 1.78 4.74
	GRC 21 4 1.42 • 10 8 2.93 • 10 7 2.60 • 10 7 1.65 • 10 8	1.1	6 1.01 1.23
	GRC 27 4 1.07 • 10 8 2.93 • 10 7 2.82 • 10 7 1.91 • 10 8	1.1	6 1.01 1.25
	GRC 32 4 9.84 • 10 7 2.93 • 10 7 2.92 • 10 7 1.85 • 10 8	1.1	6 1.01 1.22
	SRR 15 4 1.17 • 10 8 7.90 • 10 7 4.93 • 10 7 1.34 • 10 10 3.3 96	2.2 20.68
	SRR 21 3 2.39 • 10 9 7.35 • 10 7 7.09 • 10 7 5.63 • 10 8	1.8	9 1.11 1.68
	SRR 27 4 1.78 • 10 9 6.80 • 10 7 4.92 • 10 7 4.58 • 10 8	1.3	8 1.04 1.53
	SRR 32 4 1.72 • 10 9 6.34 • 10 7 4.95 • 10 7 4.01 • 10 8	1.3	7 1.03 1.48

Table 3 :

 3 Set-Min compared to Max-Min sketch. Columns E c , N c , A c are replaced by E m , N m , A m with the same meaning as theirTable 2 counterparts. USAI 11 2.57 • 10 5 1.99 • 10 5 1.24 • 10 7 2.6 63.5 2.46 6.4 USAI 15 2.49 • 10 5 2.45 • 10 5 1.72 • 10 6 1.9 13.3 1.31 1.33 USAI 21 2.38 • 10 5 2.00 • 10 5 1.76 • 10 6 1.6 14.2 1.21 1.24 MNO 15 1.43 • 10 6 1.39 • 10 6 5.74 • 10 6 1.4 5.6 1 1.02 MNO 21 1.43 • 10 6 1.41 • 10 6 1.78 • 10 7 1.1 11.2 1.03 1.31 MNO 27 1.43 • 10 6 1.11 • 10 6 1.78 • 10 7 0.9 12.1 1.01 1.2 MNO 32 1.42 • 10 6 1.11 • 10 6 1.65 • 10 7 0.9 11.3 1.01 1.17 RAI 15 7.27 • 10 6 5.65 • 10 6 5.81 • 10 7 2.1 16.8 1.08 1.38 RAI 21 7.27 • 10 6 5.73 • 10 6 3.05 • 10 7 GRC 15 2.93 • 10 7 2.78 • 10 7 4.07 • 10 8 2.9 27.8 1.78 2.67 GRC 21 2.93 • 10 7 2.60 • 10 7 1.38 • 10 8 1.1 5.5 1.01 1.08 GRC 27 2.93 • 10 7 2.82 • 10 7 1.62 • 10 8

	Name	k	T	E s	E m	N s N m	A s	A m
	SAI	11 5.50 • 10 4 5.00 • 10 4 4.15 • 10 5 1.8 13.5 1.15 1.29
	SAI	15 5.50 • 10 4 4.66 • 10 4 2.13 • 10 5 0.9	4	1.01 1.01
	SAI	21 5.50 • 10 4 5.29 • 10 4 4.02 • 10 5 0.9	7.2 1.05 1.06
						1	5.2 1.01 1.07
	RAI	27 7.27 • 10 6 6.49 • 10 6 2.84 • 10 7 1.1	4.5 1.01 1.05
	RAI	32 7.27 • 10 6 6.87 • 10 6 2.61 • 10 7 1.1	4	1.01 1.03
						1.1	6	1.01 1.09
	GRC 32 2.93 • 10 7 2.92 • 10 7 1.58 • 10 8 1.1	5.7 1.01 1.08
	SRR 15 7.90 • 10 7 4.93 • 10 7 3.31 • 10 9 3.3 61.8 2.2 7.92
	SRR 21 7.35 • 10 7 7.09 • 10 7 2.42 • 10 8 1.8	5.9 1.11 1.13
	SRR 27 6.80 • 10 7 4.92 • 10 7 2.05 • 10 8 1.3	5.2 1.04 1.06
	SRR 32 6.34 • 10 7 4.95 • 10 7 1.89 • 10 8 1.3	4.9 1.03 1.04

Table 4 :

 4 Set-Min ( = 0.01) compared to KMC and BBHash (run with γ = 1). All memory is reported in bytes. Column M kmc , M s , M bball are the memory taken by a fully functional map between k-mers and their frequencies when applying KMC, Set-Min sketch and BBHash, respectively. M bbhash is the memory of the hash function produced by BBHash without the external array of frequencies. SRR 21 2.91 • 10 10 1.61 • 10 10 1.48 • 10 9 1.10 • 10 10 SRR 27 3.73 • 10 10 1.60 • 10 10 1.48 • 10 9 1.08 • 10 10 SRR 32 4.07 • 10 10 1.46 • 10 10 1.57 • 10 9 1.04 • 10 10

	Name	k	M kmc	M s	M bbhash	M bball
	SAI	11 1.21 • 10 7	5.75 • 10 6	9.13 • 10 5 3.00 • 10 6
	SAI	15 3.80 • 10 7	1.20 • 10 6	2.06 • 10 6 5.99 • 10 6
	SAI	21 4.77 • 10 7	6.77 • 10 5	2.03 • 10 6 6.00 • 10 6
	USAI 11 1.54 • 10 7	1.49 • 10 7	1.17 • 10 6 4.61 • 10 6
	USAI 15 6.95 • 10 7	2.87 • 10 7	3.72 • 10 6 1.35 • 10 7
	USAI 21 8.53 • 10 7	3.00 • 10 7	3.91 • 10 6 1.39 • 10 7
	MNO 15 7.08 • 10 8	1.68 • 10 8	3.87 • 10 2.15 • 10 8
	MNO 21 9.80 • 10 8	3.55 • 10 7	4.66 • 10 7 2.45 • 10 8
	MNO 27 1.24 • 10 9	3.75 • 10 7	4.73 • 10 7 2.48 • 10 8
	MNO 32 1.37 • 10 9	3.84 • 10 7	4.90 • 10 7 2.36 • 10 8
	RAI	15 1.76 • 10 9	7.54 • 10 8	9.59 • 10 7 5.98 • 10 8
	RAI	21 4.37 • 10 9	6.78 • 10 8	2.16 • 10 8 1.24 • 10 9
	RAI	27 6.04 • 10 9	5.36 • 10 8	2.37 • 10 8 1.29 • 10 9
	RAI	32 6.96 • 10 9	4.79 • 10 8	2.70 • 10 8 1.38 • 10 9
	GRC 15 3.83 • 10 9	1.70 • 10 9	2.09 • 10 8 1.65 • 10 9
	GRC 21 1.86 • 10 10 1.28 • 10 9	9.32 • 10 8 6.46 • 10 9
	GRC 27 2.48 • 10 10 9.66 • 10 8	9.86 • 10 8 6.57 • 10 9
	GRC 32 2.82 • 10 10 8.38 • 10 8	1.08 • 10 9 6.54 • 10 9
	SRR 15 4.73 • 10 9	2.00 • 10 9	2.59 • 10 8 2.12 • 10 9

Here we assume k to be sufficiently large, typically k > log 4 L, where L is the data size.

http://flybase.org

ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_ GRCh38/seqs_for_alignment_pipelines.ucsc_ids/GCA_000001405.15_GRCh38_no_alt_ analysis_set.fna.gz

ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR622/SRR622461/SRR622461_1.fastq.gz . Only the file SRR622461_1 is used in this study.

Acknowledgement

The authors are grateful to the RECOMB anonymous reviewers for their helpful comments.

https://github.com/yhhshb/fress

Author Disclosure Statement

The authors declare they have no conflicting financial interests.

Funding Information

GK was partially funded by RFBR, project 20-07-00652, and joint RFBR and JSPS project 20-51-50007.