Application of the MAHDS Method for Multiple Alignment of Highly Diverged Amino Acid Sequences

The aim of this work was to compare the multiple alignment methods MAHDS, T-Coffee, MUSCLE, Clustal Omega, Kalign, MAFFT, and PRANK in their ability to align highly divergent amino acid sequences. To accomplish this, we created test amino acid sequences with an average number of substitutions per amino acid (x) from 0.6 to 5.6, a total of 81 sets. Comparison of the performance of sequence alignments constructed by MAHDS and previously developed algorithms using the CS and Z score criteria and the benchmark alignment database (BAliBASE) indicated that, although the quality of the alignments built with MAHDS was somewhat lower than that of the other algorithms, it was compensated by greater statistical significance. MAHDS could construct statistically significant alignments of artificial sequences with x ≤ 4.8, whereas the other algorithms (T-Coffee, MUSCLE, Clustal Omega, Kalign, MAFFT, and PRANK) could not perform that at x > 2.4. The application of MAHDS to align 21 families of highly diverged proteins (identity < 20%) from Pfam and HOMSTRAD databases showed that it could calculate statistically significant alignments in cases when the other methods failed. Thus, MAHDS could be used to construct statistically significant multiple alignments of highly divergent protein sequences, which accumulated multiple mutations during evolution.


Introduction
Algorithms for multiple sequence alignment (MSA) are widely used in modern molecular biology for numerous purposes, including evolutionary studies and molecular modeling. Thus, the publication describing the ClustalW multiple alignment algorithm [1] is among the top 10 most cited publications [2]. MSA is indispensable when it is necessary to determine the relatedness of polypeptide sequences containing substitutions and/or insertions/deletions (indels) of amino acid residues, and MSA algorithms are fundamental tools for the construction of phylogenetic trees, calculation of profiles for hidden Markov models (HMMs), and search for common motifs. An important application of MSA methods is the prediction of secondary and tertiary structures of polypeptides and RNA, as well as their molecular functions and interactions [3,4].
In general, the construction of an MSA is an NP-complete problem, for which no direct solution by dynamic programming methods is possible in a reasonable time, even using modern supercomputers [5]. Therefore, many heuristic methods for performing MSA in a relatively short computational time have been developed. Among them, one of the most popular is the progressive procedure based on the construction of a binary (guiding) tree in which sequences are presented as leaves [6]. Then, MSA can be built through a series of pairwise alignments added according to the branching order in the path tree. Methods that employ the progressive approach include ClustalW [1] and its variant Clustal Omega, based on HMM for protein alignments [7], MAFFT [8], and T-Coffee [9]; the latter uses a paired alignment library, which allows for the construction of more accurate multiple alignments than those produced by ClustalW or Clustal Omega.
Iterative methods such as PRRN/PRRP [10] and MUSCLE [11] have been developed with the aim to reduce errors in constructing progressive alignments, which is achieved by repeatedly rebuilding the generated MSAs. Iterative methods revert to pairwise alignments, which can optimize the objective function and improve the quality of MSA.
Some MSA algorithms such as SAM [12] and HMMER [13] use HMM. Compared to progressive methods, HMM-based methods often calculate more significant alignments, which is achieved by rebuilding the alignment each time a new sequence is added.
There are also other methods for calculating MSA, such as genetic algorithms [14], annealing modeling [15], and phylogeny-aware gap placement [16,17]. Recently, we have developed an algorithm for multiple alignments of highly diverged sequences (MAHDS), which is based on the optimization of random position weight matrix (PWM) images of MSAs [18][19][20]. To optimize random patterns in MSA calculation, the MAHDS method uses a genetic algorithm and two-dimensional dynamic programming rather than pairwise equalization; as a result, multiple alignments for highly divergent nucleotide sequences can be obtained. The application of MAHDS to random nucleotide sequences shows that it can calculate statistically significant multiple alignments of sequences, for which the number of random base substitutions per nucleotide (x) ranges from 0 to 4.4. At the same time, the other methods such as ClustalW [21], Clustal Omega [7], MAFFT [22], T-Coffee [9], Muscle [11], and Kalign [23] could do it for x from 0 to 2.4 [18], indicating that for highly divergent nucleotide sequences (x = 2. 4-4.4), the MAHDS algorithm is the only method capable of building statistically significant alignments. Previously, MAHDS has been used to construct MSAs of highly diverged promoter sequences (x = 3.6) from different genomes [19,24,25], which made it possible to identify a large number of potential promoters in rice genes [25].
In this study, we used the MAHDS method for the calculation of MSAs of protein sequences. The application of MAHDS to the alignment of random amino acid sequences with different x showed that the method could produce statistically significant MSAs for sequences with x ranging from 0 to 3.6. The alignments constructed by MAHDS and the other methods mentioned above were performed using the column score (CS) test [26][27][28], Z-score, and the benchmark alignment database (BAliBASE) (reference set 10) containing 218 protein families. Although, according to the CS criterion, the quality of the alignments built by the MAHDS method was somewhat lower than that of the other algorithms, it was compensated by greater statistical significance. Furthermore, in the alignment of 21 highly diverged protein families from the Pfam and HOMSTRAD databases, MAHDS could calculate statistically significant MSAs of more families than the other methods. Our results demonstrate that the developed MAHDS algorithm can be applied to analyze the phylogenetic relationships of evolutionary distant proteins.

Results
The performance of the MAHDS algorithm was compared with that of the MSA methods provided by EMBL-EBI resources, including T-Coffee, MUSCLE, PRANK, Clustal Omega, Kalign, and MAFFT, for which the web application programming interface is available. The following test data were used for analysis: benchmark alignment database (BAliBASE) [27][28][29], artificial sequences with certain properties, and protein families with low average sequence identity in full alignment.
The quality of alignments was evaluated using the CS [26][27][28] and Z-score criteria. In the calculation of the CS, the columns of the alignment are considered. If all remainders in the i-th column are aligned as in the reference alignment, then C i = 1; otherwise, C i = 0. The formula for calculating the CS is: where L is the number of columns in the evaluated alignment. As in order to use the CS, it is necessary to have reference alignments, the CS was applied only to evaluate the BAliBASE using its tools bali score and bali_score_reliable; the former considers all sequences, whereas the latter ignores sequences with discrepancies (annotated with SE-QERR features [30]). The MAHDS algorithm uses the following parameters: gap opening penalty (d), gap extension penalty (e), R L , and K d . R L is the multiplier parameter, which can be used to scale R 2 (Formula (2)), and K d is the equivalent of an expected E score value features [31], which defines the accuracy of determining the start and end of the local alignment (Formula (3)). Therefore, the primary task in testing MAHDS was to find a biologically relevant combination of d, e, R L , and K d .

Testing MAHDS Performance Using the BAliBASE
We used BAliBASE reference set 10, comprising 218 protein families. However, because of restraints set on the file size and the number of sequences by EMBL-EBI services, it was necessary to exclude five families (BBA0039.fasta, BBA0101.fasta, BBA0134.fasta, BBA0190.fasta, and BBA0213.fasta) for which alignment could not be performed by at least one of the compared methods. As a result, we considered the remaining 213 protein families and calculated the average CS for the MSAs as an indicator of the biological correctness, which should be maximized when choosing the parameters. At the same time, we tried to ensure that Z (Methods, Formula (7)) was close to the maximum.
We carried out a series of experiments for different values of parameters d, e, and K d at fixed R L = 5.0, which can be chosen arbitrarily, as d, e, and K d would then be adjusted accordingly. Initially, we considered e = 0.25 d as the optimal ratio, which was previously established for the application of MAHDS to DNA sequences [18], and estimated Z and CS scores for the sets of d values {8, 10,12,16,20,24,28,32,36,40,44,48, 52, 56} and K d values {−0.1, −1, −2}. However, it appeared that with e = 0.25 d there was a discrepancy between the peaks in the CS and Z; therefore, an additional series of experiments was conducted with varying e-to-d ratios (Table 1). The results indicated that the combination d = 40, e = 1, and K d = −1.0 was the optimal for protein sequence alignment and Z estimation with MAHDS. Then, we aligned the 213 protein families by the other methods and estimated the average Z (calculated as described in Section 4.6.1) and CS ( Table 2). The results showed that the MAHDS algorithm was superior to the compared methods in terms of statistical significance (Z) but inferior in terms of CS. As a rule of thumb, paired alignment is considered statistically significant at Z ≥ 8 [32]. For MAHDS, the threshold value was 6 [18] or 10 [33], depending on the nature of the sequences; however, it has been determined for the alphabet of four nucleotides and may not be relevant for protein sequences composed of 20 amino acids.
To determine the threshold of statistical significance Z t , we created random sequences. To ensure that the number of aligned sequences and their average length did not affect Z, two datasets were generated: H 1 containing 100 subsets of 100 random sequences, each 600 amino acids long, and H 2 containing 100 subsets, including 20, 40, and 40 random sequences of 600, 900, and 300 amino acids, respectively. For each H 1 and H 2 subset, we built sequence alignments using the MAHDS algorithm and calculated Z. For set H 1 , Z H 1 = 4.147 and D H 1 = 2.11, whereas for set H 2 , Z H 2 = 4.04 and D H 2 = 1.78. These results showed that the variations in sequence lengths did not affect Z, indicating that the same Z t could be used for sequences of different lengths. Therefore, we chose Z t = 10.0, which was approximately equal to the average value of Z H 1 + 3 D H 1 .

Testing MAHDS Performance Using Artificial Sequences
To compare methods for constructing multiple alignments, we created 81 sets, each containing 100 artificial sequences Des(i) generated from ancestor sequence Anc of length (L) = 600 symbols (Methods, Section 4.7). To create Des(i), the following random changes were made to Anc: the numbers of insertions or deletions were 2, 5, or 10, the indel length was 1, 5, or 20, and the number of substitutions per amino acid was 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, or 2.7. In accordance with Formula (11), the average number of substitutions per symbol between sequences in sets Des was 0.6, 1.2, 1.8, 2.4, 3.0, 3.6, 4.2, 4.8, or 5.4.
We aligned each of the 81 artificial Des sequence sets using MAHDS, T-Coffee, MUS-CLE, Clustal Omega, Kalign, MAFFT, and PRANK and calculated Z values (Tables 3-6; only x values for which Z > 0.0 are shown).  The results indicated that MAHDS could produce statistically significant alignments even if the evolutionary divergence between the aligned sequences (x) was 4.8 substitutions per amino acid; however, it applied to sequences with relatively small numbers of short indels. In other cases, MAHDS can build multiple alignments for x ≤ 3.6 ( Table 3). At the same time, most of the other tested methods could not produce statistically significant results for sequences with x > 1.2 (Tables 4-6). MUSCLE showed the best performance among the other tested methods, being able to build statistically significant alignments for sequences with x ≤ 2.4 and small numbers of short indels; however, with five indels of 20 residues each, the x value dropped to ≤ 1.2 (Table 5).

Testing MAHDS Performance on Protein Families with Low Identity
As shown above, MAHDS could construct statistically significant MSAs for artificial sequences with x ≤ 4.8, whereas the other tested algorithms could not produce statistically significant alignments at x > 2.4. In the next experiment, we verified the efficiency of MAHDS in aligning divergent sequences with a high x value using protein families from Pfam and HOMSTRAD databases with a low percentage of identity.
In total, 21 families with an identity of less than 20% were aligned using the MAHDS algorithm, and the constructed MSAs were compared with those built using T-Coffee and MUSCLE chosen because they could produce statistically significant alignments for artificial sequences with the largest x among the other tested methods. These 21 families are shown in Table 7, and the names of protein families received from Pfam are shown in Table 8.  The performance of MAHDS, T-Coffee, and MUSCLE was evaluated based on statistical significance (Z), the number of gap openings, and the total number of gaps in alignments. The results shown in Table 7 indicated that, for most analyzed families, MAHDS could build alignments of higher statistical significance than the other two methods, with the exception of PF02950, PF09624, PF09987, and PF19443 families, for which the MSAs constructed by MUSCLE had greater Z (however, it should be noted that for PF02950, both MAHDS and MUSCLE showed Z < 10.0). At the same time, only MAHDS was able to construct statistically significant alignments for PF00915, PF10846, PF10895, and PF13944 families. Furthermore, there was a tendency for MAHDS to use fewer gap openings, although the series of gaps tended to be longer than for the other two methods.

Discussion
Here, we described the application of our MAHDS method, which was previously used for the alignment of DNA sequences [18][19][20], to the construction of statistically significant alignments of highly divergent protein sequences that have accumulated multiple amino acid substitutions and/or indels. Comparison with the currently available MSA algorithms revealed that MAHDS could produce significant alignments of model sequences with low identity (x up to 4.8), whereas the other methods (T-Coffee, MUSCLE, Clustal Omega, Kalign, MAFFT, and PRANK) could only carry this out for less diverged sequences (x < 2.4). The superior performance of the MAHDS method is due to the use of the progressive algorithm instead of pairwise alignment and calculation of MSA based on patterns of random alignments. In this case, it is possible to replace N-dimensional dynamic programming, which requires significant computational resources, with two-dimensional dynamic programming. Such an approach enabled the construction of more statistically significant MSAs for 21 protein families with a low percentage of identity than could be performed by the other methods.
We made an assessment of the speed of the MAHDS. For 100 amino acid sequences with a length of 100, 400, and 700 amino acids, the running time of the algorithm was 58, 851, and 2370 s, respectively. This suggests that the calculation time grows approximately in a quadratic dependence on the length of the sequences. For sets of amino acid sequences containing 100, 400, and 700 sequences having a length of 100 amino acids, the multiple alignment building times were 58, 333, and 629 s. These results show that there is an approximately linear relationship between the time required to build an alignment and the number of sequences. For calculations, two computing nodes were used. Each computing node includes two Intel Xeon Gold 6240 CPUs. Each processor has nine physical cores, i.e., we used 36 cores.
The developed algorithm made it possible to calculate the statistical significance of any alignment and compare the main methods used in terms of the statistical significance of the alignments that they can construct for protein families. Statistical significance in our work is estimated by the Monte Carlo method based on the generation of random multiple alignments. Using such a Z-score allows us to make direct probabilistic estimates based on which we can filter out statistically insignificant multiple alignments. Such probabilistic estimates are difficult to make using currently used object functions. Such functions include sum of pairs, minimum entropy, or NorMD [34,35].
A question arises regarding when the developed MAHDS method for MSA of amino acid sequences should be applied and what results could be expected. Obviously, MAHDS would be most useful for analysis of those amino acid sequences that have accumulated multiple mutations (x > 2.4), thus facilitating the prediction of the biological role for the amino acid sequences of unknown function [36]. For this, MSA of the protein family in question could be calculated using MAHDS and then analyzed by similarity search using the profile method [37] or HMM [38]. However, the quality of the MSA can strongly influence the search for statistically significant similarities [39]. Therefore, it is important to use MAHDS to compare highly divergent sequences and calculate MSA. Probably, it could be possible to identify residues, for which the biological function is either unknown or the reliability of the annotation is questionable.
The MAHDS method could also be used to reconstruct the evolutionary history of sequences. In this respect, MAHDS can provide an advantage by constructing statistically significant alignments and tracing similarity between evolutionary distant sequences with x = 2.4-4.8, thus resolving the problem of MSA uncertainty [40] and revealing previously unknown evolutionary relationships between different protein families.
By disclosing similarities among highly diverged proteins, MAHDS can also be applied to predict various parameters of protein conformation, including fold recognition and structure-based functional activity [41][42][43], which can help in the construction of threedimensional models for various proteins [44].
We have developed a site, http://victoria.biengi.ac.ru/mahds/main (3 June 2020), where the user can build multiple alignments for any set of amino acid sequences using MAHDS. Our work was performed using resources of the NRNU MEPhI high-performance computing center.

Algorithm to Calculate Multiple Alignment
The MAHDS method was developed early for the multi alignment of promoter sequences and DNA sequences with weak similarity [18][19][20]. This method uses mathematical approaches and programs that were previously created to search for tandem repeats in various sequences [33,[45][46][47][48][49]. In the present work, the site http://victoria.biengi.ac.ru/ mahds/main (3 June 2020) was used to construct a multiple alignment by the MAHDS method for various amino acid sequences.
Let us briefly consider the mathematical algorithm that underlies the MAHDS method. This algorithm is described in detail in [18]. Let us calculate multiple alignment for a set of N sequences denoted as SI. The optimal multiple alignment (MA) for this set is characterized by a maximum of some function Ψ (maxΨ), which is calculated from MA. Then, we define the image of multiple alignment (IMA) as a function of MA: IMA = f(MA) and also consider an inverse function to f : f −1 (IMA) = MA, which allows us to construct MA from IMA and the SI sequences. Specifying functions f and f −1 makes it possible to unambiguously convert MA to its image (IMA) and vice versa, i.e., to establish one-to-one correspondence between MA and IMA.
Several existing algorithms, including progressive alignment [6], HMM [7], and others [3,4] can produce alignment MA', which is close to the optimal MA. In these methods, the construction of multiple alignment is defined by function ζ: MA' = ζ(SI). The main disadvantage of such an approach is that the alignment is based on pairwise comparisons of sequences from set SI, which precludes building of statistically significant MA' when the number of nucleotide substitutions x is greater than 2.4 [18]. A statistically significant MA' has a probability P(Ψ > Ψ 0 ) < α, where Ψ 0 is the value of function Ψ for MA'. This probability is calculated by aligning a large number of sets SR, each of which comprising a series of randomly shuffled sequences of set SI and determining the distribution of function Ψ for this alignment. The value of α can be chosen as 0.05, if the alignment is built once.
However, there is another way of constructing MA', which is not based on pairwise alignments and direct calculation of MA' = ζ(SI) but uses patterns (images) of random multiple alignments (IMAR) [18]. These images can be subjected to optimization using genetic algorithms in order to create from a random image IMAR such an image IMA m that would be the closest to IMA = f(MA). The degree of closeness between IMA m and IMA cannot be assessed, since MA and its image IMA = f(MA) are unknown at x > 2.4 for most protein families. However, we can arrive as close as possible to IMA by increasing the similarity measure between each sequence from set SI and IMA m. , which can be carried out using global two-dimensional alignment. Each such comparison would produce a maximum value of similarity function Fmax(i) (i = 1, 2, ..., N), and the sum mF = ∑ N i=1 Fmax(i) would be used as the measure of similarity between IMA m and IMA. The greater mF is obtained with the same global alignment parameters, the more similarity there is between IMA m and each sequence from set SI and the closer IMA m is to IMA. We will also use mF as an objective function when running a genetic algorithm to calculate IMA m , which then would be used to define multiple alignment: MA m = f −1 (IMA m ).
In this approach to calculate multiple alignment, we only need to define functions f and f −1 . As function f, we can take the algorithm for creating a position weight matrix (PWM), which would be the image of multiple alignment, i.e., IMA = PWM and IMA m = PWM m (described below in Section 4.2). As function f −1 , we take the algorithm for the global alignment of PWM m and each sequence from set SI (described earlier [50] and in Section 4.3 and Section 4.4). As a result, it is possible to build multiple alignment with good accuracy using PWM m . Since the sequences in set SI can have different lengths, we calculate the average length L of the sequence in the set and create a PWM in the range from 0.9L to 1.1L.
Compared to the method described in [18], here we did not combine all sequences from set SI into one large sequence to calculate global alignment with the PWM but instead aligned the PWM with each SI sequence separately, which significantly improved the quality of MSA. Furthermore, the alphabet of SI sequences (A in ) was expanded from 4 to 20 symbols to correspond to the number of amino acids so that protein sequences could be aligned. The scheme of the algorithm is shown in Figure 1.  (4)). MA m is a multiple alignment built for sequences from SI set using PWM m (see Section 4.5 and Formula (5)).
At the beginning, we have SI sequences whose number is N (step 1). In step 2, we created a set (Q) of random PWMs, each of which served as a random image IMAR. In step 3, we used a genetic algorithm to optimize each PWM from set Q and determine PWM m = IMA m . In step 4, we created multiple alignment MA m for sequences from set SI, which corresponded to the found PWM m . Finally, in step 5 we evaluated the statistical significance of the Monte Carlo alignment.

Creation of Set Q of Random PWMs
To create set Q of random PWMs, PWM rows were represented by 20 amino acids and columns by integers from 1 to L; thus, the PWM dimensions were 20 rows and L columns. A random PWM was obtained from a random amino acid sequence S 1 of length Lk (k = 10 3 ), in which the frequencies of amino acids corresponded with those in the sequences of set SI. Then, a sequence containing integers from 1 to L was generated, copied k times, and joined in tandem to yield sequence S 2 .
Then, we performed transformation of the matrices from set Q, which was necessary to ensure that the distribution function of different Fmax(i) (i = 1, 2, ..., N) was similar in the alignment of sequences from set SI with each random PWM. For normalization purposes, two restrictions were imposed on PWMs: where p 1 (i) and p 2 (i) are probabilities of amino acids in sequence S 1 and of symbols in sequence S 2. , respectively. Upon transformation, any PWM used in the algorithm must be reduced to the given R 2 and K d . The matrix transformation procedure is described in detail in [33]. However, it should be noted that R 2 cannot be set constant, since in this case the number of cells in the PWM would increase and the cell values would not remain in approximately the same order (as needed) but would rather tend to approach zero in order to correspond to Equation (1). Therefore, R 2 was specified as a function of the period and cardinality of the SI sequence alphabet: R 2 = R L A in L, where R L is the multiplier parameter, which can be used to scale R 2 .
As the sequences in set SI may have different lengths, we created sets of matrices Q(L) with length L ranging from 0.9L to 1.1L (see Section 4.1). Each Q(L) set contained 500 PWMs.

Using a Genetic Algorithm to Optimize PWMs
Then, we calculated PWM m using a genetic algorithm described in detail in [18,33] ( Figure 1, step 3). Briefly, PWMs of sets Q(L) with the size |Q| = 500 PWMs (Section 4.2) were considered as organisms. To calculate the objective function for each PWM from set Q(L), global alignment with each sequence from set SI was carried out. We calculated similarity function Fmax(l) at point (L1(l), L), where l is the sequence number in set SI, L is the number of columns in the PWM, and L1 is the length of the l-th sequence from set SI.
Then, the mF value was calculated and used as an objective function for PWM matrices from set Q(L): where k is the index of the matrix from set Q(L). The PWMs were ranked according to mF in the descending order, and the matrix with the largest mF value, mF(1), was saved, whereas two PWMs with the smallest mF values were excluded from set Q(L). Then, we created two children from the remaining PWMs; as before [18,33], the offspring were generated by gluing of two parent matrices, which were chosen randomly. However, the probability of selecting a particular matrix increased with its mF value. In addition, random mutations were introduced into 50 randomly selected PWMs (except children created at this step) by replacing one random element by a random value evenly distributed in the interval from −10 to 10. After these changes, a new Q(L) set was generated, its matrices compared with the sequences from set SI, and a new vector mF(k) (k = 1, 2, ..., 500) was obtained. The procedure of modifying set Q(L) was iterated until mF(1) ceased to increase during the last 10 cycles. As a result of the genetic algorithm for set Q(L), mF(1) (denoted as mF L (1)) and the corresponding PWM were obtained.
The iterative procedure was performed in each Q(L) set for L from 0.9L to 1.1L. Finally, we chose the L value for which mF L (1) was maximal; it was denoted as L m and the corresponding matrix as PWM m .

Global Alignment of PWMs from Set Q and Sequences from Set SI
Next, sequence SI(l) of length L1(l) was aligned with a PWM of number k from set Q(L). Sequence SI(l) was denoted as S3 and its elements as s3(i), where i = 1, 2, ..., L1(l). We also used sequence S4 with elements s4(j), where j = 1, 2, ..., L.
To construct the alignment of sequences s3(i) and s4(j), we used the global alignment algorithm [51] with an affine gap penalty function and PWM values from set Q(L) [33] instead of substitution matrices. Then, similarity function F was calculated as: where d and e are gap opening and gap extension penalties, respectively. Fmax(l) = F(L1(l),L) (Section 4.2) was used as a measure of similarity between sequences S3 and S4. We also filled in matrix F', where in each cell (i, j) we remembered the number of the cell from which we arrived at this cell using Formula (4). As a measure of similarity between set SI(l) and the PWM from set Q(L), we have chosen: The PWM m of length L m (Section 4.3) was aligned with set SI(l) using matrix F', and the results were applied to calculate multiple alignments of sequences SI(l) (l = 1, 2, . . . , N), which were evaluated according to mF(PWM m ) calculated with Formula (6).

Algorithm for Constructing Multiple Alignment
To obtain MA m for the SI sequences ( Figure 1, step 4), we used PWM m and a set of global alignments of the SI(l) sequences (l = 1, 2, ..., N,) with PWM m . For this, sequence S 4 with elements s 4 (j) (j = 1, 2, ..., L m ) was arranged horizontally, and under it the alignments of SI(l) sequences (l = 1, 2, ..., N) generated in Section 4.4 were written. If any of the SI(l) sequences had an insertion of k amino acids compared to sequence S 4 , then additional k columns were created for these amino acids in the s 4 (j) sequence (if they were not created for any other SI(l) previously). As a result, MA m was constructed.  (Figure 1, step 5), we used the Monte Carlo method. First, we generated 300 sets of random sequences Q r through random shuffling of residues in the SI sequences. For each Q r we calculated mF(PWM m ) using Formula (6), its arithmetic mean mF(PW M m ) and dispersion D(mF(PWM m )), and, finally, Z: A Z value greater than threshold Z0 indicated that the alignment obtained with PWM m was statistically significant.

Estimating the Statistical Significance of an Arbitrary MA
The algorithm described above can be applied to evaluate the statistical significance of any alignment. Let us take an MA with length K containing α sequences, each denoted as S k 5 and its elements as s k 5 (j) (where k = 1, 2, ..., α and j = 1, 2, ..., K). Let us also introduce sequence S6 with elements s6(j), j = 1, 2, ..., K representing consecutive numbers from 1 to K; S6 would be used as column numbers for MA. First, we transform MA to remove those columns where the number of amino acids is less than α/2 and that of deletions is more than α/2; as a result, multiple alignment MA' of length K' is obtained. At the same time, sequence S6 is transformed by eliminating the numbers equal to those of the columns to be removed. As a result, we obtain sequence S 6 with length K'. Then, we calculate the amino acid frequency matrix for MA' denoted as V(i,j), where i = 1, 2, . . . , 20 and j = 1, 2, . . . , K'; each element shows the number of type i residues at position j of MA'. Next, we calculate the PWM based on matrix V(i,j) as: PW M(i, j) = (V(i, j) − p(i, j))/ U p(i, j)(1 − p(i, j), where U is the total number of amino acids in MA', p(i, j) = X(i)Y(j)/U 2 , X(i) is the number of type i residues, and Y(j) is the total number of residues in the j-th column of MA'. Then, we transform PWM(i,j) according to Formulas (2) and (3).
After calculating PWM (i,j), the weight of MA can be determined. First, we compute the sum for all k = 1, 2, ..., α and j = 1, 2, ..., K: Sum = α ∑ k=1 K ∑ j=1 PW M (s k 4 (j), s 5 (j)). Next, we obtain the Del(i) vector, which shows the number of deletions with size i in MA. Finally, , which is considered as the MA weight. Then, statistical significance should be evaluated by estimating mean F MA and variance D(F MA ). For this, symbols characterizing deletions are removed from S k 5 sequences, which are then randomly shuffled to yield S k Rand sequences, where k = 1, 2, ..., α, and the statistical significance is determined for matrix PWM (i,j) and sequences S k Rand according to Formula (7), assuming that SI = S k Rand , N = K, and PWM m = PWM (i,j). Thus, we calculate the Z value for MA, which enables comparison of different MAs by their statistical significance.

Creation of Artificial Sequences to Compare Different Methods of Constructing MAs
To compare the performance of MA methods T-Coffee, MUSCLE, Clustal Omega, Kalign, MAFFT, and PRANK with that of the MAHDS algorithm, we used artificial sequences, which can be created with a given number of amino acid substitutions and indels for convenience of analysis. To generate an artificial sequence, we first created a random ancestor sequence Anc of length L and then a set of descendant sequences Des(i) (i = 1, 2, . . . , 100) by adding a given number of random substitutions and/or indels to Anc in a random order. In case of random substitutions without indels, the lengths of child sequences were the same and equal to L, and the number of substitutions in each Des(i) sequence with respect to Anc was denoted as s 1 . Let us show that the presence of s 1 random substitutions in each Des(i) sequence leads to 2s 1 random substitutions between any two Des sequences.
When one random substitution is made in sequence Anc, the probability that in sequence Des(1) with the number of substitutions s 1 a residue at position j will be changed is 1/L. Then, the probability that after s 1 replacements residue j will not be replaced is: We assume that amino acid substitutions during the creation of the Des(1) sequence occur with equal probability. For the rest of amino acids, it is likely that during the last replacement in a given position, the original amino acid will appear as substitution. The probability of falling into this category is: where |A in | is 20. Thus, the probability that residue j in sequence Des(1) will match the residue at Anc is: The events described in Formulas (9) and (10) are applicable to all situations, when a residue at position j remains unchanged after random substitutions. Now, let us consider the case when sequences Des(2) and Des(3) are independently generated from sequence Anc by making s 2 random substitutions. Then, for both sequences, Formulas (8)-(10) can be presented as: P 02 = 1 − 1 L s 2 , P 12 = 1 |A in | (1 − P 02 ), P m2 = P 02 + P 12 . Let us take P m3 as a probability that in Des(2) and Des(3) two residues could match. If P m3 = P m , then the evolutionary distance between Anc and Des(1) is equal to that between Des(2) and Des (3), which means that, on average, the same number of amino acids coincide between Anc and Des(1) and between Des(2) and Des(3). In Des(2) and Des(3), the proportion of residues unchanged compared to Anc is P 02 ; then, the proportion of residues preserved at the same positions of child sequences is P 03 = P 2 02 , since it is the probability for a residue to remain in both Des(2) and Des (3). Other amino acids matching in Des(2) and Des(3) include those that have been substituted with the same residues in both sequences, and the probability of such events is P m4 = 20 ∑ k=1 (1 − P 02 ) 2 /|A in | 2 = (1 − P 02 ) 2 /20.
Finally, there is a probability that a residue that remained from Anc in Des(2) could match the one randomly replaced in Des(3) and vice versa. The sum of the two probabilities is: P m5 = 2 20 ∑ k=1 P 02 (1 − P 02 )/|A in | 2 = 2P 02 (1 − P 02 )/20. Then, P m4 + P m5 = 1 20 1 − P 2 02 and P m3 = P 2 02 + (1 − P 2 02 )/20. If P m3 = P m , then P 2 02 + (1 − P 2 02 )/20 = P 0 + (1 − P 0 )/20, i.e., P 2 02 = P 0 . If the values for P 02 and P 0 are substituted, then: Formula (11) means that the number of s 2 random substitutions introduced in Des(i) sequences is equivalent to twice the number of s 2 mutations in pairwise comparison of Des(i) sequences relative to Anc, which happens under the condition that all child sequences in set Des(i) (i = 1, 2, ..., 100) are generated independently of each other. We used this property in the calculation of the average number of mutations between Des(i) sequences.

Conflicts of Interest:
The authors declare no conflict of interest.