Detection of tandem repeats in the Capsicum annuum genome

Abstract

In this study, we modified the multiple alignment method based on the generation of random position weight matrices (RPWMs) and used it to search for tandem repeats (TRs) in the Capsicum annuum genome. The application of the modified (m)RPWM method, which considers the correlation of adjusting nucleotides, resulted in the identification of 908,072 TR regions with repeat lengths from 2 to 200 bp in the C. annuum genome, where they occupied ~29%. The most common TRs were 2 and 3 bp long followed by those of 21, 4, and 15 bp. We performed clustering analysis of TRs with repeat lengths of 2 and 21 bp and created position-weight matrices (PWMs) for each group; these templates could be used to search for TRs of a given length in any nucleotide sequence. All detected TRs can be accessed through publicly available database (http://victoria.biengi.ac.ru/capsicum_tr/). Comparison of mRPWM with other TR search methods such as Tandem Repeat Finder, T-REKS, and XSTREAM indicated that mRPWM could detect significantly more TRs at similar false discovery rates, indicating its superior performance. The developed mRPWM method can be successfully applied to the identification of highly divergent TRs, which is important for functional analysis of genomes and evolutionary studies.

1. Introduction

In eukaryotic organisms, a considerable part of the genome is occupied by dispersed and tandem repeat (TR) sequences. In the human and plant genomes, repeats, and transposable elements constitute about 43% and 80%, respectively.¹ TRs can be found in the promoters, 3ʹ untranslated regions, exons, and introns, and it is suggested that repeats are involved in DNA compaction and evolutionary events through contribution to gene expression divergence.^2,3 Most TRs can mutate quickly, thus adding to genetic variability and promoting rapid adaptation of organisms to new environmental conditions.^4,5 TRs present in centromeric, pericentromeric, and telomeric regions^6–9 and areas of chromosome breakage can lead to chromosomal rearrangements.

TRs are generally classified into microsatellites (or short TRs [STRs]) and minisatellites with repeat lengths of 1–6 and 7–60 bp, respectively. TRs longer than 60 bp (and even longer than 1,000 bp) are also observed.¹⁰ Repeat sequences may be species-specific, and their numbers could vary considerably owing to DNA polymerase slippage and unequal recombination.¹¹ Short TRs are detected at the ends of transposable elements: SINEs, LINEs, and retrotransposons, where they affect the possibility of transposition, thus exerting genetic regulation. It is noted that more than 25 human hereditary diseases are associated with STRs,¹² including fragile X syndrome^13,14 and autism.^15,16

There are cases of correlation between the TR length and function. Thus, 3 bp period is characteristic for DNA coding regions and is associated with the triplet nature of the codons,¹⁷ whereas TRs of 10–11 bp (referred to as 10.5 bp) has been demonstrated in nucleosomal DNA.^18,19 TRs of 5–7 and 11–14 bp are characteristic for enhancers and other non-coding regions adjacent to genes, whereas long (120.9 bp) TRs are observed at the 5ʹ ends of the transcripts.²⁰

The diverse regulatory roles of TRs have spurred the development of computer programs for TR search in the genomes.²¹ Some of them can identify TRs with insertions or deletions (indels), provided that the degree of similarity between repeats is high (>50%). For example, Tandem Repeat Finder (TRF)²² searches the candidates based on the statistics of k-tuple matches, and the TR is recreated considering possible indels, when k-tuples are placed relative to each other. As the presence of indels within k-tuples complicates TR detection, in the TRStalker algorithm,²³ gapped q-grams allowing indels are tracked, resulting in the detection of long (>10 bp), more divergent repeats. T-REKS²⁴ and XSTREAM²⁵ using dynamic programming to build multiple alignments can also find TRs with indels. However, these methods cannot detect TRs with similarity < 40%.

Highly divergent (fuzzy) TRs are of particular interest because of their presence in gene regulatory regions, where they interact with transcription factors.²⁶ Fuzzy TRs can be identified with the methods that use information decomposition^21,27 and the Fourier transform (FT).^20,28,29 Such approaches have revealed that short-length (2, 3, 4, and 6 bp) repeats are characteristic for introns and exons.³⁰ However, the disadvantages of the FT and information decomposition methods are that indels are not considered.

Thus, there is currently a need for a universal method that can identify long TRs with a large number of accumulated indels. In our previous studies, we have developed the random position weight matrix (RPWM) method, which could effectively detect TRs with the average number of mutations per nucleotide (x) up to 3.2.^31,32,33 Here, we improved the RPWM method by considering the correlation between adjacent DNA bases, which allowed identification of TRs with 2–200 bp, and refined the estimation of the statistical significance of the found TRs by using the Monte Carlo method. As a result, it was possible to increase the average number of TR-containing regions detected per 10⁶ DNA bases by 1.5 times.

The modified (m) RPWM method was used to search for TRs in the pepper (Capsicum annuum) genome. Although the genome of C. annuum was sequenced in 2014,³⁴ no detailed analysis of its TRs has been performed except for the study showing that 15.64% of expressed sequence tags contain STRs (2–6 bp)³⁵ which have been subsequently used to develop various polymorphic microsatellite markers.³⁶ By applying the mRPWM method, we identified all possible TRs with a length 2–200 bp in the C. annuum genome and created a corresponding database (http://victoria.biengi.ac.ru/capsicum_tr/).

2. Materials and methods

2.1. RPWM algorithm

The pepper (Capsicum annuum) genome was used as a model to search for TRs with the improved RPWM algorithm. RPWM, which uses position weight matrices (PWMs),³³ was modified to consider the correlation of adjacent nucleotides in TR detection and more accurately determine the statistical significance of the identified TRs. In order to understand the nature of the improvement, let us review the main steps of the RPWM algorithm.

Step1. A window with length L = 650 b denoted as sequence S moved along a chromosome with a step of 10 b, and TRs of length n = 2–50 b were searched for in each S.

Step 2. A set of 500 random PWMs with dimensions of n × 4 (TR length × 4 DNA bases) was generated and denoted as Q_n; N_o = 500 is the volume of set Q_n. Each matrix Q_n(i) was transformed so that sum $∑i=1n∑j=14pwm(i,j)2$ was always constant. In this formula, pwm(i,j) is the element of the PWM on row i and column j, p₁(i) = 1/_n, and p₂(j) = N(j)/L, where N(j) is equal to the number of A, T, C, or G bases in sequence S. Sum $∑i=1n∑j=14pwm(i,j)p1(i)p2(j)$ was also kept constant for all matrices from set Q_n. This transformation was necessary to ensure that similarity function Fmax (step 3 below) had similar distribution for different n³⁷; otherwise, it would be very difficult to find the most statistically significant period and determine TR length in sequence S.

Step 3. For each set Q_n, we applied the genetic algorithm to find matrix $Qnmax$ with the best local alignment to sequence S; the maximum value of similarity function F_max was used as the object function and matrices from set Q_n – as organisms. The genetic algorithm was applied as follows. A. Each matrix from set Q_n was locally aligned to sequence S as previously described,^33,38 and F_max and local alignment coordinates in sequence S were obtained for each matrix, which allowed accurate determination of TR coordinates in sequence S; F_max values were written in vector V_n(i) (i = 1, 2,…, N_q). B. Vector V_n(i) was ranked in a descending order so that V_n(1) contained the maximum value of F_max, and random mutations were introduced into 10 randomly selected matrices from set Q_n by replacing one randomly selected cell with a random number from −10 to +10. Two matrices were randomly selected as ‘parents’ and randomly superposed to create a new ‘descendant’ matrix; then, a matrix with the smallest F_max value (contained in V_n(500)) was removed from set Q_n. Then, we went back to point A and recalculated all V_n values for the mutated matrices and the ‘descendant’; for the matrices that did not change in point B, local alignment was not performed and V_n was not recalculated. The A–B cycle was repeated until V_n(1) tended to increase. As a result, we obtained matrix V_n(1) and local alignment of sequence S for a TR of length n.

Step 4. Steps 2–3 were performed for n = 2–50 bases to calculate all V_n(1) values and determine n with the maximum V_n(1) value, which was designated as n_max(k), where coordinate k is the beginning of the window (sequence S) in the chromosome (step 1). Thus, for each window, we obtained n_max(k), $Vnmax(k)(1)$⁠, matrix $Qnmax$⁠, and local alignment of matrix $Qnmax$ and sequence S.

Step 5. Windows in sequence S were shifted by 10 bases to overlap with each other, and $Vnmax(k)(1)$ values were filtered as follows: the dependence of V(1) on the TR length (n) was excluded, resulting in $Vec(k)=Vnmax(k)(1)$⁠, and the local maxima in vector Vec(k) were chosen so that Vec(k) > V₀ and Vec(k-i) ≤ Vec(k) ≥Vec(k+i), where i ranges from 1 to 60. The selected local maxima represented the final result. Thus, for each local maximum, we obtained n_max(k), Vec(k), matrix $Qnmax$⁠, and local alignment of matrix $Qnmax$ and sequence S.

Step 6. A random sequence S was used to select threshold level V₀ for local maxima Vec(k). Steps 1–5 were repeated for this sequence, and a set of local maxima Vec_r(k) was obtained; the numbers of local maxima Vec(k) and Vec_r(k) were denoted as N_v and N_r, respectively. Then, we chose V₀ with the ratio N_r(Vec(k) > V₀)/ N_v(Vec(k) > V₀) ~ 0.027. These local maxima and the corresponding n_max(k), matrix $Qnmax$ together with local alignment of the matrix and sequence S represented the final result of the algorithm.

2.2. Modified RPWM algorithm

2.2.1 Calculation of statistical significance Z(n)

We modified RPWM³³ to consider the statistical significance of the identified TRs. Matrices Q_n(i) transformed in section 2.1 allowed comparison of V_n(1) for different n, which could be performed if F_max(n) defined by local alignment of sequence S and matrix Q_n(i) had the same distribution function. However, it is not possible to achieve complete identity of these distributions in the RPWM method. Therefore, to reduce errors in determining n_max(k), we calculated Z(n) for each period length n using the Monte Carlo method; for this, steps 4–6 in the RPWM algorithm (section 2.1) were replaced with the new ones described below.

Step 4. For each TR of length n, we applied the Monte Carlo method to estimate statistical significance Z(n) of V(1): $Z(n)=(V(1)−V(1)¯)/(D(V(1))0.5$⁠, where $V(1)¯$ and D(V(1)) are the mean value and variance, respectively, of V(1) obtained by aligning random sequences S and matrix $Qnmax$⁠.

Step 5. Steps 2–4 were performed for n from 2 to 50 to calculate all Z(n) values and determine n with the maximum Z(n) value, which was designated as n_max(k), where coordinate k is the beginning of the window (sequence S) in the chromosome (step 1). Thus, for each window, we obtained n_max(k), Z(n_max(k)), matrix $Qnmax$⁠, and local alignment of matrix $Qnmax$ and sequence S.

Step 6. Windows in sequence S were shifted by 10 bases to overlap with each other, and Z(n_max(k)) values were filtered as follows: the dependence of Z on TR length (n) was excluded, resulting in Z(n_max(k)) = Z₁(k), and the local maxima in vector Z₁(k) were chosen so that Z₁(k) > Z₀ and Z₁(k-i) ≤ Z₁(k)≥ Z₁(k+i), where i ranges from 1 to 60. The selected local maxima represented the final result. Thus, for each local maximum, we obtained n_max(k), Z(n_max(k)), matrix $Qnmax$⁠, and local alignment of matrix $Qnmax$ and sequence S.

2.2.2. Consideration of the correlation of neighbouring bases

Matrix $Qnmax$ contains n rows and 4 columns (n × 4) and does not take into account the correlation of neighbouring nucleotides in sequence S, because the columns of the matrix are not interconnected; as a result, a significant part of TRs may remain undetected. To address this problem, we created an improved version of the RPWM algorithm, which should consider the correlation of adjacent nucleotides. For illustration, let us search for TRs with the length of 4 b (n = 4) in sequence S of length L = 1,200 b. To obtain S, we first randomly selected, with a probability of 0.25, one sequence Seq(i) (i = 1–4; Seq(1) = ATCG, Seq(2) = TAGC, Seq(3) = CCAA, and Seq(4) = GGTT) and filled in bases 1–4 in sequence S; the step was repeated until all the bases in sequence S (1–1,200) were filled. Then, we determined function Z(n) using the RPWM method.³³ The results indicated that the 4 base TRs in S could be hardly seen (Fig. 1), because RPWM computed matrix M(n,4) for Z, which had n rows and 4 columns, i.e. for sequence S, at each position of repeat n, the frequency of each of the four bases (A, T, C, or G) should be the same and equal to 0.25.

An example of such sequence S shown in Fig. 2A indicated that the first column contains a number of A, T, C, and G equal to 3. The same frequencies were observed for the second, third, and fourth positions of the periods; the resulting matrix is shown in Fig. 2B. These data mean that there are no differences in base frequencies at different positions of the repeat, which makes it difficult to detect TRs in sequence S according to base frequencies in matrix M. Therefore, low Z values were obtained using the RPWM method (Fig. 1); they were not zero only because dynamic programming arranged indels in a certain way.

In contrast, after the application of the modified (m)RPWM algorithm, 4 b-long TRs were clearly visible and Z(4) had the largest value among all repeat lengths (Fig. 1), because in mRPWM, instead of matrix M(n,4) (Fig. 2B, top), we used matrix M(n,16), which had 16 columns corresponding to 16 nucleotides (Fig. 2B, bottom). The elements of matrix M(n,16) were calculated as m(v(j),i) = m(v(j),i) + 1 (where j ranged from 2 to L, i = s(j − 1) + 4(s(j) − 1), and v(j) was an element of sequence V), indicating that matrix M(n,16) considered correlations of adjacent nucleotides. Matrix M(4,16) for sequence S (Fig. 2B, bottom) was extremely heterogeneous, and many of its elements were equal to zero. Thus, the newly developed mRPWM could identify TRs of 4 b in sequence S, i.e. appeared to be more efficient in search of TRs than the original RPWM.

In mRPWM, we considered the correlation of neighbouring bases by creating a set of matrices Q_n with dimensions n × 16, which had n rows and 16 columns instead of n rows and 4 columns, as described in section 2.1. Step 1; in addition, the size of the window was increased from 650 to 1,200 b and that of the step – from 300 to 600 b. The length (n) of TRs to be detected was also increased: from 2–50 to 2–200 b. In this analysis, we searched for TRs of all sizes, including those that were multiples of 3 b. As a result, the number of detected TRs was significantly increased. Thus, RPWM could identify approximately 76,000 TRs in the genome of rice (Oryza sativa), whereas mRPWM – 908,072 TRs in that of C.annuum. Given that the sizes of the two genomes are approximately 3.75 × 10⁸ and 3 × 10⁹ b, the average numbers of identified TRs per 10⁶ b were 192 and 300, respectively, confirming the superior performance of mRPWM.

2.3. Effect of nucleotide substitutions on TR identification

Artificial sequences of 1,200 b, which contained TRs of 10 b, were generated through tandem multiplication of random 10 b segments by 120 times; then, 0, 100, 500, 1,000, 1,500, 1,800, 2,200, and 2,400 random base substitutions and 20 indels were introduced at randomly selected positions. These artificial sequences were denoted as St(x) (where x was the average number of substitutions per base between any two repeats and was equal to 0, 0.17, 0.85, 1.7, 2.5, 3.0, 3.7 and 4.0)³⁹ and used to search for TRs with RPWM and mRPWM and determine Z(10) as a function of x.

3. Results

3.1. Analysis of artificial sequences St(x)

Figure 3 shows the plot of Z(10) versus x. If we take Z₀ (section 2.1, step 6) equal to 6.0 (section 3.2), then RPWM can detect TRs with x up to 3.2, whereas mRPWM, which considers the correlation of neighbouring bases, can detect TRs with x up to 3.6, which significantly exceeds the limit for TRF and T-REKS methods: x < 1.2.³³ Thus, the performance of mRPWM in identifying highly divergent TPs is superior to those of other currently used methods.

3.2. Detection of TRs in the C. annuum genome

Next, we applied the mRPWM method to search for TRs with lengths of 2–200 bp in the 12 chromosomes of C. annuum (sequences retrieved from Ensembl Plants databank, http://plants.ensembl.org/Capsicum_annuum/Info/Index). To select the significance level Z₀, we searched for TRs in the first C. annuum chromosome and in a random sequence created by random mixing of chromosome nucleotides. The Z values for the sorted local maxima (section 2.1, steps 5 and 6) are presented in Table 1; among them, Z₀ = 6.0 was chosen as the threshold level which provided a false discovery rate (FDR) of 2.71%. The number of TRs found in each C. annuum chromosome is shown in Table 2. Overall, 908,072 TRs were detected; their density was similar in all the chromosomes and constituted ~302 TRs per 10⁶ bp. The calculations were performed on two Ryzen 9 5950X processors, and the search for TRs in the C. annuum genome took about 2 weeks.

The statistics for the 10 most common repeat lengths indicated that TRs of 3 and 2 bp were the most frequent in the C. annuum genome (Table 3), which is consistent with the data reported for other plant genomes.^31,40 Analysis of TR distribution depending on the length revealed that in addition to TRs with n = 2 and 3 bp, there were local peaks for those with n = 21, 31–33, 91–92, 111–112, and 178–187 bp (Fig. 4); the existence of lengthy TRs may be associated with chromatin compaction. It has been reported that in proteins, the most frequently observed repeat lengths are 2 and 7 amino acids,⁴¹ which correspond to 6 and 21 bp, respectively.

The performance of the newly developed mRPWM method was verified by searching for highly dissimilar TRs that could not be detected by the existing algorithms. TRs were analysed in a non-coding region of the first C. annuum chromosome (15,478,656–15,479,825 bp). The Z(n) plot indicated that a pronounced maximum (Z = 9.0) was observed for n = 7 bp and lower peaks were also detected for n = 14, 21, 28, and 35 bp (Fig. 5). Table 4 shows the alignment between a TR and an artificial periodic sequence containing 34 indels with a maximum length of 3 bp; the corresponding $Q7max$ matrix is shown in Table 5.

Supplementary material files Period_5.tif and Period_21.tif show the results for two other regions of the first C. annuum chromosome: 21,950,674–21,951,081 bp and 5,039,380–5,040,582 bp, in which periodicity of 5 and 21 bp, respectively, was detected; the corresponding matrices $Q5max$ and $Q21max$ can be found by coordinates in the database at http://victoria.biengi.ac.ru/capsicum_tr/. We also constructed multiple alignments for the three analysed regions of the first C. annuum chromosome (Supplementary material files Multi7.txt, Multi5.txt, and Multi21.txt).

Cumulatively, these results indicate that the mRPWM method can identify such TRs that cannot be detected by the other existing algorithms.

3.3. Presence of TRs in the annotated C. annuum sequences

Next, we analysed the distribution of the detected TRs in the C. annuum genes (http://ftp.ensemblgenomes.org/pub/plants/release-53/gff3/capsicum_annuum/), potential promoter sequences (PPSs; http://victoria.biengi.ac.ru/cgi-bin/dbPPS/index.cgi), and SINE repeats (http://victoria.biengi.ac.ru/sine_pepper/index/) by calculating the number of TRs that overlapped with any annotated sequence by at least 50%. The results indicated that the number of TRs intersecting with the genes was 42,210, which exceeded that of the C. annuum genes (31,600), i.e. one gene contained several TRs. The number of genes overlapping with TRs was 20,739, indicating that approximately 2/3 of the annotated C. annuum genes contained repeated sequences. Among these genes, 80% (16,656) had TRs with repeat lengths that were multiples of 3 bp (including 14,275 with the length of 3 bp), which is consistent with a previous report that over half of all genes have 3 bp TRs.^42,43

Among the 825,136 PPSs⁴⁴ and 50,077 SINEs (unpublished data) found in the C. annuum genome,⁴⁴ 277,929 and 7,444, respectively, contained TRs, among which those of 2, 3, and 21 bp long were the most prevalent, corresponding to the overall statistics of TR length distribution in the genome (Table 3).

3.4. Cluster analysis of TRs with lengths 2 and 21 bp

Among the most common TRs in the C. annuum genome (Table 3), those with n = 3 bp have been studied in detail^42,43; therefore, we focused on TRs with n = 2 and 21 bp and analysed their heterogeneity, i.e. similarity of TRs from different regions. For this, we performed cluster analysis of n × 16 matrices $Qnmax$ for n = 2 and 21 bp (section 2.1, step 6), in which matrix $Qnmax$ was considered as a point in the Euclidean space n × 16, and determined the difference between two matrices (designated as $m1(i,j)$ and $m2(i,j)$⁠) according to the Euclidean distance. As a TR-containing region can begin at any repeat position, all variants should be taken into account for cyclic shifts between matrices. In the example presented in Fig. 2A, sequence V started from the first position, which was arbitrary, because it could start from any position in the repeat (for example, 3412341234…). Therefore, to calculate Distⁿ, which is the smallest distance between two matrices, we cyclically shifted the beginning of one of them by 0, 1, 2, …, n − 1 positions so that all possible options for matching matrices were considered. In total, there were n distances between two matrices:

where k(t) = i + t − n{int((i + t − 0.5)/n)}.

The minimum was found among all values of t from 0 to n − 1. By determining Distⁿ between all pairs of matrices, we built a matrix of paired distances for all matrices $Qnmax$ obtained for a repeat of length n and then applied the hierarchical clustering algorithm. However, 156,442 and 18,194 TRs of 2 and 21 bp, respectively, resulted in very large distance matrices; therefore, to reduce matrix dimensions, we grouped 2,000 randomly selected matrices $Qnmax$ into sets W₂ and W₂₁ and performed clustering using the statistical package R and the Complete Linkage algorithm.

The results are shown in Figs 6 and 7. Distance values $Dist02$ = 17 (n = 2) and $Dist021$ = 87 (n = 21) were taken as group-forming levels. To select $Dist0n$⁠, we randomly shuffled the cells in each matrix of sets W₂ and W₂₁ to obtain sets of random matrices $WR2i$ and $WR21j$ (i, j = 1…100) and subjected them to cluster analysis, which produced 100 random matrix clusterings per set. Then, each clustering was analysed for the number of classes and matrices in each class at levels $Dist02$ = 17 and $Dist021$ = 87. For random set $WR2i$ (i = 1…100), the average number of matrices in the class (⁠$N¯2$⁠) was 43.8 and variance D(N₂) was 3.88. The statistical significance of creating classes was determined as $Zn=(Nn−N¯n)/(D(Nn))0.5$⁠. At $Dist02$ = 17, there were six classes (N₂ = 6) containing 339, 135, 293, 224, 663, and 346 matrices, respectively, and Z₂ for each class was > 20.0. For random set $WR21j$ (i = 1…100), $N¯21$ = 6.1 and D(N₂₁) = 0.009. At $Dist021$ = 87, N₂₁ = 6, the number of matrices in each class was 57, 258, 185, 124, 1028, and 348, respectively, and Z₂₁ for each class was also > 20.0. Class matrices for TRs with repeat lengths of 2 and 21 bp are given in Supplementary materials (groups2_17.txt and groups21_87.txt, respectively).

3.5. Database of C. annuum TRs

The TRs detected in the C. annuum genome are placed into a database (http://victoria.biengi.ac.ru/capsicum_tr/), which provides the following information for each TR-containing region: chromosome number, repeat length n_max, Z(n_max), coordinates in the chromosome, sequence, matrix $Qnmax$ and local alignment of the matrix. The available sequence annotation is also shown if it overlaps by at least 50% with a TR-containing region. Site http://victoria.biengi.ac.ru/splinter/login.php provides the application to search for TRs in DNA sequences using the improved algorithm mRPWM.

4. Discussion

In this study, we modified the RPWM method developed earlier³³ so that the new version, mRPWM, could detect TR regions with much longer repeats (2–200 bp) than the original RPWM (2–50 bp). mRPWM is a novel method which detects TRs using only the DNA sequence and does not require a library of previously defined TRs as is the case with HipSTR.⁴⁵

To compare the performance of mRPWM and the original RPWM, we applied them to search for TRs in the C. annuum genome. In the first C. annuum chromosome, RPWM could detect 69,546 regions with different period lengths, whereas mRPWM detected 97,727 of such regions (Table 2), which is about 40% more. For the period length from 2 to 50 bp, mRPWM detected 82,427 TRs, indicating that 15,300 regions identified by mRPWM contained repeats with the length between 51 and 200 bp.

To match TRs found by RPWM and mRPWM, we analysed their intersection, i.e. the ratio of the total length of the two TR-containing regions to the length of the smallest one (Table 6). The results revealed that 46,093 regions (67% of those detected by both methods) overlapped considerably (intersection C > 0.5), whereas 7,803 (~11%) overlapped weakly, and 15,560 (~22%) were unique, indicating consistency between the results obtained by the two RPWM methods. However, the coincidence of period lengths in the overlapping areas constituted only ~11% (8,047 regions), which could be attributed to the following reasons.

First, to estimate the statistical significance of the identified TRs, RPWM uses Vec(k), which contains F_max and which is calculated by dynamic programming (section 2.1, step 3) and must be greater than threshold value V₀.³³ The local maximum is determined using dynamic programming and the PWM, and the distribution function for F_max strongly depends on TR length n and PWM mQ,³³ which means that at the same F_max value, probability P(F > F_max) may differ for different n and mQ despite special transformations of matrix mQ.^33,37 To speed up the calculations, among n = 2–50 bp we have chosen n with the largest F_max.³³ However, it is not quite correct to do so when the distribution function for F_max depends on n; in this case, n with the lowest probability P(F > F_max) should be chosen. To estimate P(F > F_max) in this study, we used the Monte Carlo method for calculating statistical significance Z(n) (section 2.2.1, step 4), which allows a more accurate determination of the period length for TRs and which accounts for the differences in period lengths observed in the overlapping regions.

Second, in the intersecting regions, there are periods t for which t/n₀ is an integer equal to or greater than 2 (here, n₀ is the length of TRs). Vec(k) for n = t (section 2.1, step 5) may have the maximum value among all periods, whereas Z(n) (section 2.2.1, step 4) may be maximal for n = n₀ because of weak nucleotide correlation at n ≥ n₀. Therefore, the lengths of TRs identified by the RPWM and mRPWM methods would be different. As an example, we used an artificial sequence {ATCGATTCGG}₁₂₀ (1,200 nt long), in which 1,200 random substitutions were introduced at random positions (Supplementary material, file seq_period_10.txt); for this sequence, n_max(k) = 20 bases and Vec(k) = 2,480, whereas n₀ = 10 bases and Z(n₀) = 11.8 (full ranges of Z(n) and Vec(k) values depending on period length n are shown in Supplementary material files 10.new and 10.old). It appeared that RPWM and mRPWM showed different period lengths if the number of accumulated base substitutions (x) was large and that mRPWM could determine n₀ more accurately.

Third, it is possible that in some TRs a short period t is nested in a longer period n₀ without dividing the latter. In this case, RPWM can detect TRs of length n_max(k) = t, whereas mRPWM can find those with length n₀ because of base correlation. Let us illustrate it on an example of sequence {ATCGATCGATCGATCGCGG}₅₇ which simultaneously contains both short and long periods of lengths 4 and 21 bases, respectively, and in which we made 1,200 random substitutions in random places (Supplementary material file seq_period_21.txt.). As a result, n_max(k) = 4 and Vec(k) = 2,634, whereas n₀ = 21 and Z(n₀) = 9.2 (full ranges of Z(n) and Vec(k) values are shown in Supplementary material files 21.new and 21.old). Thus, mRPWM, by taking into account the correlation of neighbouring nucleotides, is able to detect much longer periods than the original RPWM.

Pepper was chosen for analysis as a popular vegetable crop widely cultured around the world, whose genetic structure and properties are of interest to agriculture, the food industry, and medicine. It is known that 81% of the C. annuum genome are repetitive sequences, of which 76.4% are transposons,⁴⁶ and more than 70% of the latter are long terminal repeat elements. The pepper genome has been previously analysed for the presence of dispersed repeats and transposable elements³⁴; however, there are no data on TRs. Using mRPWM, we identified 908,072 TRs with an average length of 917 bp, indicating that TRs occupy about 29% of the C. annuum genome. However, this number is a minimum estimate, because we targeted regions containing > 4 repeats with a length of 2–200 bp. Most TRs had a repeat length of 2 or 3 bp; the latter may be partly related to coding sequences,⁴⁷ because 2/3 of the genes contain TRs, of which 80% have a repeat length multiple of 3 bp. In contrast, 2 bp TRs are mostly found in non-coding regions and only 4.6% of them intersect with the genes. TRs with a repeat length of 21 bp are also relatively numerous in the C. annuum genome; it is suggested that weakly similar TRs with the length 10.5 bp (10.5 × 2 = 21 bp) may play a role in DNA packaging within chromatin¹⁹ or be involved in the formation of alpha helices in proteins.⁴⁸

Cluster analysis of TRs with n = 2 or 21 bp resulted in the identification of several groups and the creation of an averaged matrix for each. Such matrices could be used to find weakly similar TRs of a certain type using hidden Markov models or dynamic programming, which may be important for evolution studies.

A large number of algorithms have been developed for TR search in DNA and proteins.^{22–25,29,49} Some of them perform best in finding short perfect repeats, whereas the others are focused on weakly similar TRs with a relatively long period.⁴¹ Comparison with the other methods such as TRF,²² T-REKS,²⁴ and XSTREAM²⁵ could be useful for assessing the performance of mRPWM. However, different methods use various input parameters, which could affect the number of detected TRs⁵⁰; therefore, it is necessary to apply a common criterion of TR recognition for objective comparison. As in mRPWM we used the FDR (obtained by analyzing real and random DNA sequences) and searched for TRs in a 1,200 b window that moved along the chromosome with a step of 600 b, the same parameters were applied to compare the performance of the other algorithms with that of mRPWM in identifying TRs in the first C. annuum chromosome. Fragments from 600i + 1 to 1200 + 600i (i = 0.1…) were cut out, windows containing character N were removed, and the resultant windows were designated as set Win; then, the three methods mentioned above were used to detect TRs in sequences from set Win and to determine local maxima (section 2.1, step 6). To estimate the FDR, TRs were also searched in random sequences of set WinR obtained by random mixing of nucleotides in all windows from set Win. Various parameters of TRF and T-REKS were considered: for the former, we changed the alignment weight and for the latter, we took a different psim; the aim was to achieve FDR < 6.0%, which is not much higher than that for mRPWM (FDR = 2.71%, Table 1). XSTREAM was started with default settings. The detected TR regions were filtered to contain at least four repeats with a length of 2–200 bp in order to correspond with mRPWM search parameters. The results indicated that T-REKS, regardless of the psim parameter, could not find TRs with the FDR < 6.0%, whereas XSTREAM detected only 16,520 TR regions with more than 4 repeats at the FDR ~ 6.0%, and TRF found a total of 28,245 TRs (FDR = 2.82). mRPWM for number of repeats ≥ 4 identified 97,727 TRs at the FDR = 2.71% (Table 1). The reason for such a noticeable difference is that mRPWM does not look for similarity between individual repeats but uses an image of multiple alignments of all TRs.

We also compared the performance of mRPWM with those of GMATA and TANTAN. GMATA is a fast method to search for microsatellites in large genomes,⁵¹ and TANTAN is a program to identify low-complexity regions and STRs, which are then excluded from consideration in further sequencing studies to avoid their contribution to noise. First, we determined x values for which TRs could be detected in artificial sequences (section 2.3) and found that GMATA and TANTAN could detect TRs and low-complexity regions only for x ≤ 0.17 and x ≤ 0.85, respectively, but could not identify STRs with higher x. Then, we applied the two programs to search for TRs detected by mRPWM. For this purpose, we created a set of 10⁴ randomly selected sequences (Set_m) with length l(i) (i = 1, 2,…, 10⁴), which contained the TRs identified by mRPWM in the C. annuum genome. GMATA found only five Set_m sequences with a length greater than 50% of the corresponding l(i); in these cases, TR lengths were the same or multiples of those detected by mRPWM. In addition, GMATA found many (5,122) short sequences with an average length of 11.5 nt (< 2% of that of Set_m sequences) and n = 2–5 nt. These results indicate that GMATA is unable to identify the majority of TRs found by mRPWM and finds only their small fragments, which is likely due to the fact that the program can recognize only TRs with very few base substitutions (x ≤ 0.17). Often, GMATA is followed by HipSTR which genotypes STR alleles through comparison of the genome sequences of multiple polymorphic variants⁵²; however, HipSTR analyses only TRs found by GMATA or other de novo methods and cannot identify the repeats missed by these methods, indicating that most TRs detected by mRPWM also cannot be found by combining GMATA with HipSTR.

A similar situation was observed for TANTAN, which in Set_m identified only 8,191 short sequences with an average length of 34.6 nt and only 411 regions with a total length greater than 50% of the corresponding l(i). Thus, TANTAN, similar to GMATA, could not find most TRs detected by mRPWM and recognizes only small fragments.

We have already applied mRPWM to multiple alignments of promoter sequences⁵³ and TR search in the O. sativa genome.³³ Here, the method was improved to take into account the correlation of neighbouring nucleotides, which resulted in the increase of TR density per 10⁶ bp from 192 (O. sativa genome) to 302 (C. annuum genome), i.e. by > 1.5 times. Thus, the mRPWM method developed in this study shows superior efficiency in identifying highly divergent TRs with a wide range of repeat lengths in DNA sequences, which is important for functional and evolutionary analyses of the genomes.

References