BACKGROUND
genomic dna sequences contain numerous repetitive motifs or nucleotide sequence patterns that play various roles in the physiology of the cell, including regulation of gene expression, maintenance of dna and its structure in the cell, recombination and other forms of mutation, and even recognition of dna fragments for uptake from the extracellular space .
repetitive motifs in dna sequences can also encode functional elements in proteins and rna.
consequently, a large amount of effort has been devoted to development of computational methods for detection of sequence motifs in dna and protein sequences.
there are two major types of motif finding algorithms, namely supervised and unsupervised motif finding algorithms.
the former methods require a sample of known occurrences of the motif and utilize this information in the search for additional motif occurrences in the analyzed sequence or sequences.
the unsupervised methods, sometimes referred to as ab initio approach, do not require any prior knowledge about the motif sequences and detect novel sequence motifs that satisfy specified criteria .
in this article, we focus on unsupervised motif search.
the unsupervised motif finding algorithms can be further classified into two major groups: 1) word-based methods that mostly rely on exhaustive enumeration, i.e., counting and comparing oligonucleotide frequencies and 2) probabilistic sequence models where the model parameters are estimated from sequences.
extensive work has been done on transcription factor binding site  prediction during the past decades, driven by the obvious importance of these regulatory motifs in the organism’s physiology.
based on the type of dna sequence information used by the tfbs finding algorithm, the methods could be classified into three major classes: 1) methods that use promoter sequences from coregulated genes from a single genome  <cit> , 2) methods that use orthologous promoter sequences of a single gene from multiple species  and 3) methods combining 1) and 2)  <cit> .
as a unified portal for online discovery and analysis of sequence motifs, the meme suite web server provides various tools in finding motifs representing features such as dna binding sites and protein interaction domains  <cit> .
while tfbs receive most attention among sequence motifs the repetitive sequence motifs in genomic dna can have many other functions.
we aim to expand the range of types of sequence motifs detected by motif-finding methods by searching for spaced motifs in complete genomes, which could arise, among other mechanism, from repetitive patterns in chromosome structure.
the concept of searching for spaced sequence motifs is not new but the previous such methods generally aimed to detect tfbs, where spacing of the conserved segments of the motif is determined by the geometry of the dna-protein interaction and generally does not exceed  <dig> or 7 bp .
a more general approach implemented in helicis  <cit>  allows detection of co-localized periodically spaced motifs, such as binding sites for multiple transcription factors.
however, these methods are specifically designed for detection of tfbs in a collection of regulatory regions and are not suitable for scanning complete genomes.
in contrast, our methodology is aimed at detection of dna sequence motifs that can have diverse physiological functions, including regulatory motifs, motifs involved in dna, rna, and protein interactions, and motifs related to structural organization of dna and/or the encoded proteins.
our approach is different from previous methods in that it searches for pairs of motifs that occur at a particular distance from each other significantly more often than expected without requiring that either of the motifs on its own occurs more frequently than expected.
the goal is to detect paired motifs where the distance between the two motifs is under selective constraint but the specific sequence of the spacer is variable.
methods
the motif-finding algorithm
the goal of the algorithm is to find motif pairs that have a preferred distance of separation.
that is, they occur at a specific distance apart more often than expected by chance.
the algorithm consists of the following steps:scan the genome and count the number of times that the pairs of motifs a and b  occur at every spacer distance between  <dig> and  <dig> .
record the total number as n;
identify the distance  at which the pair has the maximum number of occurrences, record the maximum number as fmax;
test whether the maximum number of occurrences fmax is large enough to claim significance of overrepresentation by surpassing some appropriate cutoff value fcut;
repeat steps 1) – 3) for all  <dig> possible pentamer pairs; record all significant pairs with fmax > fcut;
for each significant pair, align the spacers  pair-wise.
reduce fmax by  <dig> for each pair of aligned spacers with nucleotide identity greater than 70% and remove one of the two spacers from further comparisons with other spacers;
identify and record motif pairs whose fmax is still above fcut.
record if the inverse complement is also significant with the same distance;
for motif pairs that satisfy the criteria above, investigate their distribution in the genome using pattern locator  <cit> .
step 1) consists of recording the total number of times a pair of pentamers a and b is found at the exact mutual distance d. the distance d is counted from the beginning of a to the beginning of b. pairs with exactly the same spacers are counted only once .
the motivation for discounting motif pairs with identical or very similar spacers is that they likely arise from sequence duplication of the whole segment containing the motif pair and the spacer, which is not our primary goal; we aim to detect motif pairs where the distance between the motifs arises from selective constraints rather than as a result of sequence duplication events.
the main reason we limit the motif spacer distance to 5–89 bp is due to growing computational complexity when ranges of paramters investigated are too large.
similar reasons apply to limiting the investigated motifs to pentamers; the number of possible pairs of hexamers is around  <dig> million , which is prohibitive when the methodology is to be applied to hundreds of complete genomes as in this study.
another important motivation to use pentamers is related to how prokaryotic dna-binding proteins, particularly transcription factors, interact with dna.
many prokaryotic transcription factors interact with dna in a sequence-specific manner via two alpha-helices that fit into adjacent turns of the major groove in the dna double helix  <cit> .
consequently, although the whole binding site is often ~ <dig> nucleotides long, it usually consists of two conserved segments of ~ <dig> nucleotides where the protein directly interacts with the dna bases, separated by a variable gap  <cit> .
after finding the most abundant distance d for each motif pair in step 2), statistical significance is assessed in step 3).
for a given pair of motifs, if no specific distance is preferred, one would expect a uniform distribution of spacer distances within the whole range  <dig> to 89 bp.
the problem of numbers of motif pairs found at different distances is analogous to the classical balls-in-urns problem where motif pairs represent the balls and different distances  represent the urns.
suppose n balls are allocated to k urns at random, each ball being equally likely to fall into each urn, independently of any other balls.
let fmax denote the largest number of balls in any urn.
this situation can be modeled with the uniform multinomial distribution in which there are n trials, each with k possible outcomes of equal probability.
david and barton  <cit>  showed that p can be computed as  <dig> p=n!kn×thecoefficientofλnink
subsequently, the p-value p = p =  <dig> - p.
the analogy in dna sequences is straightforward with n being the total number of motif pair occurrences with spacer distance between  <dig> and 89 bp, k being the number of different distances  and fmax being the maximum number of occurrences of a specific distance within the range.
a bonferroni correction was applied to the p-value cutoff in order to account for multiple hypothesis testing :
p-value-cutoff =  <dig>  / number of motif pairs =  <dig>  / 410 ≈ 10− <dig>
testing over multiple distances is accounted for in the uniform multinomial.
sympy, a python library for symbolic mathematics, was utilized to perform the power expansion in formula  and to extract the coefficient needed to compute p. however, using sympy to obtain the exact tail probability is computationally expensive, especially for large parameters of n and k. hence, we used an approximation  <cit>  to calculate the 100% percentile of fmax with the following equation  <dig> fα≈nk+nlogkk−loglog11−α+ <dig> n2klogk
in the formula above, fα is the 100% percentile of fmax; α is the desired type i error rate; n is the total number of occurrences of the motif pair at hand with distance between  <dig> and  <dig>  and k is the total number of distances investigated .
note that the above approximation is based on an asymptotic theorem under certain regularities  <cit> .
to verify that the approximation is sufficiently accurate with our set of parameters, comparisons between approximate and exact fcut values for different values of k were made using a fixed p-value 10− <dig> .
when k is relatively large , the equation  approximates the true values with reasonable accuracy as n gets large; however, when k is relatively small , the equation tends to overestimate the true fcut values.
since the k used in this paper is  <dig>  we decided to use exact cutoff values for n < =  <dig> while the approximated values are used for n >  <dig> fig.
<dig> comparison between approximate and exact p-value cutoffs.n is the total number of occurrences of some motif pair; fcut indicates the cutoff value corresponds to the selected p value threshold.
solid and dashed lines represent approximate values calculated by formula  and exact values computed with sympy, respectively
in step 4], the steps  <dig> to  <dig> are repeated for all possible pentamer pairs.
after all pentamer pairs are tested, the most promising candidates are selected for spacer alignment.
in step 1], only identical spacers could be recognized and eliminated but our goal is to discount all motif pairs which likely arose from recent duplication of the whole segment including the motif pair and the spacer.
consequently, sequence alignments are performed with nw-align, which is a simple and robust alignment program based on the standard needleman-wunsch algorithm  <cit> , to eliminate potential duplicated spacers.
a simple ‘greedy’ algorithm is applied to eliminate duplicated spacers.
the first pair of spacers is aligned.
if the nucleotide identity is ≥70%, then fmax is reduced by one and one of the spacers is excluded from further alignments.
additional pairs of spacers are subsequently aligned and excluded if the nucleotide identity exceeds 70%.
this is continued until no spacers are more similar than the 70% cutoff.
the value 70% was chosen as the lowest cutoff that gives a high confidence the sequence similarity is not coincidental; for example, the probability that two random sequences of ten nucleotides would match at seven positions  is less than  <dig>  and declines with increasing sequence length.
we also performed tests with 65% and 75% cutoffs, which suggested that the cutoff 65% was too conservative whereas majority pentamer pairs identified with the 75% cutoff overlapped motifs detected also with the 70% cutoff .
after this step, a motif pair whose fmax is still above fcut and whose inverse complement is also significant with the same distance is selected for further analysis.
the reason for requiring a motif pair’s inverse complement also to be significant is that dispersed sequence motifs, which are the primary target of our search, are expected to be distributed in both dna strands; consequently, both versions of the motif pair—direct and reverse—are expected to be over-represented in the genome.
the subsequent analysis of the significant motif pairs is performed with pattern locator  <cit> , which provides information on distribution of the motif pairs with respect to adjacent genes.
such information can be helpful in generating hypotheses about the motif pair’s possible function.
pattern locator also reports the percentage of copies of the motif pair found in genes and intergenic regions , which allows easy identification of motifs that are most likely to have regulatory functions  and motifs that could reflect conserved patterns in amino acid sequences of the encoded proteins.
dna sequences
annotated nucleotide sequences of complete prokaryotic genomes were downloaded from the ncbi ftp server  <cit> .
we randomly selected only one genome per species when multiple strains of the same species were available.
the final dataset included  <dig> genomes .
for genomes consisting of multiple chromosomes the analysis was performed on all chromosomes.
RESULTS
the  <dig> prokaryotic genomes, each from a separate species, were scanned with the algorithm described above.
in the  <dig> genomes,  <dig> motif pairs were identified as significant by our criteria .
the summary statistics for selected characteristics of the significant motif pairs are shown in table  <dig> table  <dig> summary statistics of  <dig> motif pairs identified in  <dig> genomes
the summary statistics is per motif.
meaning of abbreviations in the table: d, spacer length of a motif pair; initial_n, number of copies of the motif pair before alignment; reduced_n, number of copies after alignment and elimination of duplicate spacers; difference, the difference between initial_n and reduced_n; cutoff, significance cutoff ; reduced_n/cutoff, the ratio of reduced_n and cutoff, indicating the relative significance for the motif pair; gene, the percentage of each motif pair occurrences found in genes; intergenic, the percentage of each motif pair occurrences that are in intergenic regions; overlap, the percentage of each motif pair occurrences that overlap with a gene start or end.
these percentages are calculated as follows: for any given significant motif, we run a query with pattern locator, which gives the percentage of the motif occurrences that fall in gene, intergenic region or overlap with gene starts or ends.
the quantiles in the table are for these percentages over all significant motifs
distributions of some key statistics are further demonstrated in figs.
<dig> and  <dig>  the right panel in fig.
<dig> indicates that median percentage of significant pairs found in protein- coding regions is around 80%, which is close to the fraction of protein-coding region in a typical prokaryotic genome; in other words, the significant motif pairs do not show a strong overall bias for protein-coding regions or intergenic regions.
in fig.
<dig>  the first two plots show that the distributions of motif pairs are bi-modal in both gene and intergenic region.
<dig> percent of significant motif pairs occurred in genes more than 90% of the time and  <dig> percent occurred in genes less than 10% of the time, whereas fewer motif pairs are distributed approximately evenly among both genes and intergenic regions.
in other words, majority of the significant motifs show a strong preference for either genes or intergenic regions.
the same conclusion can be reached upon inspection of the raw data in additional file 3: table s <dig> fig.
<dig> distribution of key motif statistics.
frequency is number of copies of each significant motif pair; percentage is the percentage of significant motif pairs found in corresponding region
fig.
<dig> distribution of motif pairs with respect to gene and inter-genic region.
in each plot, the horizontal axis indicates the percentage of occurrences of significant pairs found in corresponding region
the distribution of the spacer length for the motif pairs is shown in fig.
<dig>  there are two peaks—one is around spacer lengths  <dig> to 20 bp and the other one is around  <dig> to 50 bp.
the reason for the first peak is partly due to a conserved palindromic motif pair , which is widely distributed among a variety of genes of mycobacterium species.
in total, approximately  <dig> motif pairs conforming to this consensus were identified in the analyzed mycobacterium genomes out of  <dig> motif pairs with spacer length  <dig> to 20 bp found in all analyzed genomes.
the peak at 40–50 bp reflects wide distribution of clustered regularly interspaced short palindromic repeats  among the analyzed genomes.
crisprs are arrays of ~30 bp perfect repeats separated by ~30–40 bp variable spacers, which give rise to many significantly overrepresented motif pairs with this spacer size  <cit>  .fig.
<dig> spacer length distribution of candidate pairs
another noteworthy statistic in table  <dig> is the difference between the initial number of copies of a significant motif pair  and the number of copies remaining after elimination of duplicate spacers .
the mean value of initial_n is  <dig>  whereas reduced_n has mean  <dig> , indicating that more than 25% of occurrences of all significant motif pairs in the genomes may have arisen from sequence duplication events.
this confirms the importance of removing the “false positives” resulting from sequence duplications when searching for candidate motif pairs whose specific relative positioning is maintain by selective constraints.
all significant motif pairs found in the analyzed genomes are reported in additional file 3: table s <dig>  in total,  <dig> genomes have at least one significant motif pair identified.
among those, rhodoferax ferrireducens t <dig> has the highest number of pairs detected,  <dig> .
thirty genomes  were selected for in-depth investigation of the identified significant motif pairs.
some of the significant motif pairs were found to be related to previously known motifs such as crispr or shine-dalgarno sequence whereas others appear to be potentially novel, previously uncharacterized sequence motifs.
majority of the identified motif pairs can be classified in one of the five categories, namely crispr-related, rho-independent transcription terminators, trna-related motifs, motifs containing the shine-dalgarno sequence, and protein-related motifs.
crispr-related motif pairs
crisprs are segments of prokaryotic dna containing periodic repetitions of ~ <dig> nucleotides.
each repetition is followed by short segments of “spacer dna”, typically 30– <dig> nucleotides long.
the spacers are short fragments of nucleotide sequence from phages or plasmids to which the bacterium had been exposed in the past  <cit> .
these fragments in conjunction with the cas proteins are subsequently used to recognize and destroy the dna of the same phages or plasmids in future encounters.
the crispr/cas system is a form of prokaryotic immune system that confers resistance to foreign genetic elements such as plasmids and phages by recognizing and digesting them  <cit> .
crisprs are found in approximately 40% of sequenced bacterial genomes and 90% of sequenced archaea  <cit> .
given its properties of regularly interspaced repeats with spacers of variable sequence but constant length, it’s not surprising that parts of crisprs are among the paired motifs detected by the algorithm.
table  <dig> lists significant motif pairs found in the  <dig> genomes subjected to in-depth analysis that are related to crispr.
results were verified by comparison with crisprdb, a comprehensive database of known crispr occurrences  <cit> .
candidatus desulforudis audaxviator mp104c contains two types of crispr sequences, gtttcaatccctcgtaggtaggctggaaac and ctttcagtccccttttcgtcgggtcggtcgctgaaac.
the first four patterns listed in table  <dig> for this genome correspond to the former whereas the last three reflect the latter crispr sequence.
the patterns consist of two pentamers near the ends of the crispr sequence and the spacers, which are variable in sequence bur have a constant length, are responsible for the fixed distance between the pentamers.
for patterns detected in clostridium thermocellum atcc  <dig>  the first two arise from repeats gtttcaattcctcataggtacgataaaaac and the rest are due to gttttatcgtacctatgaggaattgaaac.
all patterns found in corynebacterium aurimucosum atcc  <dig> are due to the crispr sequence gtgctccccgcgtaagcggggatgagcc.
similarly, all patterns detected in the other two corynebacterium species are due to the same consensus ggctcatccccgctggcgcggggagcac.table  <dig> significant motif pairs related to crispr repeats
gaaac35ggttt
aaact35gtttc
gaaac36gtttc
tgaaa39tttca
ctgaa40tttca
gtacg52tcctc
atgag53tcgta
tgagg53cgtac
atgag54cgtac
accta55tatcg
cctat56tcgta
accta56atcgt
cgggg44tcccc
cgggg45ccccg
gggga45cccgc
cgggg44tcccc
cgggg45ccccg
gggga45cccgc
cgggg44tcccc
gggga45cccgc
pattern x1x2x3x4x5dy1y2y3y4y <dig> denotes 5-mer x1x2x3x4x <dig> is followed by 5-mer y1y2y3y4y <dig> with d nucleotides apart; for non-palindromic patterns, only sequence in one dna strand is listed
note that our software is not intended for detection of crispr sequences and the detection of crispr-related sequence pairs are a byproduct of the characteristic crispr structure.
our method is designed for detection of dispersed sequence motifs, whereas crisprs consist of repeats clustered at a small number of loci in the genome.
this explains why crispr-related motif pairs were not detected in many genomes with crispr.
the main reason is the requirement that both direct and reverse complement versions of a motif pair independently satisfy the criteria for statistical significance .
many genomes contain a single crispr locus, which means that all crispr sequences are in the same orientation with respect to the direct and reverse strands.
such cripsr loci can only be detected if they contain a palindromic pair of pentamers that satisfies all other criteria.
rho-independent transcription terminators
a transcription terminator is a functional sequence located at the 3′ end of genes that mediates the transcriptional termination.
in prokaryotes, there are two types of terminators, namely rho-dependent and rho-independent terminators.
the former require rho factor to terminate the transcription process while the latter one forms a self-annealing hairpin structure to serve the purpose.
the terminator sequence generally contains a gc-rich region of dyad symmetry followed by a short poly-t tract  <cit> .
motif pairs related to transcription terminators were found in  <dig> out of the  <dig> genomes .
the results were verified by transtermhp, which is a specialized software to find rho-independent transcription terminators in bacterial genomes  <cit> .
in actinobacillus pleuropneumoniae the pattern is a high-g + c palindrome; in haemophilus influenzae the pattern detected by our method is not a palindrome, but instead includes a poly-t tract.
both features in the two patterns identified correspond to the standard structure of rho-independent terminators, which consist of a gc-rich palindrome followed by a short poly-t tract.table  <dig> motif pairs that are part of transcription terminators
caaaa15accgc
ggtct15ttgca
tgcaa15cgacc
tgcaa16aaccg
gcaaa16accgc
ttgca16agacc
gcggt16tttac
gcaaa17ccgct
cggtc17tgcaa
tttgc18gaccg
aagcg19ttgca
gcggt19gcaaa
ttttg19gaccg
cggtc38gaccg
ggtca11gtttt
cggtc12gtttt
gcggt12agttt
gcggt12tgttt
ggtca12ttttt
cggtt12gtttt
tgcgg13agttt
cggtc13ttttt
gcggt13gtttt
tgcgg13tgttt
cggta13ttttt
cggtt13ttttt
gcggt13atttt
gcggt14ttttt
gtgcg14tgttt
tgcgg14gtttt
tgcgg15ttttt
agtgc15tgttt
gtgcg15gtttt
gcggt15ttttg
agtgc16gtttt
aaaca16cactt
tgcgg16ttttg
the fact that transcription terminator-related motif pairs were detected by our program is actually unexpected because the terminators are characterized by the stem-loop secondary structure but not a specific repeated sequence.
this is also the reason why we did not detect transcription terminator-related motif pairs in other genomes.
however, the h. influenzae genome contains about  <dig> copies of the uptake signal sequence with the consensus aagtgcggt and terminators in h. influenzae are often composed of two copies of the uss in the opposite orientations that form the stem of the stem-loop structure  <cit> .
that leads to prevalence of this specific sequence in the terminator stem-loop structures and it is the reason why we are detecting terminators in h. influenzae.
analysis of the a. pleuropneumoniae genome with the aimie software  <cit>  revealed an abundant motif aagcggt, which is similar to the uss in h. influenzae and may have an analogous function  <cit> .
like the uss in h. influenzae, the a. pleuropneumoniae motif frequently occurs in dyad pairs found at the 3′ ends of genes, which probably serve as transcription terminators.
trna-related motifs
trna-related motifs were detected in four genomes –  <dig> of them are inside trna genes and correspond to parts of a conserved pattern tagagc27ggttcg near the 5′ end of the gene .
parts of this pattern are highly conserved in trna genes from diverse species  <cit> .
however, the motif pairs in a. butzleri, s. muelleri, and n. profundicola just exceed the significance cutoff of our method and they are not detected in other genomes because the number of copies is not statistically significant.table  <dig> trna-related motifs identified
the pattern in methanopyrus kandleri av <dig> differs from the other three in that it resides outside of trna genes, specifically about  <dig> base pairs downstream of the trna genes.
other trna genes  have similar patterns at the same position .
the consensus sequence of the motif is shown in fig.
<dig> fig.
<dig> sequence logo for motif pair downstream of trna genes in m. kandleri
the sequence is palindromic with positions 2– <dig> and 15– <dig> being the most conserved.
unlike the other trna-related motif pairs, this motif appears to be specific to m. kandleri and we did not find similar motifs in other genomes.
the methanogenic archaeon m. kandleri grows at temperatures of 84–110 °c and genome analyses revealed strategies to adapt to these harsh conditions.
among other adaptations, m. kandleri uses a unique trna c-to-u editing mechanism at base  <dig> for  <dig> different trna species  <cit> .
we speculate that the palindrome discovered by our software might serve as a target signal for the 3′ maturation enzyme and its presence specifically in m. kandleri may be related to the unique mechanism of trna maturation employed by this organism.
motifs related to shine-dalgarno sequence
located around 8 bp upstream of the start codon, shine-dalgarno  sequence serves as ribosome binding site in mrna of prokaryotic genomes  <cit> .
with a 6-base consensus sequence of aggagg, the sd sequence recruits and aligns ribosome to mrna in order to initiate translation  <cit> .
motif pairs related to sd sequences were found in seven of the  <dig> investigated genomes.
the left part of the motif pair corresponds to the sd sequence whereas the right part contains the start codon.
the sd-related motif pairs arise from the sd sequence being located at the same distance from the start codon in many genes .
the two campylobacter species and helicobacter hepaticus have exactly the same patterns, as do the two thermosipho species.
exiguobacterium sibiricum only differs from the thermosipho at one position ahead of start codon ; actually all  <dig> patterns are similar.
interestingly, the start codon atg in all these patterns is followed by a or aa, suggesting that this start codon context is preferred in these genomes.
we speculate that the primary reasons why sd-related motif pairs are found in only seven of the  <dig> investigated genomes is that the percentage of genes with recognizable sd sequences varies widely among different species  <cit>  and the start codon context is also variable, which leads to significant overrepresentation of exact pentamer pairs in some but not all genomes.table  <dig> genomes with shine-dalgarno sequence identified
start codon is shown in bold
protein-related motifs
the last major class of motif pairs detected are protein-related motifs, which reside mostly in protein-coding genes.
the list of such motifs found in the  <dig> genomes that were investigated in detail is provided in table  <dig>  two types of proteins stand out: one is abc transporter atp-binding protein and the other one is pe/ppe family protein.table  <dig> genomes with significant protein-related motif pairs
tacga 59acgac
ggtgg64
ggtgg61
gatga ggtgg64
ggtgg61
ggtgg61
tatca70
aactc25
cacct52
complete codons are highlighted in bold face
abc transporters are transmembrane proteins that utilize atp hydrolysis for translocation of various substrates across membranes and are widespread in all phyla from prokaryotes to humans  <cit> ..motif pairs related to abc transporter atp-binding subunits were found in  <dig> different species.
the motif pair ggtgg64gaacc is significant in leuconostoc mesenteroides and staphylococcus aureus whereas two streptococcus species have a significant motif pair ggtgg61gatga .
when combined, the motif pairs reveal an extended consensus ggtgg60gagaacc, which translates to an amino acid motif gg-x-dep, where g, d, e, and p denote glycine, aspartic acid, glutamic acid and proline, respectively, and x denotes any amino acid.
the gg is part of the characteristic abc transporter motif lsggq and the dep is part of the walker b box, which is highly conserved in the abc transporter protein family  <cit> .
the pe and ppe proteins were first reported in the genome sequence of mycobacterium tuberculosis strain h37rv and were subsequently identified in all mycobacterial species as well as some rhodococcus and nocardia farcinica genomes  <cit> .
the proteins are characterized by presence of a pe or ppe domain, respectively, which is a  <dig> amino acid residue conserved region located near the c-terminus.
the pe/ppe domain comprises a pentapeptide sequence motif gxsxg/s at the n-terminus and conserved amino acid residues ser, asp and his  <cit> .
however, the conservation of the patterns found in mycobacterium is due to conserved amino acids in the protein sequence other than this signature domain.
the pattern in a. citrulli arises from yd repeats in the protein.
motifs found in cytochrome c translate to cxxch, which is highly conserved in this protein.
potential novel motifs and other motifs of interest
while many motif pairs detected by our algorithm can be traced to known motifs with known biological functions, the method also shows potential in detecting novel motifs, such as the motif pair located downstream of trna genes in methanopyrus kandleri that might serve as target signal for trna gene 3′ end maturization enzyme.
other potential novel motifs include the motif pair cgaaa19tttcg in mycobacterium species and some other motifs located in protein coding genes that are not due to known conserved amino acids motifs .
we selected a short list of promising potential novel motifs from additional file 2: table s <dig> based on the spacer distance and the percentage of copies in intergenic regions .
the rationale is that the motifs with regulatory functions are likely to reside in intergenic regions whereas motif pairs in protein coding regions are more likely to be related to properties of the encoded proteins, which, while not without interest, are not the main target of this work.
the other criterion is that we chose motif pairs with spacer lengths different from known repeats, such as the cripspr-related motifs.
properties of each of these motif pairs based on the analysis with pattern locator  are briefly summarized in table  <dig> table  <dig> selected potential novel motifs
a significant motif pairs comprised of overlapping pentamers were combined and reported as a single motif pair
the motif pair ttaat5attaa of a. radiobacter has  <dig> copies in the genome, all of which are located in intergenic regions and mostly about 100 bp upstream of 5′ ends of genes .
that indicates their potential functions in gene regulation.
queries to prodoric  <cit> , a database of prokaryotic regulatory interactions, did not produce a match to a motif of the same or similar consensus sequence.
the motif pair cctac19taggt in d. alkenivorans has  <dig> copies with  <dig> of them located in intergenic regions.
the distances to nearest gene vary substantially from less than 50 bp to over 500 bp.
eighteen copies of the motif are located between convergent genes , which argues against a possible role in regulation of transcription initiation.
there is no obvious relationship to genes related to any particular function.
the gcttat5aagcg motif pair in a. ebreus has  <dig> copies in the genome of which  <dig> are in genes,  <dig> are in intergenic regions and  <dig> partially overlap with genes.
orientation of adjacent genes for intergenic patterns are  <dig> convergent ,  <dig> divergent  and  <dig> co-oriented.
whether the patterns are located upstream or downstream, their distances to the 5′ or 3′ ends of genes are generally small .
in an analogy to the motif pair in d. alkenivorans, frequent occurrence between convergently transcribed genes argues against a possible role in transcription initiation.
the motif gaatccat23atggattc in b. japonicum has  <dig> copies with  <dig> of them in intergenic regions and  <dig> overlapping with 5′ ends of genes.
almost all motif pairs are proximal  to the 3′ end of a gene, indicating that they might be involved in transcription termination, although the gap between the two parts of this palindromic motif  is larger than the typical loop size in the hairpin structures formed by rho-independent transcription terminators.
majority of the adjacent genes encode hypothetical proteins of unknown function.
the motif pair atagct22caaaag in p. naphthalenivorans has  <dig> copies in the genome of which  <dig> are in intergenic regions,  <dig> in genes and  <dig> overlap with a gene’s start or end.
<dig> copies are between convergently transcribed genes, again suggesting that the main role of this motif pair is not related to transcription initiation.
similar to the gcttat5aagcg motif pair in a. ebreus, the  <dig> copies in intergenic regions are often close  to adjacent genes.
the motif pair gggacag15tgtccc of b. animalis has  <dig> copies in the genome, one of which overlaps with a 5′ end of a gene and the other  <dig> are in intergenic regions.
ten copies of the motif pair feature at least one adjacent gene involved in amino acid metabolism or amino acid modifications suggesting a potential role in amino acid metabolism.
however, most copies of the motif are not proximal to 5′ ends of genes  probably arguing against a role in transcription initiation.
the attata18gtcaa motif pair of p. mobilis 18tataat) has  <dig> copies in the genome of which  <dig> are in intergenic regions, three are in genes, one partially overlaps with the 5′ end of a gene, and four overlap with the 3′ end of a gene.
notably,  <dig> copies are within 50 bp upstream of a trna gene.
four additional copies are next to ribosomal protein genes and  <dig> are next to genes annotated as transcription regulators.
the motif conforms to the consensus − <dig> and − <dig> sites of bacterial promoters and the precise adherence of these specific genes to the promoter consensus probably reflects the importance of these genes and their efficient transcription.
the motif pair ttgac19ataat in c. kluyveri also matches the consensus − <dig> and − <dig> site sequences.
it has  <dig> copies of which  <dig> are in intergenic regions,  <dig> are in genes  <dig> is overlap with gene.
over 1/ <dig> of the copies found in intergenic regions are within 50 bp from the 5′ end of genes.
discussion
ab initio detection of candidate functional sequence motifs in complete genomes traditionally relies on word-counting approaches, which follow from a reasoning that selective constraints on a functional sequence motif could lead to statistically significant excess  of the motif occurrences in the genome.
consequently, such methods aim to detect specific words  in dna sequence that occur significantly more or less often than expected in some null model, generally based on a markov chain representation of the sequence .
other commonly used unsupervised motif-finding methods utilize a probabilistic model of the motif, which allows formulating the motif-finding problem as an optimization task.
the optimization is typically solved by one of the markov chain monte carlo  class of methods .
however, the latter methods are generally applied to a collection of relatively short sequences  and are not well suited for application to complete genomes.
our approach is similar to word-counting methods but differs from methods described above in aiming to detect pairs of motifs that exhibit a significant preference for occurrence at a specific distance from each other rather than unexpectedly high occurrence relative to some null model of a random sequence.
another important feature of our algorithm is exclusion of motif pairs that feature identical or similar spacers.
this step is designed to filter out motif pairs that arise from duplications of sequence fragments containing the motif pair as well as the spacer sequence, such as those resulting from proliferation of transposons, duplicated genes, or other forms of  sequence repeats.
consequently, our method specifically targets motif pairs that are likely maintained by selective constraints acting on the specific relative distance of the motifs, unlike motifs that arise from sequence duplication and may not be functional per se.
the current implementation of the method that was tested in this work is limited to pairs of pentanucleotides and spacer lengths up to  <dig> nucleotides.
these limitations are partially driven by considerations of computational complexity of applying the software for analysis of hundreds of complete genomes.
however, the pentanucleotides are also a reasonable choice with respect to dna-protein interactions, where the sequence-specific bonds form between the exposed amino acid side chains and the base pairs in the dna major groove, often leading to patterns of four to five adjacent conserved positions  <cit> .
moreover, pairs of motifs longer than five nucleotides can still be detected as pairs of overlapping pentanucleotides .
investigating locations of the detected motif pairs relative to genes and other known functional elements can often provide hints regarding their possible physiological roles or at the minimum help exclude some possible functions.
for that purpose, our motif detection pipeline is combined with pattern locator, which provides basic statistics on locations of sequence patterns in the genome  <cit> .
detailed investigation of motif pairs found in  <dig> selected genomes allowed us to link some of the motif pairs to known biological functions, such as rho-independent transcription terminators, crispr elements, and combinations of the transcription factor binding site and translation initiation site , which occur at specific distance from each other .
a number of conserved motif pairs were found in genes and we speculate that many of such patterns can arise from constraints on the properties of the encoded proteins, although only some of these patterns are related to previously known protein sequence motifs .
a novel motif pair detected in methanopyrus kandleri is of particular interest due to its location downstream of trna genes .
this motif is unlikely to be detected by motif-finding software designed for detection of transcription factor binding sites because this particular motif is located outside the regulatory regions where transcription factor sites generally reside.
the word-counting methods that are more suitable for scanning complete genomes for overrepresented sequence motifs are also unlikely to detect such a motif due to their low sensitivity to spaced motif pairs separated by a variable spacer.
due to the unique mechanism of trna maturation in m. kandleri  <cit> , we speculate that this motif could be involved in the trna maturation process, possibly as a binding site for the maturation enzyme.
it is important to note that the motif pairs listed in tables  <dig>   <dig>   <dig>   <dig> and  <dig> represent only a fraction of the significant motif pairs detected by our methodology, namely those that we were able to link to a likely physiological function.
the full set of motif pairs detected in  <dig> genomes  contains novel motifs that may not be detectable by other methods.
due to the conservative criteria for the detection of motif pairs and filtering of motif pairs with conserved spacers, we proffer that many of the motif pairs listed in additional file 2: table s <dig> arise from selective constraints on the relative positions of the motifs.
such motif pairs could play important physiological roles and at the same time escape detection by standard motif-finding techniques.
a subset of selected candidate novel motifs  was subjected to additional analysis that allows narrowing down a possible range of functions of these motifs.
we conclude that method for detection of significant motif pairs described above is capable of detection of novel functional elements that are unlikely to be discovered by other method and we offer the paired-motif methodology as a tool to complement mainstream motif-finding methods for detection of functional elements in complete genomes.
CONCLUSIONS
we present a novel motif-finding method based on detection of pairs of sequence motifs with statistically significant preference for a specific distance from each other.
the method is suitable for application to long dna sequences including complete genomes.
application to  <dig> prokaryotic genomes yielded a list of  <dig> significant motif pairs  including many that are related to sequence motifs of known functions as well as novel sequence motifs.
the complete list of all significant motif pairs found in the  <dig> genomes is available in additional file 2: table s <dig>  we propose the methodology as a suitable complement of the set already existing techniques for discovery of functionally significant sequence motifs.
additional files:
additional file  <dig>
table s <dig>  comparison of significant pentamer pairs detected with different sequence similarity cutoffs used to eliminate duplicated spacers.
description of data: comparison of significant pentamer pairs detected with sequence similarity cutoffs 65%, 70%, and 75% in three genomes.
additional file 2:
table s <dig> .
significant motif pairs detected in  <dig> analyzed genomes.
description of data: the table contains  <dig> motif pairs in  <dig> genomes that identified as significant by our approach.
additional file 3:
table s <dig>  summary statistics for each analyzed genome.
description of data: the table contains mean statistics for variable of interests and distribution of spacer length at the genome level in which at least one significant motif pair is detected.
additional file 4:
table s <dig>  list of the  <dig> genomes selected for detailed analysis.
description of data: the table contains the  <dig> genomes that are selected to investigate in detail in this paper.
abbreviations
bpbase pairs
sdshine-dalgarno
tfbstranscription factor binding site

 genomic dna sequences contain numerous repetitive motifs or nucleotide sequence patterns that play various roles in the physiology of the cell, including regulation of gene expression, maintenance of dna and its structure in the cell, recombination and other forms of mutation, and even recognition of dna fragments for uptake from the extracellular space. repetitive motifs in dna sequences can also encode functional elements in proteins and rna. consequently, a large amount of effort has been devoted to development of computational methods for detection of sequence motifs in dna and protein sequences. there are two major types of motif finding algorithms, namely supervised and unsupervised motif finding algorithms. the former methods require a sample of known occurrences of the motif and utilize this information in the search for additional motif occurrences in the analyzed sequence or sequences. the unsupervised methods, sometimes referred to as ab initio approach, do not require any prior knowledge about the motif sequences and detect novel sequence motifs that satisfy specified criteria. in this article, we focus on unsupervised motif search. the unsupervised motif finding algorithms can be further classified into two major groups: 1) word-based methods that mostly rely on exhaustive enumeration, i.e., counting and comparing oligonucleotide frequencies and 2) probabilistic sequence models where the model parameters are estimated from sequences. extensive work has been done on transcription factor binding site prediction during the past decades, driven by the obvious importance of these regulatory motifs in the organism’s physiology. based on the type of dna sequence information used by the tfbs finding algorithm, the methods could be classified into three major classes: 1) methods that use promoter sequences from coregulated genes from a single genome, 2) methods that use orthologous promoter sequences of a single gene from multiple species and 3) methods combining 1) and 2). as a unified portal for online discovery and analysis of sequence motifs, the meme suite web server provides various tools in finding motifs representing features such as dna binding sites and protein interaction domains. while tfbs receive most attention among sequence motifs the repetitive sequence motifs in genomic dna can have many other functions. we aim to expand the range of types of sequence motifs detected by motif-finding methods by searching for spaced motifs in complete genomes, which could arise, among other mechanism, from repetitive patterns in chromosome structure. the concept of searching for spaced sequence motifs is not new but the previous such methods generally aimed to detect tfbs, where spacing of the conserved segments of the motif is determined by the geometry of the dna-protein interaction and generally does not exceed or 7 bp. a more general approach implemented in helicis allows detection of co-localized periodically spaced motifs, such as binding sites for multiple transcription factors. however, these methods are specifically designed for detection of tfbs in a collection of regulatory regions and are not suitable for scanning complete genomes. in contrast, our methodology is aimed at detection of dna sequence motifs that can have diverse physiological functions, including regulatory motifs, motifs involved in dna, rna, and protein interactions, and motifs related to structural organization of dna and/or the encoded proteins. our approach is different from previous methods in that it searches for pairs of motifs that occur at a particular distance from each other significantly more often than expected without requiring that either of the motifs on its own occurs more frequently than expected. the goal is to detect paired motifs where the distance between the two motifs is under selective constraint but the specific sequence of the spacer is variable. methods the motif-finding algorithm the goal of the algorithm is to find motif pairs that have a preferred distance of separation. that is, they occur at a specific distance apart more often than expected by chance. the algorithm consists of the following steps:scan the genome and count the number of times that the pairs of motifs a and b occur at every spacer distance between and. record the total number as n; identify the distance at which the pair has the maximum number of occurrences, record the maximum number as fmax; test whether the maximum number of occurrences fmax is large enough to claim significance of overrepresentation by surpassing some appropriate cutoff value fcut; repeat steps 1) – 3) for all possible pentamer pairs; record all significant pairs with fmax > fcut; for each significant pair, align the spacers pair-wise. reduce fmax by for each pair of aligned spacers with nucleotide identity greater than 70% and remove one of the two spacers from further comparisons with other spacers; identify and record motif pairs whose fmax is still above fcut. record if the inverse complement is also significant with the same distance; for motif pairs that satisfy the criteria above, investigate their distribution in the genome using pattern locator. step 1) consists of recording the total number of times a pair of pentamers a and b is found at the exact mutual distance d. the distance d is counted from the beginning of a to the beginning of b. pairs with exactly the same spacers are counted only once. the motivation for discounting motif pairs with identical or very similar spacers is that they likely arise from sequence duplication of the whole segment containing the motif pair and the spacer, which is not our primary goal; we aim to detect motif pairs where the distance between the motifs arises from selective constraints rather than as a result of sequence duplication events. the main reason we limit the motif spacer distance to 5–89 bp is due to growing computational complexity when ranges of paramters investigated are too large. similar reasons apply to limiting the investigated motifs to pentamers; the number of possible pairs of hexamers is around million, which is prohibitive when the methodology is to be applied to hundreds of complete genomes as in this study. another important motivation to use pentamers is related to how prokaryotic dna-binding proteins, particularly transcription factors, interact with dna. many prokaryotic transcription factors interact with dna in a sequence-specific manner via two alpha-helices that fit into adjacent turns of the major groove in the dna double helix. consequently, although the whole binding site is often ~ nucleotides long, it usually consists of two conserved segments of ~ nucleotides where the protein directly interacts with the dna bases, separated by a variable gap. after finding the most abundant distance d for each motif pair in step 2), statistical significance is assessed in step 3). for a given pair of motifs, if no specific distance is preferred, one would expect a uniform distribution of spacer distances within the whole range to 89 bp. the problem of numbers of motif pairs found at different distances is analogous to the classical balls-in-urns problem where motif pairs represent the balls and different distances represent the urns. suppose n balls are allocated to k urns at random, each ball being equally likely to fall into each urn, independently of any other balls. let fmax denote the largest number of balls in any urn. this situation can be modeled with the uniform multinomial distribution in which there are n trials, each with k possible outcomes of equal probability. david and barton showed that p can be computed as p=n!kn×thecoefficientofλnink subsequently, the p-value p = p =  - p. the analogy in dna sequences is straightforward with n being the total number of motif pair occurrences with spacer distance between and 89 bp, k being the number of different distances and fmax being the maximum number of occurrences of a specific distance within the range. a bonferroni correction was applied to the p-value cutoff in order to account for multiple hypothesis testing: p-value-cutoff =  / number of motif pairs =  / 410 ≈ 10− testing over multiple distances is accounted for in the uniform multinomial. sympy, a python library for symbolic mathematics, was utilized to perform the power expansion in formula and to extract the coefficient needed to compute p. however, using sympy to obtain the exact tail probability is computationally expensive, especially for large parameters of n and k. hence, we used an approximation to calculate the 100% percentile of fmax with the following equation fα≈nk+nlogkk−loglog11−α+ n2klogk in the formula above, fα is the 100% percentile of fmax; α is the desired type i error rate; n is the total number of occurrences of the motif pair at hand with distance between and and k is the total number of distances investigated. note that the above approximation is based on an asymptotic theorem under certain regularities. to verify that the approximation is sufficiently accurate with our set of parameters, comparisons between approximate and exact fcut values for different values of k were made using a fixed p-value 10−. when k is relatively large, the equation approximates the true values with reasonable accuracy as n gets large; however, when k is relatively small, the equation tends to overestimate the true fcut values. since the k used in this paper is we decided to use exact cutoff values for n < = while the approximated values are used for n >  fig. comparison between approximate and exact p-value cutoffs.n is the total number of occurrences of some motif pair; fcut indicates the cutoff value corresponds to the selected p value threshold. solid and dashed lines represent approximate values calculated by formula and exact values computed with sympy, respectively in step 4], the steps to are repeated for all possible pentamer pairs. after all pentamer pairs are tested, the most promising candidates are selected for spacer alignment. in step 1], only identical spacers could be recognized and eliminated but our goal is to discount all motif pairs which likely arose from recent duplication of the whole segment including the motif pair and the spacer. consequently, sequence alignments are performed with nw-align, which is a simple and robust alignment program based on the standard needleman-wunsch algorithm, to eliminate potential duplicated spacers. a simple ‘greedy’ algorithm is applied to eliminate duplicated spacers. the first pair of spacers is aligned. if the nucleotide identity is ≥70%, then fmax is reduced by one and one of the spacers is excluded from further alignments. additional pairs of spacers are subsequently aligned and excluded if the nucleotide identity exceeds 70%. this is continued until no spacers are more similar than the 70% cutoff. the value 70% was chosen as the lowest cutoff that gives a high confidence the sequence similarity is not coincidental; for example, the probability that two random sequences of ten nucleotides would match at seven positions is less than and declines with increasing sequence length. we also performed tests with 65% and 75% cutoffs, which suggested that the cutoff 65% was too conservative whereas majority pentamer pairs identified with the 75% cutoff overlapped motifs detected also with the 70% cutoff. after this step, a motif pair whose fmax is still above fcut and whose inverse complement is also significant with the same distance is selected for further analysis. the reason for requiring a motif pair’s inverse complement also to be significant is that dispersed sequence motifs, which are the primary target of our search, are expected to be distributed in both dna strands; consequently, both versions of the motif pair—direct and reverse—are expected to be over-represented in the genome. the subsequent analysis of the significant motif pairs is performed with pattern locator, which provides information on distribution of the motif pairs with respect to adjacent genes. such information can be helpful in generating hypotheses about the motif pair’s possible function. pattern locator also reports the percentage of copies of the motif pair found in genes and intergenic regions, which allows easy identification of motifs that are most likely to have regulatory functions and motifs that could reflect conserved patterns in amino acid sequences of the encoded proteins. dna sequences annotated nucleotide sequences of complete prokaryotic genomes were downloaded from the ncbi ftp server. we randomly selected only one genome per species when multiple strains of the same species were available. the final dataset included genomes. for genomes consisting of multiple chromosomes the analysis was performed on all chromosomes.  the prokaryotic genomes, each from a separate species, were scanned with the algorithm described above. in the genomes, motif pairs were identified as significant by our criteria. the summary statistics for selected characteristics of the significant motif pairs are shown in table  table summary statistics of motif pairs identified in genomes the summary statistics is per motif. meaning of abbreviations in the table: d, spacer length of a motif pair; initial_n, number of copies of the motif pair before alignment; reduced_n, number of copies after alignment and elimination of duplicate spacers; difference, the difference between initial_n and reduced_n; cutoff, significance cutoff; reduced_n/cutoff, the ratio of reduced_n and cutoff, indicating the relative significance for the motif pair; gene, the percentage of each motif pair occurrences found in genes; intergenic, the percentage of each motif pair occurrences that are in intergenic regions; overlap, the percentage of each motif pair occurrences that overlap with a gene start or end. these percentages are calculated as follows: for any given significant motif, we run a query with pattern locator, which gives the percentage of the motif occurrences that fall in gene, intergenic region or overlap with gene starts or ends. the quantiles in the table are for these percentages over all significant motifs distributions of some key statistics are further demonstrated in figs.  and the right panel in fig.  indicates that median percentage of significant pairs found in protein- coding regions is around 80%, which is close to the fraction of protein-coding region in a typical prokaryotic genome; in other words, the significant motif pairs do not show a strong overall bias for protein-coding regions or intergenic regions. in fig.  the first two plots show that the distributions of motif pairs are bi-modal in both gene and intergenic region. percent of significant motif pairs occurred in genes more than 90% of the time and percent occurred in genes less than 10% of the time, whereas fewer motif pairs are distributed approximately evenly among both genes and intergenic regions. in other words, majority of the significant motifs show a strong preference for either genes or intergenic regions. the same conclusion can be reached upon inspection of the raw data in additional file 3: table s fig. distribution of key motif statistics. frequency is number of copies of each significant motif pair; percentage is the percentage of significant motif pairs found in corresponding region fig. distribution of motif pairs with respect to gene and inter-genic region. in each plot, the horizontal axis indicates the percentage of occurrences of significant pairs found in corresponding region the distribution of the spacer length for the motif pairs is shown in fig.  there are two peaks—one is around spacer lengths to 20 bp and the other one is around to 50 bp. the reason for the first peak is partly due to a conserved palindromic motif pair, which is widely distributed among a variety of genes of mycobacterium species. in total, approximately motif pairs conforming to this consensus were identified in the analyzed mycobacterium genomes out of motif pairs with spacer length to 20 bp found in all analyzed genomes. the peak at 40–50 bp reflects wide distribution of clustered regularly interspaced short palindromic repeats among the analyzed genomes. crisprs are arrays of ~30 bp perfect repeats separated by ~30–40 bp variable spacers, which give rise to many significantly overrepresented motif pairs with this spacer size.fig. spacer length distribution of candidate pairs another noteworthy statistic in table  is the difference between the initial number of copies of a significant motif pair and the number of copies remaining after elimination of duplicate spacers. the mean value of initial_n is whereas reduced_n has mean, indicating that more than 25% of occurrences of all significant motif pairs in the genomes may have arisen from sequence duplication events. this confirms the importance of removing the “false positives” resulting from sequence duplications when searching for candidate motif pairs whose specific relative positioning is maintain by selective constraints. all significant motif pairs found in the analyzed genomes are reported in additional file 3: table s in total, genomes have at least one significant motif pair identified. among those, rhodoferax ferrireducens t has the highest number of pairs detected,. thirty genomes were selected for in-depth investigation of the identified significant motif pairs. some of the significant motif pairs were found to be related to previously known motifs such as crispr or shine-dalgarno sequence whereas others appear to be potentially novel, previously uncharacterized sequence motifs. majority of the identified motif pairs can be classified in one of the five categories, namely crispr-related, rho-independent transcription terminators, trna-related motifs, motifs containing the shine-dalgarno sequence, and protein-related motifs. crispr-related motif pairs crisprs are segments of prokaryotic dna containing periodic repetitions of ~ nucleotides. each repetition is followed by short segments of “spacer dna”, typically 30– nucleotides long. the spacers are short fragments of nucleotide sequence from phages or plasmids to which the bacterium had been exposed in the past. these fragments in conjunction with the cas proteins are subsequently used to recognize and destroy the dna of the same phages or plasmids in future encounters. the crispr/cas system is a form of prokaryotic immune system that confers resistance to foreign genetic elements such as plasmids and phages by recognizing and digesting them. crisprs are found in approximately 40% of sequenced bacterial genomes and 90% of sequenced archaea. given its properties of regularly interspaced repeats with spacers of variable sequence but constant length, it’s not surprising that parts of crisprs are among the paired motifs detected by the algorithm. table  lists significant motif pairs found in the genomes subjected to in-depth analysis that are related to crispr. results were verified by comparison with crisprdb, a comprehensive database of known crispr occurrences. candidatus desulforudis audaxviator mp104c contains two types of crispr sequences, gtttcaatccctcgtaggtaggctggaaac and ctttcagtccccttttcgtcgggtcggtcgctgaaac. the first four patterns listed in table  for this genome correspond to the former whereas the last three reflect the latter crispr sequence. the patterns consist of two pentamers near the ends of the crispr sequence and the spacers, which are variable in sequence bur have a constant length, are responsible for the fixed distance between the pentamers. for patterns detected in clostridium thermocellum atcc the first two arise from repeats gtttcaattcctcataggtacgataaaaac and the rest are due to gttttatcgtacctatgaggaattgaaac. all patterns found in corynebacterium aurimucosum atcc are due to the crispr sequence gtgctccccgcgtaagcggggatgagcc. similarly, all patterns detected in the other two corynebacterium species are due to the same consensus ggctcatccccgctggcgcggggagcac.table significant motif pairs related to crispr repeats gaaac35ggttt aaact35gtttc gaaac36gtttc tgaaa39tttca ctgaa40tttca gtacg52tcctc atgag53tcgta tgagg53cgtac atgag54cgtac accta55tatcg cctat56tcgta accta56atcgt cgggg44tcccc cgggg45ccccg gggga45cccgc cgggg44tcccc cgggg45ccccg gggga45cccgc cgggg44tcccc gggga45cccgc pattern x1x2x3x4x5dy1y2y3y4y denotes 5-mer x1x2x3x4x is followed by 5-mer y1y2y3y4y with d nucleotides apart; for non-palindromic patterns, only sequence in one dna strand is listed note that our software is not intended for detection of crispr sequences and the detection of crispr-related sequence pairs are a byproduct of the characteristic crispr structure. our method is designed for detection of dispersed sequence motifs, whereas crisprs consist of repeats clustered at a small number of loci in the genome. this explains why crispr-related motif pairs were not detected in many genomes with crispr. the main reason is the requirement that both direct and reverse complement versions of a motif pair independently satisfy the criteria for statistical significance. many genomes contain a single crispr locus, which means that all crispr sequences are in the same orientation with respect to the direct and reverse strands. such cripsr loci can only be detected if they contain a palindromic pair of pentamers that satisfies all other criteria. rho-independent transcription terminators a transcription terminator is a functional sequence located at the 3′ end of genes that mediates the transcriptional termination. in prokaryotes, there are two types of terminators, namely rho-dependent and rho-independent terminators. the former require rho factor to terminate the transcription process while the latter one forms a self-annealing hairpin structure to serve the purpose. the terminator sequence generally contains a gc-rich region of dyad symmetry followed by a short poly-t tract. motif pairs related to transcription terminators were found in out of the genomes. the results were verified by transtermhp, which is a specialized software to find rho-independent transcription terminators in bacterial genomes. in actinobacillus pleuropneumoniae the pattern is a high-g + c palindrome; in haemophilus influenzae the pattern detected by our method is not a palindrome, but instead includes a poly-t tract. both features in the two patterns identified correspond to the standard structure of rho-independent terminators, which consist of a gc-rich palindrome followed by a short poly-t tract.table motif pairs that are part of transcription terminators caaaa15accgc ggtct15ttgca tgcaa15cgacc tgcaa16aaccg gcaaa16accgc ttgca16agacc gcggt16tttac gcaaa17ccgct cggtc17tgcaa tttgc18gaccg aagcg19ttgca gcggt19gcaaa ttttg19gaccg cggtc38gaccg ggtca11gtttt cggtc12gtttt gcggt12agttt gcggt12tgttt ggtca12ttttt cggtt12gtttt tgcgg13agttt cggtc13ttttt gcggt13gtttt tgcgg13tgttt cggta13ttttt cggtt13ttttt gcggt13atttt gcggt14ttttt gtgcg14tgttt tgcgg14gtttt tgcgg15ttttt agtgc15tgttt gtgcg15gtttt gcggt15ttttg agtgc16gtttt aaaca16cactt tgcgg16ttttg the fact that transcription terminator-related motif pairs were detected by our program is actually unexpected because the terminators are characterized by the stem-loop secondary structure but not a specific repeated sequence. this is also the reason why we did not detect transcription terminator-related motif pairs in other genomes. however, the h. influenzae genome contains about copies of the uptake signal sequence with the consensus aagtgcggt and terminators in h. influenzae are often composed of two copies of the uss in the opposite orientations that form the stem of the stem-loop structure. that leads to prevalence of this specific sequence in the terminator stem-loop structures and it is the reason why we are detecting terminators in h. influenzae. analysis of the a. pleuropneumoniae genome with the aimie software revealed an abundant motif aagcggt, which is similar to the uss in h. influenzae and may have an analogous function. like the uss in h. influenzae, the a. pleuropneumoniae motif frequently occurs in dyad pairs found at the 3′ ends of genes, which probably serve as transcription terminators. trna-related motifs trna-related motifs were detected in four genomes – of them are inside trna genes and correspond to parts of a conserved pattern tagagc27ggttcg near the 5′ end of the gene. parts of this pattern are highly conserved in trna genes from diverse species. however, the motif pairs in a. butzleri, s. muelleri, and n. profundicola just exceed the significance cutoff of our method and they are not detected in other genomes because the number of copies is not statistically significant.table trna-related motifs identified the pattern in methanopyrus kandleri av differs from the other three in that it resides outside of trna genes, specifically about base pairs downstream of the trna genes. other trna genes have similar patterns at the same position. the consensus sequence of the motif is shown in fig.  fig. sequence logo for motif pair downstream of trna genes in m. kandleri the sequence is palindromic with positions 2– and 15– being the most conserved. unlike the other trna-related motif pairs, this motif appears to be specific to m. kandleri and we did not find similar motifs in other genomes. the methanogenic archaeon m. kandleri grows at temperatures of 84–110 °c and genome analyses revealed strategies to adapt to these harsh conditions. among other adaptations, m. kandleri uses a unique trna c-to-u editing mechanism at base for different trna species. we speculate that the palindrome discovered by our software might serve as a target signal for the 3′ maturation enzyme and its presence specifically in m. kandleri may be related to the unique mechanism of trna maturation employed by this organism. motifs related to shine-dalgarno sequence located around 8 bp upstream of the start codon, shine-dalgarno sequence serves as ribosome binding site in mrna of prokaryotic genomes. with a 6-base consensus sequence of aggagg, the sd sequence recruits and aligns ribosome to mrna in order to initiate translation. motif pairs related to sd sequences were found in seven of the investigated genomes. the left part of the motif pair corresponds to the sd sequence whereas the right part contains the start codon. the sd-related motif pairs arise from the sd sequence being located at the same distance from the start codon in many genes. the two campylobacter species and helicobacter hepaticus have exactly the same patterns, as do the two thermosipho species. exiguobacterium sibiricum only differs from the thermosipho at one position ahead of start codon; actually all patterns are similar. interestingly, the start codon atg in all these patterns is followed by a or aa, suggesting that this start codon context is preferred in these genomes. we speculate that the primary reasons why sd-related motif pairs are found in only seven of the investigated genomes is that the percentage of genes with recognizable sd sequences varies widely among different species and the start codon context is also variable, which leads to significant overrepresentation of exact pentamer pairs in some but not all genomes.table genomes with shine-dalgarno sequence identified start codon is shown in bold protein-related motifs the last major class of motif pairs detected are protein-related motifs, which reside mostly in protein-coding genes. the list of such motifs found in the genomes that were investigated in detail is provided in table  two types of proteins stand out: one is abc transporter atp-binding protein and the other one is pe/ppe family protein.table genomes with significant protein-related motif pairs tacga 59acgac ggtgg64 ggtgg61 gatga ggtgg64 ggtgg61 ggtgg61 tatca70 aactc25 cacct52 complete codons are highlighted in bold face abc transporters are transmembrane proteins that utilize atp hydrolysis for translocation of various substrates across membranes and are widespread in all phyla from prokaryotes to humans..motif pairs related to abc transporter atp-binding subunits were found in different species. the motif pair ggtgg64gaacc is significant in leuconostoc mesenteroides and staphylococcus aureus whereas two streptococcus species have a significant motif pair ggtgg61gatga. when combined, the motif pairs reveal an extended consensus ggtgg60gagaacc, which translates to an amino acid motif gg-x-dep, where g, d, e, and p denote glycine, aspartic acid, glutamic acid and proline, respectively, and x denotes any amino acid. the gg is part of the characteristic abc transporter motif lsggq and the dep is part of the walker b box, which is highly conserved in the abc transporter protein family. the pe and ppe proteins were first reported in the genome sequence of mycobacterium tuberculosis strain h37rv and were subsequently identified in all mycobacterial species as well as some rhodococcus and nocardia farcinica genomes. the proteins are characterized by presence of a pe or ppe domain, respectively, which is a amino acid residue conserved region located near the c-terminus. the pe/ppe domain comprises a pentapeptide sequence motif gxsxg/s at the n-terminus and conserved amino acid residues ser, asp and his. however, the conservation of the patterns found in mycobacterium is due to conserved amino acids in the protein sequence other than this signature domain. the pattern in a. citrulli arises from yd repeats in the protein. motifs found in cytochrome c translate to cxxch, which is highly conserved in this protein. potential novel motifs and other motifs of interest while many motif pairs detected by our algorithm can be traced to known motifs with known biological functions, the method also shows potential in detecting novel motifs, such as the motif pair located downstream of trna genes in methanopyrus kandleri that might serve as target signal for trna gene 3′ end maturization enzyme. other potential novel motifs include the motif pair cgaaa19tttcg in mycobacterium species and some other motifs located in protein coding genes that are not due to known conserved amino acids motifs. we selected a short list of promising potential novel motifs from additional file 2: table s based on the spacer distance and the percentage of copies in intergenic regions. the rationale is that the motifs with regulatory functions are likely to reside in intergenic regions whereas motif pairs in protein coding regions are more likely to be related to properties of the encoded proteins, which, while not without interest, are not the main target of this work. the other criterion is that we chose motif pairs with spacer lengths different from known repeats, such as the cripspr-related motifs. properties of each of these motif pairs based on the analysis with pattern locator are briefly summarized in table  table selected potential novel motifs a significant motif pairs comprised of overlapping pentamers were combined and reported as a single motif pair the motif pair ttaat5attaa of a. radiobacter has copies in the genome, all of which are located in intergenic regions and mostly about 100 bp upstream of 5′ ends of genes. that indicates their potential functions in gene regulation. queries to prodoric, a database of prokaryotic regulatory interactions, did not produce a match to a motif of the same or similar consensus sequence. the motif pair cctac19taggt in d. alkenivorans has copies with of them located in intergenic regions. the distances to nearest gene vary substantially from less than 50 bp to over 500 bp. eighteen copies of the motif are located between convergent genes, which argues against a possible role in regulation of transcription initiation. there is no obvious relationship to genes related to any particular function. the gcttat5aagcg motif pair in a. ebreus has copies in the genome of which are in genes, are in intergenic regions and partially overlap with genes. orientation of adjacent genes for intergenic patterns are convergent, divergent and co-oriented. whether the patterns are located upstream or downstream, their distances to the 5′ or 3′ ends of genes are generally small. in an analogy to the motif pair in d. alkenivorans, frequent occurrence between convergently transcribed genes argues against a possible role in transcription initiation. the motif gaatccat23atggattc in b. japonicum has copies with of them in intergenic regions and overlapping with 5′ ends of genes. almost all motif pairs are proximal to the 3′ end of a gene, indicating that they might be involved in transcription termination, although the gap between the two parts of this palindromic motif is larger than the typical loop size in the hairpin structures formed by rho-independent transcription terminators. majority of the adjacent genes encode hypothetical proteins of unknown function. the motif pair atagct22caaaag in p. naphthalenivorans has copies in the genome of which are in intergenic regions, in genes and overlap with a gene’s start or end. copies are between convergently transcribed genes, again suggesting that the main role of this motif pair is not related to transcription initiation. similar to the gcttat5aagcg motif pair in a. ebreus, the copies in intergenic regions are often close to adjacent genes. the motif pair gggacag15tgtccc of b. animalis has copies in the genome, one of which overlaps with a 5′ end of a gene and the other are in intergenic regions. ten copies of the motif pair feature at least one adjacent gene involved in amino acid metabolism or amino acid modifications suggesting a potential role in amino acid metabolism. however, most copies of the motif are not proximal to 5′ ends of genes probably arguing against a role in transcription initiation. the attata18gtcaa motif pair of p. mobilis 18tataat) has copies in the genome of which are in intergenic regions, three are in genes, one partially overlaps with the 5′ end of a gene, and four overlap with the 3′ end of a gene. notably, copies are within 50 bp upstream of a trna gene. four additional copies are next to ribosomal protein genes and are next to genes annotated as transcription regulators. the motif conforms to the consensus − and − sites of bacterial promoters and the precise adherence of these specific genes to the promoter consensus probably reflects the importance of these genes and their efficient transcription. the motif pair ttgac19ataat in c. kluyveri also matches the consensus − and − site sequences. it has copies of which are in intergenic regions, are in genes is overlap with gene. over 1/ of the copies found in intergenic regions are within 50 bp from the 5′ end of genes. discussion ab initio detection of candidate functional sequence motifs in complete genomes traditionally relies on word-counting approaches, which follow from a reasoning that selective constraints on a functional sequence motif could lead to statistically significant excess of the motif occurrences in the genome. consequently, such methods aim to detect specific words in dna sequence that occur significantly more or less often than expected in some null model, generally based on a markov chain representation of the sequence. other commonly used unsupervised motif-finding methods utilize a probabilistic model of the motif, which allows formulating the motif-finding problem as an optimization task. the optimization is typically solved by one of the markov chain monte carlo class of methods. however, the latter methods are generally applied to a collection of relatively short sequences and are not well suited for application to complete genomes. our approach is similar to word-counting methods but differs from methods described above in aiming to detect pairs of motifs that exhibit a significant preference for occurrence at a specific distance from each other rather than unexpectedly high occurrence relative to some null model of a random sequence. another important feature of our algorithm is exclusion of motif pairs that feature identical or similar spacers. this step is designed to filter out motif pairs that arise from duplications of sequence fragments containing the motif pair as well as the spacer sequence, such as those resulting from proliferation of transposons, duplicated genes, or other forms of sequence repeats. consequently, our method specifically targets motif pairs that are likely maintained by selective constraints acting on the specific relative distance of the motifs, unlike motifs that arise from sequence duplication and may not be functional per se. the current implementation of the method that was tested in this work is limited to pairs of pentanucleotides and spacer lengths up to nucleotides. these limitations are partially driven by considerations of computational complexity of applying the software for analysis of hundreds of complete genomes. however, the pentanucleotides are also a reasonable choice with respect to dna-protein interactions, where the sequence-specific bonds form between the exposed amino acid side chains and the base pairs in the dna major groove, often leading to patterns of four to five adjacent conserved positions. moreover, pairs of motifs longer than five nucleotides can still be detected as pairs of overlapping pentanucleotides. investigating locations of the detected motif pairs relative to genes and other known functional elements can often provide hints regarding their possible physiological roles or at the minimum help exclude some possible functions. for that purpose, our motif detection pipeline is combined with pattern locator, which provides basic statistics on locations of sequence patterns in the genome. detailed investigation of motif pairs found in selected genomes allowed us to link some of the motif pairs to known biological functions, such as rho-independent transcription terminators, crispr elements, and combinations of the transcription factor binding site and translation initiation site, which occur at specific distance from each other. a number of conserved motif pairs were found in genes and we speculate that many of such patterns can arise from constraints on the properties of the encoded proteins, although only some of these patterns are related to previously known protein sequence motifs. a novel motif pair detected in methanopyrus kandleri is of particular interest due to its location downstream of trna genes. this motif is unlikely to be detected by motif-finding software designed for detection of transcription factor binding sites because this particular motif is located outside the regulatory regions where transcription factor sites generally reside. the word-counting methods that are more suitable for scanning complete genomes for overrepresented sequence motifs are also unlikely to detect such a motif due to their low sensitivity to spaced motif pairs separated by a variable spacer. due to the unique mechanism of trna maturation in m. kandleri, we speculate that this motif could be involved in the trna maturation process, possibly as a binding site for the maturation enzyme. it is important to note that the motif pairs listed in tables  and represent only a fraction of the significant motif pairs detected by our methodology, namely those that we were able to link to a likely physiological function. the full set of motif pairs detected in genomes contains novel motifs that may not be detectable by other methods. due to the conservative criteria for the detection of motif pairs and filtering of motif pairs with conserved spacers, we proffer that many of the motif pairs listed in additional file 2: table s arise from selective constraints on the relative positions of the motifs. such motif pairs could play important physiological roles and at the same time escape detection by standard motif-finding techniques. a subset of selected candidate novel motifs was subjected to additional analysis that allows narrowing down a possible range of functions of these motifs. we conclude that method for detection of significant motif pairs described above is capable of detection of novel functional elements that are unlikely to be discovered by other method and we offer the paired-motif methodology as a tool to complement mainstream motif-finding methods for detection of functional elements in complete genomes.  we present a novel motif-finding method based on detection of pairs of sequence motifs with statistically significant preference for a specific distance from each other. the method is suitable for application to long dna sequences including complete genomes. application to prokaryotic genomes yielded a list of significant motif pairs including many that are related to sequence motifs of known functions as well as novel sequence motifs. the complete list of all significant motif pairs found in the genomes is available in additional file 2: table s we propose the methodology as a suitable complement of the set already existing techniques for discovery of functionally significant sequence motifs. additional files: additional file table s comparison of significant pentamer pairs detected with different sequence similarity cutoffs used to eliminate duplicated spacers. description of data: comparison of significant pentamer pairs detected with sequence similarity cutoffs 65%, 70%, and 75% in three genomes. additional file 2: table s. significant motif pairs detected in analyzed genomes. description of data: the table contains motif pairs in genomes that identified as significant by our approach. additional file 3: table s summary statistics for each analyzed genome. description of data: the table contains mean statistics for variable of interests and distribution of spacer length at the genome level in which at least one significant motif pair is detected. additional file 4: table s list of the genomes selected for detailed analysis. description of data: the table contains the genomes that are selected to investigate in detail in this paper. abbreviations bpbase pairs sdshine-dalgarno tfbstranscription factor binding site 
unlike existing methods for motif discovery, the motifs themselves are not required to occur at unusually high frequency but only to exhibit a significant preference to occur at a specific distance from each other. in the present implementation of the method, motifs are represented by pentamers and all pairs of pentamers are evaluated for statistically significant preference for a specific distance. a number of methods have been developed for discovery of such sequence motifs with a primary focus on detection of regulatory motifs and particularly transcription factor binding sites. an important step of the algorithm eliminates motif pairs where the spacers separating the two motifs exhibit a high degree of sequence similarity; such motif pairs likely arise from duplications of the whole segment including the motifs and the spacer rather than due to selective constraints indicative of a functional importance of the motif pair. dna sequences contain repetitive motifs which have various functions in the physiology of the organism. the method was used to scan complete prokaryotic genomes for novel sequence motifs. we present a new motif-finding technique that is applicable to scanning complete genomes for sequence motifs. we conclude that our approach to detection of significant motif pairs can complement existing motif-finding techniques in discovery of novel functional sequence motifs in complete genomes. genomesequence motifsmotif-findingbacteriaarchaeadna sequence repeatshttp://dx.doi.org/ /100000153division of biological infrastructuredbi-0950266issue-copyright-statement© the author 2017 the results from analysis of genomes suggest that the method detects previously known motifs that are expected to be found as well as new motifs that are unlikely to be discovered by traditional motif-finding methods. selected motif pairs were subjected to further investigation and in some cases their possible biological functions were proposed.
 dna sequences contain repetitive motifs which have various functions in the physiology of the organism. a number of methods have been developed for discovery of such sequence motifs with a primary focus on detection of regulatory motifs and particularly transcription factor binding sites. most motif-finding methods apply probabilistic models to detect motifs characterized by unusually high number of copies of the motif in the analyzed sequences.  we present a novel method for detection of pairs of motifs separated by spacers of variable nucleotide sequence but conserved length. unlike existing methods for motif discovery, the motifs themselves are not required to occur at unusually high frequency but only to exhibit a significant preference to occur at a specific distance from each other. in the present implementation of the method, motifs are represented by pentamers and all pairs of pentamers are evaluated for statistically significant preference for a specific distance. an important step of the algorithm eliminates motif pairs where the spacers separating the two motifs exhibit a high degree of sequence similarity; such motif pairs likely arise from duplications of the whole segment including the motifs and the spacer rather than due to selective constraints indicative of a functional importance of the motif pair. the method was used to scan complete prokaryotic genomes for novel sequence motifs. some motifs detected were previously known but other motifs found in the search appear to be novel. selected motif pairs were subjected to further investigation and in some cases their possible biological functions were proposed.  we present a new motif-finding technique that is applicable to scanning complete genomes for sequence motifs. the results from analysis of genomes suggest that the method detects previously known motifs that are expected to be found as well as new motifs that are unlikely to be discovered by traditional motif-finding methods. we conclude that our approach to detection of significant motif pairs can complement existing motif-finding techniques in discovery of novel functional sequence motifs in complete genomes. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords genomesequence motifsmotif-findingbacteriaarchaeadna sequence repeatshttp://dx.doi.org/ /100000153division of biological infrastructuredbi-0950266issue-copyright-statement© the author 2017 