An Improved multi-objective genetic algorithm based on orthogonal design and adaptive clustering pruning strategy

Two important characteristics of multi-objective evolutionary algorithms are distribution and convergency. As a classic multi-objective genetic algorithm, NSGA-II is widely used in multi-objective optimization fields. However, in NSGA-II, the random population initialization and the strategy of population maintenance based on distance cannot maintain the distribution or convergency of the population well. To dispose these two deficiencies, this paper proposes an improved algorithm, OTNSGA-II II, which has a better performance on distribution and convergency. The new algorithm adopts orthogonal experiment, which selects individuals in manner of a new discontinuing non-dominated sorting and crowding distance, to produce the initial population. And a new pruning strategy based on clustering is proposed to self-adaptively prunes individuals with similar features and poor performance in non-dominated sorting and crowding distance, or to individuals are far away from the Pareto Front according to the degree of intra-class aggregation of clustering results. The new pruning strategy makes population to converge to the Pareto Front more easily and maintain the distribution of population. OTNSGA-II and NSGA-II are compared on various types of test functions to verify the improvement of OTNSGA-II in terms of distribution and convergency.

INTRODUCTION
I N real production and life, most of the actual problems require comprehensive optimization of multiple objectives, which are often highly complex and non-linear. Multi-Objective Evolutionary Algorithms (MOEAs) are very suitable for solving such problems and have become one of the mainstream algorithms for solving multi-objective optimization problems [1]. The purpose of the multi-objective evolutionary algorithm is to make the solution set approaches the Pareto ideal frontier of the problem, and the distribution is wide and uniform, which makes the distribution and convergency become the performance index to evaluate algorithms [2].Distribution and convergency are of great significance to solving multi-objective optimization problems. Well distribution can provide more reasonable and effective options for decision makers, and well convergency can solve practical problems more accurately. In recent years, some multi-objective evolutionary algorithms have been proposed on evolutionary computing. Among  them, the most representatives are: Zitzler and Thiele proposed the SPEA [3], Kim et al. proposed on the basis of SPEA2 [4], Srinivas and Deb proposed non-dominated sorting genetic algorithm NSGA [5], others on the basis of NSGA-II [6], and [7] Pareto Envelope-based Selection Algorithm and PESA-II [8] proposed by Corne. NSGA-II is one of the most widely used multi-objective evolutionary algorithms. Its characteristic is the determination of individual fitness values are based on the Pareto dominance relationship and density information between individuals [9]. However, there are shortcomings of such a fitness calculation in improper maintenance of distribution and convergency. Wen Shihua et al. retained some representative individuals to improve NSGA-II according to the distance measurement method [10]. This method only considers the influence of crowding distance on the distribution of population, and does not take the perspective of the convergency and distribution problems of same kind of individuals caused by the existence of similar individuals and inferior individuals. And the influence of the distribution of initial populations on subsequent evolutionary processes is not considered. Aim at these two deficiencies of NSGA-II, this paper proposes an improved convergency and distribution retention strategy. Before the evolution, the multi-objective orthogonal experiment was used to initialize the population by using the non-dominated sorting and crowding distance evaluation to prevent the population from falling into local convergence or slow convergence due to random initialization, and to avoid the unevenness of individuals inside the initial population which leads to the lack of distribution. In the evolutionary process, the algorithm clusters each evolution, dynamically adjust the strength of the pruning within the class according to the similarity, and maintain the population by adaptively removing the appropriate number of in-class and non-dominated sorting and intraclass individuals with poor crowding distance and a small number of individuals far from the frontier face, speeding up the population convergence while maintaining the distribution. The comparison experiments between OTNSGA-II and NSGA-II show that the algorithm has good distribution and convergency.

RELATED CONCEPTS
The following describes the definition of multi-objective optimization problems from a rigorous mathematics perspective. The mathematical definition of the multi-objective optimization (MOP) problem that is commonly adopted and generally accepted in the field of multi-objective optimization is as follows: Definition 1 (MOP): A general MOP consists of n decision variable parameters, k objective functions, and m constraints. The objective function, constraints, and decision variables are functional relationships. The optimization objectives are as follows: e(x) = (e 1 (x), e 2 (x), ..., e m (x)) ≤ 0 among x = (x 1 , x 2 , ..., x n ) ∈ X y = (y 1 , y 2 , ..., y k ) ∈ Y Where x represents the decision vector, y represents the target vector, X represents the decision space formed by the decision vector x, Y represents the target space formed by the target vector y, and the constraint e(x) ≤ 0 determines the feasible range of values of the decision vector.
When there are multiple objective functions, the concept of "optimal solution" produces new changes. When solving a multi-objective problem, it is actually seeking a set of equilibrium solutions instead of a single global optimal solution. The concept of this widely adopted optimal solution was proposed by Francis Ysidro Edgeworth early in 1881. Subsequently, the famous French economist and sociologist Vilfredo Pareto promoted the concept in 1896. From the perspective of economics, he transformed multiple incomparable objects into single indicators for optimal solution, involving the concept of multi-objectives. Pareto first proposed the concept of vector optimization, which is Pareto optimal that is used widely nowadays.
The difficulty of MOP is that in most cases, the targets may be conflicting and the improvement of one target may cause the performance of other targets to be reduced. It is generally impossible to achieve multiple targets at the same time. Otherwise, it is not multi-objective optimization. The ultimate solution to the MOP can only be compromising and balancing among the objectives, so that the objective functions can be optimized as far as possible. Therefore, the optimal solution of MOP is essentially different from the optimal solution of single-objective optimization problem. In order to solve MOP correctly, the concept of solution must be defined. Definition 2 (feasible solution set): The feasible solution set X f is defined as a set of decision vectors x satisfying the constraint condition e(x) in equation 1 The expression of the target space corresponding to the feasible region of X f is: Equation 3 represents that all xs in the feasible solution set Xf form a subspace in the target space through the optimization function mapping, and the decision vectors of the subspace belong to the feasible solution sets. For the minimization problem, it can be easily converted to the above maximization problem to solve.
The feasible solution set of the single-objective optimization problem can determine the merit relationship and degree of the scheme through the unique objective function f. For the MOP problem, the situation is different, because in general, the decision vector in X f cannot be sorted completely, but only a certain index can be sorted, that is, partially sorted [11]. To this end, define the mathematical relationshipbetween decision vectors =, >, ≥ as follows: Definition 4: For vector u, v: Similarly, the following mathematical relations can be defined: <,≤.
Then define the mathematical relationship between the decision vectors: Definition 5 (Pareto win): For decision vectors a, b: If the vector a, b has a Pareto superior relationship, the values of all the objective function of one vector are greater than the values corresponding to the other vector. Now define the concept of the Pareto optimal solution: Definition 6 (Pareto optimal solution): The decision vector x ∈ X f is called the Pareto optimal solution if and only if: The Pareto optimal solution is also known as non-inferior solutions or Non-dominated solutions. If the decision vector x is a Pareto optimal solution, this means that x cannot improve any of the objective functions without weakening at least one other objective function. For a Pareto optimal solution in a given target space Y, its corresponding point in the decision space is called efficient or non-inferior.
One point in X f is efficient if and only if its image is non-dominant in Y. Different from the Single-objective Optimization Problem, there is usually no unique optimal solution in MOP problem, but a Pareto optimal solution set. The concept of Pareto optimal solution set is as follows: Definition 7 (Pareto optimal solution set): For a given MOP problem, the Pareto optimal solution set A is defined as: In most cases, in the multi-objective optimization problem, only the Pareto optimal solution exists [12], and the optimal solution similar to the single-objective optimization does not exist. Multi-objective optimization problems, mostly with many Pareto optimal solutions, its Pareto optimal solution is only an acceptable non-inferior solution or satisfactory solution. If there are so-called optimal solutions for a multi-objective optimization problem, they must be a Pareto optimal solution consisting of only these optimal solutions and no other solutions. Therefore, the Pareto optimal solution is A reasonable solution set for multiobjective optimization problems [13] [14]. For the practical application problem, one or more solutions must be selected from the Pareto optimal solution of the multi-objective optimization problem as the optimal solution for the multiobjective optimization problem according to the degree of understanding of the problem and the personal preference of the decision-makers [15] [16]. Therefore, the first step and key to solving the multi-objective optimization problem is to find as many Pareto optimal solutions as possible. Two goals for multi-objective optimization are as follows: (1) Find a set of solutions that are as close as possible to the Pareto optimal domain.
(2) Find a set of solutions that are as different as possible and evenly distributed.
The first goal is necessary in any optimization work. It is not advisable to converge to a set of solutions that are not close to the true Pareto optimal solution set. Only when the solution converges to approximate true optimal solutions can we guarantee their approximate optimal properties. In addition to converge to the approximate optimal domain, the solutions obtained must also be distributed as evenly as possible throughout the Pareto optimal domain. Only by finding a set of solutions that are sufficiently diverse and covering as much as possible can each of them include a number of different effective solutions for decision makers to choose.

The problem of NSGA-II population initialization
In the NSGA-II method, the initial population is formed by using the random initialization method, that is, the value of the solution is randomly selected within the range of the values of the independent variables of each dimension, thereby forming a new individual, which will result in the distribution of randomly generated individual sets is not uniform, and most of the individuals are not close to the frontier, making the early evolutionary exploration easy to fall into the local optimum, which is not conducive to converge to the Pareto frontier, and also causes the lack of distribution. The randomly initialized population is shown in Figure 1: In Figure 1, we can see that the initial population distribution generated by this random initialization method is not uniform enough, some individuals are too concentrated, some individuals are too separated; and most individuals are far from Pareto ideal frontier, which is not conducive to the maintenance of convergency and distribution in the subsequent evolution process inevitably.

Method for setting initial population in fault multiobjective
Orthogonal experimental design, as a method that can solve multi-factor and multi-level experimental problems effectively [17], uses orthogonal tables to arrange a few experiments to find the best or better experimental conditions [18] [19] [20], so it is widely used in solving the problem of optimization [20]. Suppose an experimental system has F factors, each of each with Q levels. If a comprehensive combination test is performed, a experiments are required. While using orthogonal table L M (Q F ) , just choose M combinations to test. Where L M (Q F ) represents an orthogonal table with F factors and Q levels, L represents the Latin square, M represents the number of horizontal combinations, there are M rows, each row represents a horizontal combination, and M is generally much smaller than [21]. For the four-factor and three-level problems, if all the horizontal combinations are fully tested, a total of 3 4 = 81 tests are required; and if an orthogonal table is used, such as formula 8, only 9 tests are required, which can greatly reduce the number of experiments.
For convenience of representation, it is noted that L M (Q F ) = [a i,j ] M ×F , where a ij represents the horizontal value of the j th factor of the ith combination , if j th satisfies j = 1, 2, (Q 2 − 1) ÷ (Q − 1) + 1, (Q 3 − 1) ÷ (Q − 1) + 1, ...(Q j − 1) ÷ (Q − 1) + 1 , the j th column is called the basic column, other columns are non-basic columns. Orthogonal table creation process is as follows [22] [23]: In order to make the population converge quickly to the Pareto frontier and distribute evenly, the initial population should be covered as much as possible to cover the entire feasible domain and close to the Pareto ideal frontier in the population initialization process of the multi-objective optimization problem. The multi-objective orthogonal experiment of fault is a multi-factor and multi-level mathe-Algorithm 1 Construction of Orthogonal Array Algorithm 1: Step 1: Calculate the minimum value J that satisfies 2: 9: Step 3: Construct the basic columns: 10: for k = 1 to J do do 11: j = Q k−1 Q−1 + 1 12: for i = 1 to i < Q J do do 13: end for 15: end for 16: Step 4:Construct a non-basic column: 17: for k = 2 to J do do 18: for s = 1 to j − 1 do do 20: for t = 1 to Q − 1 do do 21: a j+(s−1)/(q−1)+t = (a s × t + a j )mod Q 22: end for 23: end for 24: end for 25: Step5 matical experiment method from the perspective of Pareto multi-objective optimization. According to the orthogonal table, some representative points with uniform dispersion and uniformity are selected from the comprehensive experiment to carry out experiments. At the same time, for multi-objective problems, it is necessary to calculate the non-dominated sorting layer according to the relationship between multiple targets, and select a part of the individuals in the front layer as an alternative set. And in this alternative set, the method of evaluating the individual based on the non-dominated sorting stratification and the crowding distance is selected to better adapt to the initial population scene in which the multiple targets are mutually restrained while maintaining the convergence and distribution of the population. . The good distribution and convergency of the initial population ensure the diversity of the population and the richer model and close to the optimal solution set, so that the algorithm can converge at a faster speed in the global scope and has a uniform distribution and a wide coverage. When using the fault multi-objective orthogonal experiment method to initialize the population, if the feasible solution space [l, u] is large, in order to improve the search efficiency and precision, the feasible solution space [l, u] is divided into S subspaces (For the segmentation method, see Algorithm 2 in this section), then construct the orthogonal table L M (Q F ). According to the orthogonal table, for each subspace, the proposed SOC operator is used for cross operation, and a discretized initial individual set is generated, so that the candidate set is selected therefrom and the initial population is composed by selection. 1: Step 1: Select the s th dimension of the feasible solution space, where s satisfies the following formula:

Algorithm 2 Subspace Segmentation
The OTNSGA-II design initializes the population by multi-objective orthogonal experiments by using faults and non-dominated sorting and crowding distances for individual selection. Using this method of tomographic multiobjective orthogonal experiments, the population is initialized before evolution, so that the initialized individuals can be evenly distributed at different positions close to the Pareto front, preventing the population from falling into local convergence or convergence is too slow due to random initialization, and avoiding the uneven distribution of the initial population individuals leads to the loss of distribution and thus is not conducive to extending to the entire Pareto front as much as possible. The population initialized by the fault multi-objective orthogonal experiment is shown in Figure 2.
In Figure 2, we can see that the initial population of the fault multi-objective orthogonal experiment is evenly distributed, and most individuals are close to the Pareto ideal frontier, which not only maintains good initial distribution and convergence, but also lays the foundation for the continued maintenance of the above two excellent characteristics in the subsequent evolution process. From this algorithm flow that maintains the convergency and distribution of the initial population, it can be seen that: In the process of designing the fault multi-objective orthogonal experiment to initialize the population, according to the multi-factor and multi-level characteristics, the distributed Algorithm 3 self-adaptive orthogonal crossover (SOC) 1: Step 1:Let P 1 = (P 1,1 , P 1,2 , ..., P 1,N ), P 2 = (P 1,1 , P 1,2 , ..., P 1,N ) be the two parent individuals involved in the cross operation. The feasible solution space determined by P 1 and P 2 is [lparent, uparent].
where set J is the set of positions where the components of low similarity in the parental individuals p 1 and p 2 are located, and t is the number of components of low similarity in p 1 and p 2 , and θ0 is a given positive real number close to zero. The vector k holds the position of the component with low similarity in P 1 , P 2 , that is, the position at which the factor is divided. Let x be any individual in 3: Let k 0 = 0, then the Q levels of the i th factor f i can be expressed as: 4: Step 3: According to the algorithm 1, construct the Based on the orthogonal table L M (Q F ), orthogonal experimental design is performed for the t factors determined in equation (12) and the Q levels corresponding to each factor determined in equation (13), and M child generation individuals will be generated: population with uniform distribution and neatness is generated, which ensures excellent initial distribution. Selecting the individuals with higher non-dominated stratification in the discretized population into the candidate set maintains the characteristics of these individuals close to the Pareto ideal frontier and improves the initial convergency performance of the population. In the alternative set, according to the non-dominated sorting stratification and the crowding distance calculation, the better individuals are selected to form the initial population, so that the maintenance of Algorithm 4 Orthogonal initialization population 1: Step 1: Divide the feasible solution space [L, U ] in each dimension of the optimization problem into S subspaces, and the segmentation method is shown in Algorithm 2 2: Step 2: Take the horizontal number of the orthogonal table as Q 0 , and cross-operate each subspace based on the SOC operator (algorithm 3) to generate a new discretized population P , where Q 0 is the discretized horizontal number; 3: Step 3: Calculate the value fit of the target of each individual in the population P in the problem to be optimized, expressed as the value of the i th individual on the i th objective function; 4: Step 4:Perform calculation of non-dominated sorting stratification according to the superiority and inferior relationship between each target value fit for each individual in the population P , and mark the level c i of each individual, which is represented as the level at which the i th individual is located; 5: Step 5:According to the hierarchical ranking, the individual layers are sequentially acquired from the first layer into another set until the number of individuals in the set is greater than or equal to 4 * n.( If the individual number of population P is less than 4 * n, take the individual collection of the population.) and the number of different layers in the set is m. This set is an alternative set, where n is the number of individuals in the population; 6: Step 6:Calculate the congestion distance d i as the value of the i th individual congestion distance for each individual in the candidate set by using the candidate set as the ensemble space; 7: Step 7:In the alternative set, the non-dominated sorting level c i and the crowding distance d i of each individual are taken as two targets, and the non-dominated relationship between the individuals is evaluated to select the optimal first n individuals to generate the initial population P0.
the population distribution is comprehensively considered while ensuring the proximity to the ideal frontier. However, the way in which the algorithm NSGA-II randomly initializes the population leads to uneven distribution of individuals and most individuals are far away from the ideal frontier, which is not conducive to the convergency and distribution.

Problems in the calculation of crowding distance in NSGA-II
NSGA-II ranks each individual hierarchically by the non-dominant relationship of each target between individuals and takes the number of layers as the quantitative fitness value of the individual, and calculates the crowding distance between individuals in the population. Then each individual is ranked according to the number of layers in which the individual is located and the calculated crowding distance [24], and selects a certain number of individuals from the front to form the next generation of new populations. The crowding distance of each individual is the sum of the distance differences between the individual and the two individuals adjacent to each other in each target dimension [25] [26]. As shown in Fig. 4-1, two targets f1 and f2 are provided, and the crowding distance of the individual i is the sum of the length and the width of the dotted quadrilateral in the figure. Let d[i]distance be the crowding distance of the individual i, and d [i].m is the function value of the individual i on the target m, then the crowding distance of the individual i in Figure 3 is: In general, when there are r targets, the crowding distance of individual i is: The method of calculating the crowding distance of an individual is as shown in Algorithm 5. However, using the above crowding distance calculation method to select individuals may result in some individuals with good distribution to be eliminated, while others with poor distribution may be retained, which will lead to the concentration of individuals in some areas, and the sparseness of others in other areas [27] [28], which is not conducive to the maintenance of population distribution; and a small number of individuals away from the frontier may also be retained, which will reduce the possibility that the next generation of new populations will continue to evolve to produce new superior individuals, which is not conducive to populations. As shown in Figure 4, for a generation of solution sets in the evolution process, individuals a and b have similar characteristics, showing that the two positions are close and far from other points. Using NSGA-II method to calculate the crowding degree of individuals a and b, the values are relatively close and relatively large, indicating that individuals a and b are in the same set of similar features, then the points in this collection are likely to be eliminated or retained at the same time, resulting in individuals in the region being too sparse or too dense, and the ideal situation is to retain a better performing individual in the set as a representative point with the feature set. The same is true for individuals c and d. Individuals e, f , and g are more evenly distributed, but the crowding distance obtained by the NSGA-II method will be small, and it is likely to be eliminated during the selection process, resulting in the loss of individuals in this area and the ideal situation is to retain one or two excellent representative points to maintain the distribution characteristics of similar collections of this feature. It can be seen that the calculation method of the crowding distance will affect the retention of the population distribution. As shown in Figure 5, according to the calculation method of NSGA-II, individual a will enter the next gen-eration for subsequent evolution. However, for the characteristic individual a, unlike other individuals who are far from the Pareto front, the individual a tends to produce some poorly performing individuals in the next generation of evolution, which is detrimental to the rapid convergence of the population.It can be seen that this calculation method also affects the convergency performance of the population.

Fig. 6: Schematic diagram of a generation of evolutionary solution sets
A concrete example is further explained in Figure 6: in this figure, because the rankings of individuals a and b in the non-dominated sorting are ranked last and away from the Pareto frontier, its existence is not conducive to the generation of good individuals in the subsequent evolution process, resulting in the entire population converge too slowly; although individuals c, d, e, and f are in the same sorting layer, because their characteristics are the same and the crowding distance is too low, they belong to the same set of individuals with similar characteristics, which may result in being eliminated or retained at the same time. If it is eliminated at the same time, the point set of this segment is empty, which is not conducive to the expansion to the entire frontier of Pareto; If too many individuals with similar characteristics are retained at the same time, it is not conducive to the combination of new individuals in the process of evolution, which leads to the excessive distribution of these individuals, which hinders the distribution maintenance of the solution set.

Improved pruning strategy
In the method of the algorithm OTNSGA-II, the individuals in the population are first clustered, and the appropriate number of non-dominated sorting and individuals with poor crowding distance are adaptively trimmed according to the similarity within the class. In this way, it is ensured that some points with poor ranking and congestion distance are removed from each set with similar features, and most of the non-dominated sorting points and the representative points with better crowding distance are retained, so that the point sets with the same characteristics will not be the same which is not conducive to further expansion of the solution set, and it is easier to combine excellent individuals; filtering a few point sets far away from the front surface is beneficial to the next generation of evolution to quickly converge to the Pareto front. In this way, the distribution and convergency of the population are maintained from the above two aspects. There are three basic issues that need to be addressed in this way: (1) the calculation of the average similarity within the class; (2) the setting of the number of individuals within the class; (3) the choice of the pruned individual. The calculation of the average similarity within the class: This paper uses K-means to cluster the population of each iteration result, calculate the similarity for any two bodies in the class, and then calculate the average of all two body similarities as the intra-class average similarity, the formula is as follows: Where i, j represent an individual, M is the number of individuals within the class, and p ij represents the similarity between two individuals i, j. Set the number of individuals in the reserved class: 1 − δ * the average similarity within each class, and the average similarity within the class is 0 1 (In practice, it is usually 0 0.8, because there are few cases where the average similarity within the class is 1, that is, all the individuals in the class repeat.) In the case of guaranteeing the main flow of the algorithm, adaptively trimming a small number of inferior individuals based on the similarity within the class, retaining most of the representative points with the similar feature set , the formula is as follows: p k refers to the intra-class average similarity in the above formula 17. Due to the characteristics of the solution set distribution of different types of functions and the different dimensions of each target, a large number of experiments have found that the appropriate adaptive pruning parameters δ are 0.12-0.15. The parameter should not be too large or too small. If it is too large, it will trim out too many individuals so that some excellent individuals can be trimmed together; if too small, the effect of removing inferior individuals with higher similarity is not obvious. The method of selecting the pruned individuals: in the same class, retaining the above number of individuals based on a comparison rule that considers the non-dominated sorting firstly and considers the crowding distance secondly, that is, the selection is based on the non-dominated sorting hierarchical ranking and the value of the crowding distance. For two Individuals, if the non-dominated sorting hierarchy is different, individuals with larger levels are selected; if the non-dominated sorting level is the same, individuals with larger crowding distances are selected.
For example, as shown in Figure 7, for a certain generation of the solution set in the evolution process, taking the coordinate values in each dimension of the individual and the non-dominated hierarchical ranking and the crowding distance as features, K-means clustering is performed on the solution set and is divided into five categories (areas in the five graphs A, B, C, D, and E). Through the calculation of the average similarity within the class, some individuals with the same non-dominated hierarchical ranking and poor crowding distance in the same category are filtered out. The higher the intra-class similarity, the more individuals with similar features need to be removed. The average similarity within a class is proportional to the number of individuals removed within the class, inversely proportional to the number of individuals retained within the class. Taking different intra-class average similarities as a metric, the appropriate numbers of individuals are removed from each class to achieve self-adaptive effect. Specifically, in the A class, because the hierarchical ranking of the individual a in the non-dominated sorting is the last and away from the Pareto frontier, its existence is not conducive to the combination of new individuals with characteristics that maintain the characteristics of population convergency in the subsequent evolution process. On the contrary, it increases the probability of deviating from the Pareto frontier individuals, resulting in the convergence of the entire population is too slow, the average similarity in this set is relatively small through the adaptive clustering pruning strategy, only the individual a away from the Pareto front surface is removed; in the D class, the individual b has the same situation as the individual a, and is removed; in the C class, there are 5 individuals in the class, and the individual g has a relatively large crowding distance, while for the four individuals c, d, e, and f , they are in the same nondominated sorting hierarchy, the features are very similar and the crowding distance is too low, which may result in being eliminated or retained at the same time. If it is eliminated at the same time, the point set in this area will be empty, which is not conducive to the expansion of the entire frontier of Pareto; if it is retained at the same time, it will make the similar individuals in this area too much and the distribution is too dense, which is not conducive to the emergence of new and better individuals with more maintainable population distribution characteristics during the evolution process, which hinders the distribution of the population. The adaptive clustering pruning strategy calculates that the average similarity in this set is relatively large, and the two individuals with poor performance, c and d, are removed together (the crowding distance of the individual g is relatively large and belongs to the outstanding individual and will be retained); due to the low intra-class average similarity, the individual collections in class B and class E are all retained by calculation and are not pruned.
In summary, the adaptive clustering pruning strategy adaptively removes some non-dominated sorting stratification and individuals with poor crowding distance or away from the ideal frontier for the average similarity of individual sets of similar features. Therefore, the convergency and distribution of the population are well maintained. Step 5: Substitute the similarity p ij between two bodies in the class into the formula (17), and calculate the average similarity P k in each set; 6: Step6. Substitute the average similarity P k in each set into the formula (18), and calculate the number n k of individuals in each class that need to be retained; 7: Step 7: According to the calculated number of individuals n k retained by each class, trim the (n U k −n k ) individuals with poor non-dominated hierarchical ranking and crowding distance and the individuals away from the ideal frontier, where n U k is the number of individuals in the k th set; 8: Step 8: Reorganize the n k better individuals retained after these pruning into the new population P new and continue the next generation of evolution.
It can be seen from the algorithm flow that in the process of population maintenance, the population distribution is adaptively maintained based on the aggregation characteristics of each set in the population. According to the average similarity of individual sets with similar features, an adaptively pruned a small number of individuals with similar features and non-dominated hierarchical ranking and poor crowding distance, retaining most representative individuals. At the same time, a small number of individuals far from the frontier surface were filtered, which accelerated the rapid convergence of the population to the Pareto frontier during the evolution process. However, the NSGA-II crowding distance solving method is used to complete the selection operation, so that there are a large number of individuals with similar features in the iterative process, in which individuals in each set divided by similar features have a high probability of being homogenous which is not conducive to the maintenance of population distribution; And retaining individual individuals away from the Pareto front, is not conducive to the maintenance of population convergency.

Experimental environment and test functions
In this experiment, the 26 test functions of SCH, FON, POL, KUR, ZDT series, DTLZ series and UF series of multi-objective optimization problems are used to test the algorithm effect. The program runs under Windows 7 environment of 2.6 GHz CPU and 8G memory.These multiobjective test functions with different characteristics such as continuity, discontinuity, convexity, and concavity are easy to test the optimization effect of the algorithm under different types of problems. Parameter setting: In NSGA-II and OTNSGA-II (in this paper algorithm), population size popsize = 100, running algebra ngen = 250, both algorithms use real number coding, crossover probability P m = 0.9, mutation probability P c = 0.1, To ensure that the comparison experiment was carried out under the same conditions.

Performance evaluation
The two basic evaluation metrics for multi-objective optimization are convergency and distribution. The calculation of the solution set evaluation in multi-objective optimization is an important and complex problem [29], and the researchers have proposed some effective methods [?]. For the convergency evaluation, Generational Disdance (GD) [33] can be used to measure the degree of convergence between the final solution set of the algorithm iteration and the Pareto ideal frontier. The formula is as follows: In equation 19, n is the number of individuals in the solution set, and d i is the minimum of the Euclidean distance from each individual to each point in the Pareto ideal frontier. The smaller the value of GD, the closer the solution set is to the Pareto ideal frontier. If GD = 0, the solution of the algorithm is on the Pareto ideal frontier, which is also the most ideal experimental result. For the distribution evaluation, the SP [34] proposed by Schott can be used to measure the uniformity of the final solution set distribution of the algorithm iteration. The calculation formula is as follows: Where d i represents the distance between adjacent points parallel to the Pareto ideal front, di = min(|f 1i (x) − f 1j (x)|+|f 2i (x)−f 2j (x)|), (i, j = 1, 2, ..., n) ,d is the average of all d i . SP = 0 when the final solution set of the algorithm iteration is completely evenly distributed in the target space. In order to comprehensively evaluate the convergency and distribution of the final solution set of the algorithm iteration, the inverse generation distance (inverted generational distance, referred to as IGD) [35] can be used as a performance evaluation index. IGD is a measure of the distance between the true Pareto frontier and the final solution set of the algorithm iteration. The smaller the value calculated by the index value, the better the convergency and distribution performance of the final solution set of the algorithm iteration, and the closer to the true Pareto frontier and evenly distributed. Let P be the real Pareto optimal solution set, and A is the final solution set of the evolutionary algorithm iteration, then the formula is as follows: where,Dist i = min

Experiment and data analysis
In the aforementioned experimental environment, using the two different algorithms of NSGA-II and OTNSGA-II (the algorithm of this paper), comparative experiments were performed on the 26 different types of test functions. The experimental results are shown in Table 1, Table 2, Table 3: According to the above experimental results, from the comparison of the convergency indicators in Table 1, we can find that the convergency of the 26 functions of SCH, FON, POL, KUR, ZDT series, DTLZ series, UF series has been improved, even the 14 functions of SCH, KUR, ZDT1, ZDT2, ZDT3, ZDT6, DTLZ1, DTLZ3, UF1, UF2, UF3, UF5, UF6, UF10 are increased by one order of magnitude or two orders of magnitude. This is because OTNSGA-II uses the method of fault multi-objective orthogonal experiment to initialize the population and select the individual set close to the Pareto ideal front from the evenly distributed individuals as the initial population. This initial advantage of more agile characteristics accelerates the evolutionary trend of subsequent Pareto ideal frontiers; and in the process of evolution, the evolutionary results of each generation are clustered according to the principle of similar individual characteristics within the same class, and adaptively trimming individuals that are farther away from the Pareto ideal frontier in each type of individual set, so that there are more opportunities for the remaining good individuals to combine to produce better individuals, moving closer to Pareto ideals. The direction of the frontier evolved.
We can find that in the comparison of the distribution indexes in Table 2  series and UF series has been improved, Even on the five functions of ZDT3, ZDT6, DTLZ3, UF1, UF3, it is increased by one order of magnitude or two orders of magnitude. It is only slightly worse on the four functions of FON, POL, DTLZ2 and UF9. This is because OTNSGA-II uses the method of fault multi-objective orthogonal experiment to initialize the population, so that the population is evenly distributed in the spatial range closer to the Pareto ideal frontier, maintaining a good initial distribution for the subsequent evolution of the population and laiding the foundation for the population to continue to maintain this good distribution during subsequent evolution; And in the process of evolution, the appropriate number of intra-class non-dominated sorting and individuals with poor crowding distance are trimmed by clustering adaptively, so that most of the individual sets with similar features retain most of the well-represented representative points, which will not make the individual collections of similar features co-exist, thus maintaining a better distribution of the population. However, because the four functions of FON, POL, DTLZ2, and UF9 are characterized by the fact that the Pareto ideal frontier is rather lengthy, In the process of initializing the population by the multi-objective orthogonal experiment of the fault, the method of discretizing the individual by nondominated sorting and crowding distance is easy to lose some representative points to the maintenance of distribu- tion, and in the adaptive clustering pruning strategy, it is easy to cut off some representative points too much, which will cause a slight loss of distribution, resulting in a slightly worse performance on these four functions.
We can see that in the comparison of the comprehensive performance indicators in Table 3,the convergence of the 26 functions of SCH, FON, POL, KUR, ZDT series, DTLZ series and UF series has been improved, even in FON, ZDT2, ZDT3. The 11 functions of ZDT6, DTLZ1, UF1, UF2, UF3, UF5, UF7, and UF10 are increased by one order of magnitude or two orders of magnitude. This is because OTNSGA-II uses the method of fault multi-objective or-thogonal experiment to initialize the population and select the individual set close to the Pareto ideal front from the uniformly distributed individuals as the initial population, which not only ensures a good initial distribution of the population, but also has the initial convergence of individual individuals in the population close to the Pareto ideal frontier. And adaptively trimming the appropriate number of non-dominated sorts and individuals with poor crowding distance and individuals far from the Pareto ideal frontier in each type of individual set, which ensures that most of the individual collections with similar characteristics retain most of the better representative points, and do not survive the same life, maintaining the distribution of the population; but also filters out a small number of inferior individuals with large deviations, making it easier for the retained good individuals to combine new and better individuals, so that the population accelerates to converge on the Pareto ideal frontier and maintains the convergency of the population.
Taking the function UF3 as an example, the detailed changes of the GD, SP, and IGD evaluation index values of the two algorithm populations with algebra are observed, as shown in Figure 8, 9, 10. In Figure 8, we can see from the comparison of the GD values representing the initial population (0 th generation) that the OTNSGA-II design fault multi-objective orthogonal experiment initializes the population so that the generated initial population is close to the ideal frontier. Therefore, the initial convergency of the population is maintained, which is conducive to maintaining this good convergency advantage and accelerating the population convergency in the subsequent evolution process; and the GD value remains stable after 20 generations, which is due to due to the use of adaptive clustering pruning strategy in the evolution process, so that some individuals deviating from the optimal solution set are eliminated, prompting the population to converge rapidly to the Pareto ideal frontier. Before the evolution of NSGA-II (0 th generation), the GD value has fallen behind, and it has remained basically stable after 200 years of evolution and is still in a backward state. This is due to the random initialization of the population, which makes the initial abundance of the population not perform well, and because of the retention of individuals away from the frontiers during evolution, the constant participation in the evolution and other individual genes slows down the rapid convergence of the entire population towards the Pareto ideal frontier. In Figure 9, it can be seen from the comparison of the SP values representing the initial population (0 th generation) that the initial population of the OTNSGA-II fault multiobjective orthogonal experiment can make the initial population distribution uniform and maintain the initial distribution of the population, which is conducive to the subsequent evolution to continue to maintain this good distribution; And it can be seen from the figure that there is a period of high SP value between the 20 th and 40 th generations of the pre-evolution, which is due to the fact that crossvariation between individuals makes the population highly dense. Using the adaptive clustering pruning strategy, the individual collections with similar features will not be the same, leaving most of the representative points in the same set, and removing a few inferior points with similar features and poor performance. The SP value decreased significantly after the 40 th generation, which is due to the fact that the use of this pruning strategy makes the distribution of the population significantly effective. However, before the evolution of NSGA-II (0th generation), the SP value has lagged behind, and it has gradually fallen behind and become more serious in the evolution process. This is due to the random distribution, which makes the initial distribution poor, and with the continuous evolution of the evolution process, especially after the 50 th generation, individuals with similar characteristics appear in large numbers and are continuously combined into new individuals with similar characteristics, resulting in high SP values for evaluation distribution, which is not conducive to the maintenance of population distribution. In Figure 10, we can see from the comparison of IGD val-ues representing the initial population (0 th generation) that the OTNSGA-II fault multi-objective orthogonal experiment initializes the population to make the initial population distribution uniform and close to the Pareto ideal frontier, and maintains good convergency; And the IGD value remained stable after 20 generations. This is due to the use of the adaptive clustering pruning strategy, which adaptively trims out the non-dominated sorting stratification and the poor crowding distance of individuals with similar characteristics in the same individual set and retains most of the outstanding representative points, so that individuals in similar collections do not simultaneously present or simultaneously pruned, and a small number of individuals away from the front surface are removed, so that the retained excellent individuals produce better individual combinations, thus maintaining the distribution of the population while maintaining the convergency of the population. On the contrary, the IGD value before the evolution of NSGA-II (0 th generation) has fallen behind, and it remained stable after the 220 th generation and remained in a backward state. This is because random initialization makes the initial population distribution uneven and most individuals are far away from the frontier surface, so that they do not have good initial distribution and convergence. In the process of evolution, a large number of individuals are combined with each other to make the population too dense in some spaces, which leads to the simultaneous retention or simultaneous elimination of individuals with similar characteristics, and the existence of some individuals far from the frontier is not conducive to the combination of excellent individuals. So, in the evolution process, it can not maintain good distribution and convergency.

CONCLUSIONS
By analyzing the random population initialization method and the inter-individual crowding distance calculation method in NSGA-II, aiming at the problem that the population distribution and convergency of the algorithm are not properly maintained, this paper proposes an improved algorithm NSNSGA-II of NSGA-II, which uses the method of non-dominated sorting and crowding distance evaluation to set up multi-objective orthogonal experiment to initialize the population. And the iterative results of each generation are clustered by K-means, and some inferior individuals in the class are pruned, and the convergency and distribution of the evolutionary population are maintained. The use of fault multi-objective orthogonal initialization in OTNSGA-II makes the initial individual distribution in the population uniform and close to the Pareto ideal front, maintaining good initial distribution and initial convergence; In the evolution process, the adaptive clustering pruning strategy is used to adaptively prune the appropriate number of individuals with poor crowding distance and non-dominated sorting and individuals far from the Pareto ideal frontier within the same class according to the degree of aggregation within the class, and Prevent individual collections of similar features from being eliminated or retained at the same time, and make the excellent genes easy to retain and develop in the subsequent evolution process, speed up the understanding of the convergence of the set to


extfieldsofstudy:
Computer Science
Medicine

s2fieldsofstudy:
Computer Science
Biology

text:
Detection of False-Positive Deletions from the Database of Genomic Variants

Next generation sequencing is an emerging technology that has been widely used in the detection of genomic variants. However, since its depth of coverage, a main signature used for variant calling, is affected greatly by biases such as GC content and mappability, some callings are false positives. In this study, we utilized paired-end read mapping, another signature that is not affected by the aforementioned biases, to detect false-positive deletions in the database of genomic variants. We first identified 1923 suspicious variants that may be false positives and then conducted validation studies on each suspicious variant, which detected 583 false-positive deletions. Finally we analysed the distribution of these false positives by chromosome, sample, and size. Hopefully, incorrect documentation and annotations in downstream studies can be avoided by correcting these false positives in public repositories.

Introduction
A genomic variant is an alteration of the DNA sequence of an organism. Since an organism's DNA sequence encodes the genetic instructions used in its development, any alteration of this sequence may cause genetic abnormalities or even fatality. According to their sizes, genomic variants are classified into small-scale variants, such as single nucleotide polymorphisms (SNP) and indels, and large-scale variants, namely, structural variations (SV), including copy number variations (CNV), insertions, deletions, inversions, segmental duplications, and translocations [1]. Various complex diseases have been reported to be associated with genomic variants in human genomes [2].
Prior to next generation sequencing (NGS), cytogenetic techniques, such as fluorescence in situ hybridization (FISH) and array comparative genomic hybridization (aCGH), were employed to detect SV. However, due to their relatively low genomic resolution (e.g., microscopic scale (Mbp) for FISH, and submicroscopic scale (kbp) for aCGH) [3], most medical and biological research teams have migrated their platforms to NGS, which can provide base-pair level resolution.
Several approaches have been proposed to detect SV from NGS data. Generally, these approaches can be classified into two categories: paired-end read mapping (PEM)-based or depth of coverage (DOC)-based approaches [4]. For PEMbased approaches, if the span of a pair of mapped reads is longer/shorter than a specified cutoff related to the insert size of the sequencing library, a deletion/insertion can be identified [5], whereas for DOC-based approaches, if the local depth of reads is significantly larger/smaller than the global DOC, a duplication/deletion can be identified [3,6,7]. PEMbased approaches have advantage of detecting balanced SVs (inversion) and unbalanced SVs (deletion and insertion) of relatively small sizes, whereas DOC-based approaches are good at detecting unbalanced SVs (CNV) of relatively large sizes. Besides PEM and DOC, several other supplementary signatures, such as split read mapping, have been combined into PEM and DOC to improve detection performance, leading to integrative models [8][9][10].
Despite PEM's advantage in balanced SV detection, the majority of detection approaches use DOC as the primary signature to identify CNVs [11]. However, the DOC signature is biased due to two main factors: GC content [12] and mappability [13]. BioMed Research International (i) Since G and C form a triple hydrogen bond (whereas A and T form a double bound), theoretically the melting temperature of GC-rich segments is approximately 2 ∘ C higher than that of AT-rich segments [14]. As a result, when the sequencing protocol involves polymerase chain reaction (PCR), GC-rich segments and AT-rich ones are unevenly amplified [15], yielding the correlation between DOC and GC content [16,17].
(ii) Due to the complexity of the human genome, there are regions in which sequenced short reads cannot be uniquely mapped, e.g., repeated regions, such as retrotransposons (LINE and SINE) [13]. Mappability was introduced to measure the uniqueness of such regions using a score that ranges from 0 to 1, corresponding to highly repeated and unique regions, respectively. From the definition, it is clear that the DOC is correlated with mappability.
In human genomes, both GC content and mappability are distributed unevenly along chromosomes, and, therefore, they introduce biases into DOC. Several methods have been developed to correct these two biases [13,18,19].
However, due to the overlook of biases introduced by GC content and mappability, some DOC-based SV studies contain false detection. From the mechanism of PEM and DOC signatures, it is clear that PEM is less affected by GC content or mappability than DOC, and, therefore, PEM can be used to detect false positives. In this paper, we used this idea to verify the entries in the database of genomic variants (DGV). We hope that incorrect documentation and annotations can be avoided in downstream studies by correcting the false positives in this public repository and other related ones such as EMBL-EBI's Genomic Variants archive (DGVa) and NCBI's dbVar.

Samples and Data.
The database of genomic variants (DGV, http://dgv.tcag.ca/dgv/app/home) [20] provides a comprehensive summary of structural variation (SV) in the human genome. In DGV, SVs are defined as genomic alterations that involve segments of DNA with length larger than 50 bp. In DGV, the 6.4 millions of variants represent collections from 55 thousand healthy control samples in 72 studies. DGV provides a curated catalogue of genomic variations in the human genome, which was integrated into EMBL-EBI's Genomic Variants archive (DGVa) and NCBI's dbVar. Therefore, this database is of tremendous importance to investigators whose study interest is about genomic variance, which is also the focus of the current study.
The 1000 Genomes Project (http://www.1000genomes .org/ [21]) is a well-known international collaborative NGS project, which aims to sequence the genomes of approximately 2500 people from 25 populations around the world. In our study, the BGZF compressed sequence alignment/map (SAM) data files (BAM) of most samples were downloaded from the website of this project as the primary dataset.
Sequence read archive (SRA, https://www.ncbi.nlm.nih .gov/sra) is the NCBI database that stores raw sequencing data obtained from NGS technology. The aim of SRA is to make NGS data available to researchers to both improve reproducibility and enable new discoveries. This database includes data from most common sequencing platforms and most NGS studies. The BAM file of some specific samples in our study was not available, so the SRA files were downloaded as the primary dataset.

Methods.
The steps conducted in the study are shown as follows, and the pseudocode is listed in Pseudocode 1 to illustrate the logical structure of these steps.
(1) The latest spreadsheet of supporting variants was downloaded from the DGV website. The GC content and mappability profile of each chromosome of hg18 were downloaded from the readDepth website (https://github.com/chrisamiller/readdepth) [18], which is an NGS-based CNV detection software package.
(2) GC content profiles were smoothed with LOESS [22], and segments with size larger than 500 bp and an average GC content lower than ℎ 1 = 0.26 or greater than ℎ 2 = 0.59 were obtained as suspicious regions.
(3) Mappability profiles were smoothed with LOESS, and segments with size larger than 500 bp and an average mappability lower than ℎ 3 = 0.92 were added to suspicious regions.
(4) All supporting variants in DGV were resolved oneby-one to collect the fields needed in the current study, including variant accession ID, chromosome, genomic location (starting and ending loci), variant subtype, reference, method, and samples used.
(5) Variants associated with the sequencing method and loss or deletion subtype, size smaller than 10 kbp, and nonempty sample fields were filtered for further analysis.
(7) GDV variants that overlapped (F-score greater than 0.9) with suspicious regions were identified as suspicious variants. Here are some notes that should be addressed: (i) A supporting variant represents a variant called in a single sample/individual, which can also be described as sample level variant [20].
(ii) One sample is the minimal requirement to verify a specific supporting variant, so the sample field should be nonempty.
(iii) The terms deletion and loss are equivalent in the database [20].
(iv) Duplicated variants are defined in the sense of the same chromosome ID and genomic location.
(v) The corresponding BAM or SRA file of a suspicious variant was retrieved according to the sample and reference field.

Validation.
The validation of variants is based on the PEM signature. First, for each suspicious variant, we extracted all the read pairs in which both ends were mapped within the region of interest (ROI) from the corresponding BAM file. To provide an adaptive zoom, the ROI is defined as the genomic region that extends both upstream and downstream with 1 kbp plus half of the variant length. Next, the F-score [23] is employed to measure the overlapping quality between the span of a suspicious variant and that of a mapped read pair. The F-scores quantify the overlap quality between two spans, with values ranging from 0 to 1 (see Figure 1, which demonstrates several typical scores). A small value close to 0 means a bad overlap, whereas a high value close to 1 means a good overlap. The F-score is calculated as follows: for a test span, if it has no overlap with the reference span, the F-score is set to 0; otherwise, = 2( /( + )), where is the precision (percentage of the test span that overlaps with the reference span) and is the recall (percentage of the reference span that overlaps with the test span).
In our study, mapped read pairs with F-scores larger than 0.7 were selected, and the average and sum of the mapping quality of all selected pairs were calculated. A suspicious variant was identified as a true positive if the average and sum were above 30 and 90, respectively; otherwise, it was classified

Results
We identified a total of 1923 suspicious variants, which cluster in 7 samples from 8 studies (see Figure 4(b) and   Figure 2 shows two typical suspicious variants from NA18507. Figure 2(a) shows the variant with accession essv4528478, whose genomic location is chr1:210537566-210539442 (green bar in all panels). It is shown that both GC content (the magenta curve in the upper left panel) and mappability (the red curve in the middle left panel) profiles at the ROI are abnormal, and hence they yield a valley in the DOC profile (the blue curve in the lower left panel). As a result, a deletion variant was detected in the DOC profile. However, the PEM signature (the right panel) contains no mapped read pairs (the horizontal lines) that overlap with the green bar, suggesting that this suspicious variant is a false positive. In contrast, Figure 2(b) shows the variant with accession essv4968609 (genomic location chr1:154793140-154795767), whose GC content, mappability, and DOC profiles show similar behaviours to those in (a). However, the PEM signature contains well-mapped read pairs (three black lines in the centre) that overlap with the green bar, suggesting a true positive. Figure 3(a) shows the distribution of the GC content of the human genome (hg18) with a bin size of 100 bp. Based on this distribution, the threshold values ℎ 1 = 0.26 and ℎ 2 = 0.59 are used to determine the estimated extreme GC content regions, such that both the left and right tail areas cover 5% of the whole distribution. Figure 3(b) shows the distribution of the mappability of the human genome (hg18) with a bin size of 100 bp. Since exact 1 and 0 mappability values occupy a large portion of the distribution (68% and 7%, respectively), these two values were excluded from the distribution. Based on the remaining values, the threshold value ℎ 3 = 0.92 was chosen such that the right area covers 20% of the distribution.  Size of SVs / log 10 Figure 4(c) shows the distribution of suspicious variants and false positives with respect to the size. Three modalities are shown: the left one with size smaller than 1 kbp (3), the middle one with size between 1 kbp and 5 kbp (3.7), and the right one with size larger than 5 kbp, which have proportions of 36%, 46%, and 18%, respectively. These results indicate that most false positives are small-or medium-sized variants. By fitting each modality with a Gaussian curve, the means of the three modalities are 660 bp (2.8), 2.2 kbp (3.3), and 7.2 kbp (3.9) for the left, middle, and right modalities, respectively.

Conclusion and Discussion
We proposed an approach to detect GC content and mappability related to false positives from the database of genomic variants. The proposed approach utilized the PEM signature, whose presence is necessary and provides evidence for true positives. 583 false positives were detected by conducting a validation study on the database of genomic variant. The results can avoid incorrect documentation and annotations in downstream studies.
We excluded variants with sizes larger than 10 kbp in this study, and the reasons are as follows: for most NGS alignment/mapping tools, the maximal insert size is limited to thousands of base pairs; e.g., the default values of both thea parameter in BWA and -X parameter in Bowtie 2 are 500 bp, and a very large insert size degrades mapping performance. However, large variants (deletions and inversions) do require a large insert size. As a result, there is a conflict between the mapping quality and the maximal size of detectable variants, and a tradeoff has to be taken with caution. Therefore, we confined the maximal variant size to 10 kbp.
In this study, since there are several software packages/algorithms (e.g., smoothing, segmentation, mapping, etc.) and parameters (F-score, thresholds ℎ 1 , ℎ 2 , ℎ 3 , etc.) that are used in the method, the robustness is an important issue. A global optimization of all parameters is not conducted due to the larger number of parameters, but we set each software/algorithm to its recommended setting and tune each parameter separately to a reasonable value (e.g., th1, th2, and th3). We used the F-score to identify whether a segment overlaps with another segment. From Figure 1, we can see that a threshold value of 0.7 is appropriate to determine that two segments are almost overlapping with each other, so we used 0.7 in the validation step to determine an overlap status. When it is increased to 0.75, 644 false positives are detected, and when it is decreased to 0.65, 546 false positives are detected. Therefore, the results are roughly robust with respect to this parameter. In step 7, we used a large threshold value of 0.9 in order to narrow down the total number of suspicious variants to be validated. We also used three threshold values ℎ ( = 1, 2, 3) to identify suspicious regions. In Figure 3(a), we chose ℎ 1 = 0.26 and ℎ 2 = 0.59 so that both the left and right tail areas covered 5% (or 10% in total) of the whole distribution, whereas in (b) we chose the threshold value ℎ 3 = 0.92 so that the right area covered 20% of the distribution. Since the distribution of GC content is close to a Gaussian distribution, tail areas with 10% are appropriate. However, the distribution of mappability is far from a Gaussian distribution, and values other than 1 are unfavourable. When we adopted the strategy used for GC content to mappability, i.e., the left tail area covering 10%, the resultant threshold value was ℎ 3 , which was too tight to identify suspicious variants. Therefore, we chose 20% for the right area, which yielded a much looser threshold.
There are two limitations of the current study. First, since the PEM signature is more straightforward for studying deletion than other types of variants, the current study focuses on only this type of variant. In future studies, we hope to extend the spectrum of variants being studied. Second, this study focused on only the database of genomic variants as a pilot study to validate our method. In future works, we hope to conduct large-scale validations on other well-known public repositories related to structural variants, such as dbVar (https://www.ncbi.nlm.nih.gov/dbvar), which is NCBI's database of human genomic structural variation, and the Database of Genomic Variants archive (DGVa, https://www.ebi.ac.uk/dgva), which is an EMBL-EBI database that archives publicly available genomic structural variants of all species.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.