Fast and global detection of periodic sequence repeats in large genomic resources

Abstract Periodically repeating DNA and protein elements are involved in various important biological events including genomic evolution, gene regulation, protein complex formation, and immunity. Notably, the currently used genome editing tools such as ZFNs, TALENs, and CRISPRs are also all associated with periodically repeating biomolecules of natural organisms. Despite the biological importance of periodically repeating sequences and the expectation that new genome editing modules could be discovered from such periodical repeats, no software that globally detects such structured elements in large genomic resources in a high-throughput and unsupervised manner has been developed. We developed new software, SPADE (Search for Patterned DNA Elements), that exhaustively explores periodic DNA and protein repeats from large-scale genomic datasets based on k-mer periodicity evaluation. With a simple constraint, sequence periodicity, SPADE captured reported genome-editing-associated sequences and other protein families involving repeating domains such as tetratricopeptide, ankyrin and WD40 repeats with better performance than the other software designed for limited sets of repetitive biomolecular sequences, suggesting the high potential of this software to contribute to the discovery of new biological events and new genome editing modules.


INTRODUCTION
Significant roles of repetitive DNA and protein sequences have been widely reported in both eukaryotes and prokary-otes. Transposable DNA elements are thought to be among the most important evolutionary driving factors that have been expanding within and between species' genomes via their copy-and-paste or cut-and-paste mechanisms (1,2). These repetitive elements induce large-scale genomic rearrangements and transcriptional regulation. Nonmobile short tandem repeat DNA sequences are also key elements inducing structural genome evolution in prokaryotic species. Tandem DNA repeats induce mispairing and slippage between repetitive units during DNA replication and drive genomic contraction and expansion (3). Intramolecular crossover DNA recombination is also promoted between tandem repeat regions of the genome (4). Some of these events are known to be reversible and lead to genomic phase variation, allowing cells and species to rapidly adapt to changing environments without having to undergo irreversible mutations (5).
Repetitive protein domains often serve as structural binding modules that stably interact with biopolymers. They are widely involved in protein folding and interactions in various biological processes. Tetratricopeptide repeats (TPRs) (6) and ankyrin (ANK) (7) repeats are large protein repeat families that are conserved from prokaryotes to eukaryotes. The repeat units of TPRs and ANK repeats are 34 and 33 amino acids (aa) long, respectively, and both are composed of a helix-turn-helix structure. These repetitive domains have been reported to mediate interactions with other proteins and RNAs and play important roles in cell cycle control, transcriptional regulation, translational inhibition, and protein translocation (6,7). WD40 repeat is another large protein repeat family found in both prokaryotic and eukaryotic species, but its functions are particularly well known in eukaryotes (8). A WD40 repeat is composed of seven-bladed ␤-propellers, where each propeller is around 40 aa long, involving four anti-parallel ␤-sheets, and serves as a scaffold for protein interaction. Accordingly, WD40 proteins coordinate multi-protein complex formation and underlie diverse biological functions such as signal transduction, transcriptional regulation, cell cycle control, chemotaxis, autophagy and apoptosis (8,9).
The structural code for proteins in general remains largely unclear and there have been major challenges in engineering these repeat protein modules to develop synthetic binding reagents for biomedical and nanotechnology applications (10)(11)(12). Most protein repeat sequences are 'imperfect' or 'degenerated,' where each repetitive unit contains variable amino acid residues and the degrees of repeat imperfectness vary widely. Some of these variable residues determine binding to specific biomolecules and deciphering this code for DNA binding has been extremely beneficial in the development of genome editing technologies (13,14). Several transcriptional regulators involve tandem protein repeats with specific periodicities to wrap around the double-stranded DNA helix. Cys2His2 zinc fingers (C2H2 ZNFs) are the most common DNA-binding motif found in eukaryotic transcription factors (15). C2H2 ZNFs represent periodic protein repeats that make tandem contact to targeting DNA sequence. The repeat unit size ranges from 28 to 30 aa and the variable amino acid residue pattern in each unit defines its binding to a specific DNA triplet (16). Similarly, transcription activator-like effectors (TALEs) of the type III secretion system encoded in the plant pathogenic bacteria of the Xanthomonas genus also have repeating domains (17). They are virulence proteins that bind to the host plant genomic DNA and hijack its gene expression system. The periodicity of the repeat unit ranges from 33 to 35 aa, where the combination of two variable amino acid residues at the 12th and 13th positions of the repeat sequence has a one-to-one relationship with a specific mononucleotide. By fusing DNA cleavage domains such as FokI endonuclease to C2H2 ZNFs and TALEs, the genome editing tools zinc finger nucleases (ZFNs) and TALE nucleases (TALENs), respectively, have been developed, both of which enable highly specific targeted DNA cleavage. Other effector proteins have also been fused to C2H2 ZNFs and TALEs to regulate gene expression and chromosomal structures in various organisms (18,19).
The CRISPR-Cas systems have become the most widely used genome editing technologies in recent years (20). As indicated by their name, clustered regularly interspaced short palindromic repeats (CRISPRs) are widely encoded in prokaryotic genomes (21). The unique characteristics of these CRISPRs and CRISPR-associated (Cas) proteins in bacterial and archaeal immunity have been rapidly identified recently (22). In the immunization process, a fragment of defined length from invading phage or plasmid DNA is incorporated into the 5 end of a CRISPR locus with a constant motif sequence. Accordingly, the periodic interspaced repeats of CRISPRs have been derived by continuous cycles of this immunization process. In the immunity process, an RNA originating from the immunized DNA is transcribed and processed and guides Cas protein(s) to its complementary sequence of exogenous DNA for cleavage and degradation. Harnessing different Cas proteins and RNAs involved in the immunization/immunity processes of different CRISPR-type families, various genome editing technologies have been established (20,23,24). Cas9 with double-stranded DNA cleavage activity from the type II CRISPR system has been widely used for targeted gene disruption and targeted fragment knock-in in various organisms including mammals. Similar to ZFNs and TALENs, nuclease-deficient Cas9 (dCas9) or mutant Cas9 nickase (nCas9) fused to effector proteins such as transcription factors, deaminases, and fluorescent proteins have been used for various applications such as gene silencing (25), activation (26), base editing (27), and chromosomal labeling (28).
Since periodically repeating DNA and protein sequences have diverse and important roles in biology, a simple and optimistic hypothesis could be proposed that new genome editing modules can be discovered from other periodic repeats in large-scale genomic resources. However, there is no universal software that captures various types of periodic repeats from large-scale genomic datasets in an unsupervised manner (Table 1). For example, RepeatMasker is one of the most commonly used tools to detect interspersed DNA repeats and low-complexity DNA sequences (29). However, this software screens only DNA sequences against a database of reported elements and does not evaluate repeat periodicity. Previous software programs developed for de novo searches of repetitive biomolecular sequences also have certain limitations. Tandem Repeat Finder is one of the first types of software to screen tandem and low-complexity DNA repeats without prior knowledge (30), but is incapable of capturing highly degenerated or interspaced DNA sequences or protein repeats. RECON (31) and RepeatScout (32) also screen only DNA sequences, focus only on interspersed repeats regardless of periodicity, and exclude tandem or low-complexity repeats. PRAP captures both tandem and interspersed repeats, but screens only DNA sequences (33). Although the recently developed software XSTREAM (34) and T-REKS (35) search for both tandem and highly degenerate repeats from DNA and protein sequences, both are ineffective at capturing interspersed or interspaced repeats including CRISPRs. With the recent interest in genome editing, several software packages such as CRISPRFinder (36), CRISPRDetect (37) and Anno-TALE (38) have been developed to capture genome-editingassociated sequences. However, such specialized software does not have the potential to discover novel genome editing modules.
The previously developed software focuses on limited types of repeat sequences for specific biological targets, but it seems that any software that combines the abilities of the previous software packages for any type of repetitive sequence would give an ambiguous and large set of sequences, which would require substantial effort for further curation and validation. However, none of the above-mentioned software screens repetitive sequences based on sequence periodicity that commonly appears in many significant biological processes. This could be a strong constraint in screening to obtain a set of biomolecular sequences with high potential for expanding our biological knowledge and developing new biotechnologies. Accordingly, we have been motivated to develop simple and fast software called SPADE (Search for Patterned DNA Elements) that globally captures such periodically repetitive biomolecular sequences in large genomic datasets mainly based on a simple evaluation of kmer periodicity.

Genomic resources
The GenBank files for the 7006 complete prokaryotic genomes (downloaded on 31 March 2017) and the human reference genome version GRCh38.p10 were downloaded from the NCBI RefSeq genomes FTP server (ftp://ftp.ncbi. nlm.nih.gov/genomes/refseq/).

SPADE screening phase 1: detection of highly repetitive regions
In SPADE, each entry sequence is first scanned by a sliding window of a given size w to roughly detect highly repetitive regions (HRRs). Let W i and k i be the sliding window at sequence position i of the entry sequence and its leftmost k-mer sequence, respectively. At every sliding window position i, the number of k i within W i (n i ) is cumulatively counted for every position in the left-most k i sequence region. In the same sliding window, 1 is also counted for every position of the other k i sequence regions. Let c i be the resulting cumulative k-mer score at position i of the entry sequence. After scanning the query sequence with the sliding window, cumulative k-mer peak areas (CKPAs) for which the peak heights are s or more are extracted. From the sequence regions for all of the CKPAs (c > 0), the broadest possible regions consisting of multiple gaps of size g or less were extracted as highly repetitive regions (HRRs). We adopted w = 1000, k = 10, s = 20, and g = 300 for nucleotide sequence and w = 300, k = 3, s = 6, and g = 50 for protein sequence as default parameters of the software.

SPADE screening phase 2: evaluation of periodicity
Let h be the size of a given HRR. For each HRR surrounding region of size h ± m, SPADE generates a position-period matrix (PPM) using a similar sliding window of size h. Let W i and k i be the sliding window at sequence position i and its left-most k-mer sequence, respectively. When multiple k i sequence regions are detected in W i , the number of k i regions (n i ) is cumulatively counted for all of the corresponding row-column cells of the first two k-mer regions, where row represents distance between two identical k-mers and column represents sequence position. When n i > 2, from the second k i sequence regions, this procedure is iteratively repeated except that the number added to each cell is 1. After scanning by the sliding window, the highest peak period d in the column sum distribution of the resulting PPM is identified. All values in a sub-PPM of rows from [d × 0.8] to [d × 1.2] are then added up and divided by that from 1 to half the column size of the PPM to produce the periodicity score. HRRs with periodicity scores of p and more are redefined as periodic repeat regions (PRRs). We set m = 1000 and P = 0.5 for nucleotide sequence and m = 300 and P = 0.3 for protein sequence as default parameters.

SPADE screening phase 3: identification of repetitive motifs
From each PRR with sequence period d, the k-mer sequence that has contributed the most to the sequence periodicity is extracted as k seed . When multiple k-mer sequences are extracted as the k-mers contributing the most, the left-most k-mer in the PRR is selected as k seed . Starting from all of the k seed sequences found in the PRR, SPADE obtains sequence fragments of size d. The extracted sequences are then aligned by multiple sequence alignment using MAFFT version 7.22 (39) to identify their consensus sequence motif. For each sequence position of the alignment result, the information content of appearing letters (b, bit) and the frequency of alignment gaps (f) are calculated using the Python WebLogo 3.6.0 package (40). After removing positions with f of more than q from the alignment result, letter consistency l of every position is calculated by b × f. The positions of the alignment result are then treated as circular since they are for periodic repeats and punctuated by removing the longest continuous nonconsensus region (l < u) of more than r letters. When this punctuation does not happen, the sequence alignment result is linearized as it was before.
To map the repeat motif to the PRR sequence, a representative sequence is obtained by taking the most frequent letter in each position of the alignment result. When the representative sequence is shorter than k-mer, the identical sequence regions are scanned in the PRR and annotated as repeat units. Otherwise, the representative sequence is mapped using BLAST+ version 2.6 (41) with the blastn-short (for nucleotide) or blastp-short (for protein) option and alignment length threshold of 50% to the query length or E-value of 0.01 or less. The hit regions in the PRR are then used to construct a sequence logo profile using Python WebLogo 3.6.0 package (40). From the sequence logo profile, a repeat motif sequence is generated by the most frequent letters, where a highest letter frequency of less than 60% is masked with '*'. q = 0.5, u = 0.8, and r = 5 were set as default parameters of the software.

Protein secondary structure prediction
For each visualized protein repeat motif sequence, the confidence score for ␣-helix structure or ␤-sheet structure was calculated using PSIPRED version 3.3 (42). For each PRR detected by SPADE, PSIPRED was initially used to predict all possible secondary structure motifs with the confidence score at each amino acid residue position. We then calculated the average confidence score for each motif at every position in the repeat sequence unit.

Evaluation of performance for detecting CRISPRs
From the entire periodic DNA repeats detected by SPADE, we extracted CRISPR candidates with interspace sizes of 25-60 bp and repeating periods of 58-81 bp. The interspace size parameters and the minimum threshold for the repeating period (interspace size plus repetitive sequence size) were set with reference to the CRISPRFinder screens for CRISPR candidates with interspace size being 25-60 bp and repetitive sequence size being 23-55 bp, but our maximum threshold for the repeating period was defined empirically based on the reported RefSeq CRISPRs. Region overlap agreement (ROA) between two given regions was calculated by dividing the size of the overlapping region by the combined size of the two regions. Recall and precision of the recapturing RefSeq CRISPRs were evaluated for each ROA threshold.

Evaluation of performance for detecting tandem protein repeats
From the 7006 prokaryotic genome resources, we screened the positive reference set (PRS) proteins for TALE, TPR, ANK repeat and WD40 repeat families using HMMER ). A total of 100 000 randomly picked prokaryotic proteins and the entire human proteome were screened for Pfam-A domain families version 31.0. Among those that do not have more than one copy of any Pfam domain with an E-value of less than 1.0e-10, we randomly selected 10 000 prokaryotic proteins and 10 000 human proteins as negative reference sets ProNRS10K and HuNRS10K. The performance of the software programs SPADE, XSTREAM and T-REKS was estimated using the recall of PRS proteins and the false positive rate (FPR) in ProNRS10K (for prokaryotic protein repeats) or HuNRS10K (for human protein repeats). The positive likelihood ratio (PLR) was calculated by dividing recall by FPR. Each software was used with its default parameters. Similar analysis was also performed by restricting the detected repeat unit sizes to within the range of expected sizes for different repeat families (34 ± 5 aa, 34 ± 5 aa, 33 ± 5 aa, 42 ± 5 aa, and 28 ± 5 aa for TALE, TPR, ANK repeats, WD40 repeats, and C2H2 ZNF, respectively). Note that, owing to the size filtering, FPRs varied for different repeat families, even when the same negative reference set was used. These measurements were also repeated with PRS proteins prepared using different criteria.

Overview of SPADE
We implemented SPADE to efficiently screen periodically repeating sequences as follows (Supplementary Figure S1). The software first automatically extracts multiple sequence entries from an input file (GenBank or FASTA format) and identifies the sequence type (DNA or protein) for each entry. Each entry sequence is scanned by a sliding window to count k-mers and highly repetitive regions are extracted.
The sequence periodicity of each highly repetitive region is then evaluated based on a position-period matrix that cumulatively plots the distance between the same neighboring k-mers and their sequence positions (see 'Materials and Methods'). From each periodic sequence region, the periodic sequence units are queried for a multiple alignment to identify repetitive motifs. A representative motif sequence is then aligned back to the entry sequence to annotate the periodically repeating units. Finally, the annotations for the detected periodic repeats are added to the input information and output in the GenBank format with options to visualize k-mer density, position-period matrix, repetitive unit loci with neighboring genes, and motif sequence logo for each periodic repeat.

Periodic repeats in a CRISPR-encoding genome
Using SPADE, we exhaustively searched for periodic DNA and protein sequences in the 7006 complete prokaryotic genomes that were available in the NCBI RefSeq database. The default parameter set was used for the entire analysis of this study. In the Streptococcus thermophilus LMD-9 genome, 7 periodic DNA repeats and 27 periodic protein repeats were detected, including 2 previously annotated CRISPR loci ( Figure 1A, Supplementary Table S1). The repeat periods of the annotated CRISPRs were both 66 bp and their detected repeat motif sequences were identical to the reported motifs ( Figure 1B). We also found an unannotated interspaced repeat region containing four repeats with a period of 72 bp, in which the repeat motif and interspace sequences were all 36 bp long ( Figure 1C). While type II-A Cas genes were found in the neighboring regions of the reported CRISPRs, a type III-A Cas gene cluster was found in the adjacent region of the unannotated repeat, suggesting a functional type III-A CRISPR system in this genome. The other periodic DNA repeats were all short tandem repeats with a period size of 1-7 bp that were commonly found in prokaryotic genomes ( Figure 1D). Among the 27 periodic protein repeats, 24 were short tandem repeats with periodicity of 10 aa or less. The other three included a peptidoglycan-binding protein (three repeats with a 17aa period) and a subtilisin-like serine protease (three repeats with a 32-aa period), both of which were annotated to involve protein repeats, and a nucleotide exchange factor (four repeats with a 14-aa period), which was annotated to involve two ␣-helices.

Performance in detecting CRISPRs
We then measured the performance of SPADE in detecting CRISPRs, the annotation criteria of which are standardized in the NCBI prokaryotic genome annotation pipeline (43). From the entire 161 465 periodic DNA repeats detected in the 7006 prokaryotic genomes, we obtained 8168 genomic regions with a repeat period size and interspace size of 58-81 and 25-60 bp, respectively (Supplementary Table S2). These parameters were partly derived from CRISPRFinder, the most commonly used tool for CRISPR annotations in recent genomic resources (36,44,45), and partly defined empirically based on the reported CRISPRs (see 'Materials and Methods'). We confirmed that the distribution of period-interspace size combinations for the defined parameter space had good agreement with that for the 6354 reported CRISPRs in the RefSeq database ( Figure 1D and E). We then compared the performance of SPADE with the commonly used CRISPR annotation software tools CRISPRFinder and CRISPRDetect in capturing RefSeq CRISPRs. In the same genomic datasets, 8033 regions were detected by CRISPRFinder and 5924 regions were detected by CRISPRDetect (Supplementary Table S2). Precision and recall were decreased along with region overlap agreement (ROA) with reported RefSeq CRISPR regions for all of the software tools but recalls by SPADE were the highest for ROA of up to 98% where the recalls by CRISPRDetect were markedly lower than those by the other tools ( Figure  1F). On the other hand, CRISPRDetect outperformed pre-  Figure 1F).
At 50% ROA, SPADE, CRISPRFinder and CRISPRDetect captured 6181, 6093 and 5548 RefSeq CRISPR regions, respectively, where 5420 were captured by all of the software tools ( Figure 1G). We defined 5548 RefSeq CRISPR regions captured by CRISPRDetect with a minimized false positive rate as a high confidence gold standard CRISPR set. Amongst this set, 5520 and 5445 regions were captured by SPADE and CRISPRFinder, respectively. Given that precisions and recalls of CRISPRFinder were higher than those by SPADE for ROA of more than 98%, we concluded that SPADE was slightly better than CRISPRFinder at roughly capturing CRISPRs, but not at the single-base resolution. In sum, although SPADE was not specifically designed for CRISPR annotation, its performance for capturing CRISPRs with simple size thresholds was at least comparable to the most commonly used CRISPR prediction software tools.

Periodic repeats in a TALE-encoding genome
In the Xanthomonas oryzae pv. oryzae (Xoo) PXO83 genome encoding TALE genes and TALE pseudogenes (38), 49 DNA repeats and 194 protein repeats were detected by SPADE (Figure 2A and B, Supplementary Table S3). All of the reported TALEs were recaptured with a repeating period of 34 aa and variable residues at the 12th and 13th amino acid residues and two ␣-helices in each repeat unit, which were all consistent with the reported features of TALE. We also detected an annotated large CRISPR locus where a highly constant motif of 31 bp was repeated periodically 86 times each with an interspace sequence of around 34 bp ( Figure 2C).
Among the other 46 periodic DNA repeats, 40 were short tandem repeats with a period of 10 bp or less, including 25 heptamer repeats that were previously suggested to contribute to phase variation in the Xanthomonas genus (46). Three short tandem DNA repeats were found in intergenic regions, one with a period of 12 bp and two with a period of 14 bp. Another short tandem DNA repeat region was found in the middle of an ABC transporter-encoding gene with a period of 16 bp, which is relatively prime to 3, the protein coding frame size ( Figure 3A), and another longer sequence with a period of 60 bp was also found to encode a hypothetical gene in less than half of its region ( Figure 3B). Furthermore, we found a large periodic DNA region from the genomic position of 3 559 997 to 3 563 142 (3144 bp long) with an average period of ∼787 bp (Supplementary Figure S2). Following a transposase-encoding gene, this region involved three different hypothetical genes, each of which was in a different repeat unit. Interestingly, all of the three repeats partially overlapping with proteincoding regions were found to be widely conserved in the Xanthomonas genus with different numbers of repeats, but the coding gene architecture had markedly diverged evolutionarily ( Figure 3C-E), indicating that phase variations of protein-coding patterns for these regions rapidly occurred after speciation by genomic contraction and expansion via the repetitive sequences.
Except for TALEs, ∼88% (156 out of 178) of protein repeats were composed of short tandem repeats with a repeat unit size of 10 aa or less (Supplementary Table S3). The other repetitive proteins included three chemotaxisassociated proteins with different periods of 27, 46 and 90 aa, a DNA topoisomerase I, a TolB-like protein known to involve non-WD40 ␤-propellers, and transporters and a hypothetical protein involving six repeats with a large unit size of 215 aa. Notably, another type III secretion system effector protein of the Xanthomonas host infection process was found to have repetitive peptide units, suggesting another function of pathogenic periodic protein structure in hijacking the host plant system (Supplementary Figure S3).

Performance in detecting TALEs and C2H2 ZFNs
As C2H2 ZNFs are the most widely used transcription factors in the human genome, we also examined whether SPADE can capture human C2H2 ZNFs. When a C2H2 ZNF encoded on human chromosome 7p22.1 was scanned by SPADE, 20 degenerative repeats of ∼28 aa were detected with two cysteine and two histidine residues conserved at specific positions, like typically reported C2H2 ZNFs (Figure 4). We then assessed the performance of SPADE in detecting TALEs and C2H2 ZNFs. Using the protein domain search software HMMER, we obtained positive reference sets (PRSs) for TALE and human C2H2 ZNF from the prokaryotic genomic dataset and the human proteome, respectively, so each PRS protein contained three or more of the corresponding Pfam motifs (see 'Materials and Methods'). We also prepared 10 000 prokaryotic proteins and 10 000 human proteins that did not have any Pfam motif more than once as negative reference sets (NRSs) ProNRS10K and HuNRS10K, respectively (see 'Materials and Methods'). Using SPADE, repetitive sequences of any period were detected in 328 out of 331 TALE PRS proteins (99.1%) and 3079 out of 4084 human C2H2 ZNF PRS proteins (75.4%), while 192 ProNRS10K proteins (1.9%) and 1269 HuNRS10K proteins (12.7%) were positive (Figure 5A, Supplementary Table S4). When the detected positives were filtered by maximum repeat unit size per protein (maxRUSPP) to be within ±5 aa from the expected average repeat unit size (34 aa for TALE and 28 aa for C2H2 ZNF), the recall of TALE PRS stayed the same (99.1%) and the recall of human C2H2 ZNF PRS was 58.9%, while the false positive rate (FPR) of TALE estimated using ProNRS10K and the FPR of C2H2 ZNF estimated using HuNRS10K were greatly decreased to 0.03% and 0.21%, respectively ( Figure 5B and C). This simple size limitation improved positive likelihood ratios (PLRs) of the prediction from 51.6 to 3303.1 (64.0-fold) for TALE and from 5.9 to 280.7 for human C2H2 ZNF (47.2-fold).
Comparison with other software capturing tandem protein repeats SPADE successfully detected the other degenerate tandem protein repeats widely spread in prokaryotes, including TPRs, ANK repeats, and WD40 repeats ( Figure 5D-F). The secondary structure prediction of these degenerated repeats also properly captured their reported structural  motifs. In each of the repeat sequence motifs identified for a TPR and an ANK repeat, both of which have been reported to have helix-turn-helix structures, we observed two ␣-helical loops ( Figure 5D and E). Four ␤-strands were also captured in a repeat sequence motif of WD40, consistent with its ␤-propeller structure ( Figure 5F). As XSTREAM (34) and T-REKS (35) have been widely used to explore tandem protein repeats in an unsupervised manner in recent studies (47,48), we next performed a benchmark comparison of SPADE, XSTREAM, and T-REKS in detecting TPRs, ANK repeats, and WD40 repeats, in addition to TALEs and human C2H2 ZNFs. For TPRs, ANK repeats, and WD40 repeats, PRSs were prepared as described above for TALE. ProNRS10K and HuNRS10K were again used as NRSs for detecting repeats in prokaryotic and human protein families, respectively.
T-REKS performed the best in recall for detecting repetitive sequences regardless of repeat unit size, except for WD40, in which SPADE performed the best (Figure 5A, Supplementary Table S4). However, T-REKS also demonstrated the highest FPRs in both ProNRS10K and HuNRS10K datasets. When the overall prediction performance was estimated by PLR, SPADE performed the best in every repeat type (between 1.02-fold and 3.39-fold compared with the second-best software XSTREAM for all repeat types). We also found that the maxRUSPPs detected by SPADE were distributed with peaks at 34, 33 and 42 aa for TPR, ANK repeats, and WD40 repeats, respectively, all of which were the reported typical unit sizes for these protein repeats ( Figure 5B). This was not the case for all of the repeats detected by XSTREAM and T-REKS. XSTREAM captured wider ranges of repeat unit sizes for every repeat type and T-REKS tended to capture shorter tandem repeats for the subpopulation of positive reference proteins for TPRs, ANK repeats, and WD40 repeats. Filtering the detected positives by maxRUSPP to be within ±5 aa from the expected average repeat sizes, the recall performance of SPADE was the best for all repeat types, whereas the FPRs of the three software packages were all minimized to below 0.005 in all of the repeat types ( Figure 5C, Supplementary  Table S4). (Note that the performances could not be compared using PLR as many FPRs for different protein families were zero.) These observations were maintained when the positive reference protein sets were prepared differently (Supplementary Figure S4).

DISCUSSION
No software program have been developed that can universally screen for periodic DNA and protein repeats; the only available software tools are those that screen for reported motifs or certain types of periodic repeats. Nevertheless, the performance of SPADE capturing CRISPRs was on par with the commonly used CRISPR prediction software tools and outperformed XSTREAM and T-REKS in the sensitivity for capturing various tandem protein repeats, regardless of the degree of consensus in the repeat unit motifs. SPADE also captured TALEs and ZNFs in a highly specific and unsupervised manner, indicating its potential to contribute to the discovery of new genome editing modules from large genomic and/or metagenomic resources. This is supported by the fact that we found that a non-TALE     Figure S3). We also captured bacterial homologs of pentatricopeptide repeats (PPRs) that are involved in translational regulation in plants (Supplementary Figure S5). As the binding code of PPR to RNA has re-cently been deciphered, it has been suggested as a potential programmable RNA editing and regulating module (49).
The majority of the periodic repeats detected in the 7006 prokaryotic genomes still need further investigation. We detected many short tandem DNA and proteins repeats. In particular, tandem heptamer DNA repeats were the most abundant in intergenic regions of a wide range of prokaryotic species ( Figure 1D). However, there has been no clear clue about the function of this globally existing prime number periodicity in genomic DNA. We also found various interspaced repeats that had clear sequence periodicities with no CRISPR annotation or neighboring Cas gene. They included many tRNA operons in various prokaryotes, as reported previously (Supplementary Figure S6), but the others remain to be explored. Genomic expansion and contraction have been thought to occur at the tandem repeat sequences, leading to phase variation. Even after excluding corresponding protein repeats, the repeat periods of both tandem and interspaced DNA repeats showed particular abundance for these in multiples of three. Furthermore, some genes were encoded in part of a repeat unit of a large tandem repeat region (Supplementary Figure S2). As seen in the Xanthomonas genus ( Figure 3E), these findings suggest the roles of tandem repeats in de novo gene birth or gene death. We also found many tandem DNA repeats within (or partially within) protein-coding regions, some of which  repeats. (B) Distribution of maximum repeat unit sizes per protein (maxRUSPPs) detected by each software for each protein family. Each bold black bar and each gray box denote the reported typical repeat unit size for the corresponding protein category and the ±5 aa range from the reported repeat unit size, respectively. (C) Recalls and FPRs of the different software after filtering maxRUSPPs of the detected positives to be within ±5 aa from the expected repeat unit size. From the maxRUSPP distributions of ProNRS10K (negative control for prokaryotic protein repeat families) and HuNRS10K (negative control for ZNFs), FPRs for different protein repeats of different expected repeat unit sizes were estimated. (D-F) Example protein repeats detected by SPADE for TPR (D), ANK repeat (E), and WD40 repeat (F). The heat map under each repeat motif sequence logo represents the confidence scores for ␣-helical structure (red heat map) or ␤-sheet structure (blue heat map) at each amino acid residue position.
were indicated to have contributed to functional phase variation of protein-coding patterns ( Figure 3A-D).
The default parameter set of SPADE robustly captured most of the important biological sequences tested in this study with higher precision than did the other software. We further analyzed various kinds of simulated repeats using SPADE and confirmed that the default k-mer parameters for DNAs and proteins performed the best and precisely captured the periodicities of various degenerate tandem and interspaced repeats with up to ∼10% and ∼30% sequence noise for DNA and proteins, respectively, regardless of the repeat unit size and interspace size (Supplementary Figure  S7). SPADE was implemented using Python and can be executed with MAFFT and BLAST on MacOS X and Linux operating systems, and on Windows Subsystem for Linux (WSL) of Windows 10. It automatically detects the input file and sequence types and outputs results in the Gen-Bank file format, which can be further input into other software programs, with various visualizations as represented in the figures. Accordingly, here we propose that SPADE is fast and user-friendly software based on a simple algo- rithm to globally capture periodic biomolecular sequences. Although we mainly focused on measuring the performance of this software predominantly using prokaryotic genomes in this study, further wide-ranging investigations of these periodically repeating sequences together with screening of eukaryotic and metagenomic resources could lead to the discovery of new biological events and genome editing tools.

DATA AVAILABILITY
SPADE is open-source software under GNU General Public License v3.0, which is available at https://github.com/ yachielab/SPADE.