Abstract:
Methods, systems, and computer-readable media are provided which can identify and provide local variations in regions of similarity among two or more data sets. These data sets may be represented as sequences such as, e.g., genomic sequences or words in a text. The local variations in similarity levels can be provided by selecting an initial prior distribution relating the data sets, organizing the first data set into windows and the remaining data sets into blocks, using the priors to sample one or more sets of words from the first data set, computing a similarity curve from exact and inexact matches for these words and, if convergence of results is not achieved, computing a new set of priors and repeating the sampling and computation of similarity curves. The computations can be performed using an amount of computational time that is linearly proportional to the size of the data sets.

Description:
CROSS-REFERENCE TO RELATED APPLICATION(S) 
     This application claims priority from U.S. Patent Application No. 60/674,009, filed Apr. 22, 2005, the entire disclosure of which is incorporated herein by reference. 
    
    
     GOVERNMENTAL SUPPORT 
     The research leading to the present invention was supported, at least in part, by the National Science Foundation under Grant number CCR-0325605, and by DARPA under Grant no. F30602-01-2-0556. Thus, the U.S. government may have certain rights in the invention. 
    
    
     FIELD OF THE INVENTION 
     The present invention relates to an efficient method, system and software arrangement for comparing two or more long strings of characters with “almost common” subsequences interspersed. Examples of such strings of characters may include genome sequences of nucleotides, shorter sequences of amino acids, natural-language texts, natural or computer-generated texts or hypertexts, weblogs or “blogs,” network logs, scripts or object code encoding computer procedures that may have been successively revised, plagiarized texts, etc. 
     BACKGROUND INFORMATION 
     Comparisons of large sets of data to detect regions of similarity among them has many practical applications. For example, natural-language word sequences appearing in stories, articles, term papers, etc., may be compared to determine whether material has been copied or plagiarized from one to another. Computer code listings can be compared to assess whether or not code in one listing has been copied from another listing. In the fields of biology and bioinformatics, sequences of genes, nucleotides, amino acids, etc., that appear in different biological structures may be compared to detect if homology exists that can indicate similar biological functions or that the sequences are related to a common origin. 
     Sequence alignment tools have incorporated many innovations to greatly advance detailed comparative genomics studies. For example, recent segmental duplications in mammalian genomes (with identity level &gt;90%) can be detected using BLAST or other sequence alignment tools. 
     However, many of these tools used in genomic applications may use exact or inexact k-mers as homology seeds for local alignment extension. Thus there may be a compromise between sensitivity and computational efficiency when they are used to detect homologous segments, segmental duplications, or homology-based phylogenic distances, as homology levels become lower. 
     To improve sensitivity, many of these tools rely on exhaustive searches of exact matches with short mers or inexact matches with longer mers. This approach can yield a large number of false-positives, which may require subsequent post-processing filtering, which can be expensive. Alternatively, more stringent search criteria (e.g., longer mers with more exact matches) can be used to improve efficiency. However, these algorithms may fail to detect low-homology regions such as, for example, ancient duplication events. In order to detect less-recent duplications, orthologous genes have been used as “anchors” to map out duplication blocks. However, this approach may not be suitable for identifying duplications that are not subject to a strong selection process, e.g., sequences containing only non-coding regions. 
     SUMMARY OF THE INVENTION 
     One object of the present invention is to provide a system, a method, a software arrangement, and a computer-accessible medium (hereinafter collectively referred to as “PRIZM”) for efficiently and robustly detecting or identifying variations in similarity among two or more sets of data, which may be in the form of sequences. Similar regions between two sequences can be efficiently determined even when the similarity level drops to a value around 65%. In genomics applications, algorithms may be employed that incorporate background knowledge about genome evolution. For example, information relating to expected statistical distributions of similarity variations among the sets of data can be used. This information can be domain-dependent, and it can be based on the nature of the data sets being analyzed, or on a statistical model. This information may be in the form of a statistical prior distribution, a hierarchal prior model, a uniform distribution, a similarity curve, or a homology curve. A variety of priors may be used such as, e.g., a noninformative improper prior, an exponential prior, a Gamma prior, or a prior based on a Juke-Cantor one-parameter evolution model or a Kimura two-parameter evolution model. These priors may reflect statistical models based on, e.g., mutations or edit operations. The results obtained may be generally unaffected by the choice of a particular prior, while the computational efficiency may be only mildly affected by this choice. The time required to perform the computations can be linearly proportional to the size of the data sets. 
     In certain exemplary embodiments of the present invention, PRIZM can process multiple sequences, and the similarity regions between a given sequence and a set of other sequences may be determined. A sequence of rearrangement events and a rough estimation of the time of these rearrangements may also be computed. In addition, exemplary embodiments of the present invention can be applied to analyze other domains by establishing appropriate priors. 
     In further exemplary embodiments of the present invention, a Bayesian scheme can be used to refine parameters and priors used to converge similarity results to a desired degree of accuracy. A boosting technique may be used in conjunction with the Bayesian scheme. 
     In exemplary embodiments of the present invention, regions of predetermined similarity and/or local alignment data may be identified with respect to two or more data sets. 
     In still further exemplary embodiments of the present invention, similarity may be detected in comparative genomics, genome sequences of nucleotides, shorter sequences of amino acids, natural-language texts, natural or computer-generated texts or hypertexts, network logs, scripts or object code encoding computer procedures that may have been successively revised, code listings for webpages, weblogs, sets of time-series data, speech data, image data, video data, plagiarized texts, etc. 
     In certain exemplary embodiments of the present invention, regions of predetermined similarity and/or variations in similarity across regions can be determined using exact and/or inexact matches. Inexact matches can include matches that include one error or two errors between compared regions or subregions of the sets of data. Subregions of the sets of data can be selected for comparison based on expected statistical distributions of similarity variations. Subregions can be organized as sets of blocks and/or windows, and the blocks can be stored in the form of, e.g., a hash table or a suffix array. 
     In exemplary embodiments of the present invention, an iterative process may be used to obtain accurate determinations of variations in similarity regions and/or regions of predetermined similarity among data sets. Results of a calculated comparison may be used to modify or refine parameters or models associated with information relating to expected statistical distributions of similarity variations, and the comparisons may again be computed using the modified or refined information. This iterative process can be continued until specified condition is met. The condition can include, e.g., a preselected number of iterations or total computational time, or a convergence of the similarity calculations to a uniform result to a preselected degree of accuracy. 
     These and other objects, features and advantages of the present invention will become apparent upon reading the following detailed description of embodiments of the invention, when taken in conjunction with the appended claims. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       Further objects, features and advantages of the invention will become apparent from the following detailed description taken in conjunction with the accompanying figures showing illustrative embodiments of the invention, in which: 
         FIG. 1  is a general flowchart of certain exemplary embodiments of the present invention; and 
         FIG. 2  is a schematic diagram of a system in accordance with certain exemplary embodiments of the present invention. 
     
    
    
     DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS 
     The method, system and software arrangement according to an exemplary embodiment of the present invention is capable of generating a function called a “similarity curve” based on two sequences. The value of this function at a location on a first sequence can be selected as the highest local similarity level of a short region straddling that location, when compared with a second sequence. Certain characteristics related to this curve such as, e.g., its structure, a distribution of similarities across the entire first sequence, a correlation between the curve and a composition of the underlying regions (e.g., coding vs. noncoding, CG-content, stability or flexibility), etc. can provide information describing how the two sequences are related in terms of, for example, evolutionary distances, patterns of conservation, and/or mechanisms of evolution and selection that have acted since their divergence. 
     A location can be provided in a first sequence, G 1 [α], and a short window W 1,α =G 1 [α, α+m−1] can be specified starting at α. This window can be compared with a word of equal length from the second sequence starting at G 2 [β], e.g., W 2,β =G 2 [β, β+m−1]. The similarity level may be defined for the locations G 1 [α] and G 2 [β] and a window as 
     
       
         
           
             
               
                 h 
                 
                   ( 
                   2 
                   ) 
                 
               
               ⁡ 
               
                 ( 
                 
                   
                     
                       G 
                       1 
                     
                     ⁡ 
                     
                       [ 
                       α 
                       ] 
                     
                   
                   , 
                   
                     
                       G 
                       2 
                     
                     ⁡ 
                     
                       [ 
                       β 
                       ] 
                     
                   
                   , 
                   m 
                 
                 ) 
               
             
             = 
             
               
                 1 
                 m 
               
               ⁢ 
               
                 
                   ∑ 
                   
                     γ 
                     = 
                     0 
                   
                   
                     m 
                     - 
                     1 
                   
                 
                 ⁢ 
                 
                   { 
                   
                     
                       
                         1 
                       
                       
                         
                           
                             
                               G 
                               1 
                             
                             ⁡ 
                             
                               [ 
                               
                                 α 
                                 + 
                                 γ 
                               
                               ] 
                             
                           
                           = 
                           
                             
                               G 
                               2 
                             
                             ⁡ 
                             
                               [ 
                               
                                 β 
                                 + 
                                 γ 
                               
                               ] 
                             
                           
                         
                       
                     
                     
                       
                         0 
                       
                       
                         otherwise 
                       
                     
                   
                 
               
             
           
         
       
     
     The highest similarity level for the sequence G 1  at position a as computed with respect to G 2  can be denoted by h(α), which may be expressed as: 
     
       
         
           
             
               h 
               ⁡ 
               
                 ( 
                 α 
                 ) 
               
             
             = 
             
               
                 
                   
                     max 
                   
                 
                 
                   
                     
                       1 
                       ≤ 
                       β 
                       &lt; 
                       
                         
                            
                           
                             G 
                             2 
                           
                            
                         
                         - 
                         m 
                       
                     
                   
                 
               
               ⁢ 
               
                 
                   h 
                   
                     ( 
                     2 
                     ) 
                   
                 
                 ⁡ 
                 
                   ( 
                   
                     
                       
                         G 
                         1 
                       
                       ⁡ 
                       
                         [ 
                         α 
                         ] 
                       
                     
                     , 
                     
                       
                         G 
                         2 
                       
                       ⁡ 
                       
                         [ 
                         β 
                         ] 
                       
                     
                     , 
                     m 
                   
                   ) 
                 
               
             
           
         
       
     
     The “similarity curve” for the first sequence, G 1  with respect to the second sequence G 2  may then be defined as:
 
 h:G   1   [α . . . α+m− 1]→[0,1]
 
α           h(α)

     Low-similarity alignments may be detected without a quadratic number of pair-wise inter-sequence inexact searches, even when verification of a low-homology match is relatively easy. That is, the correctness of a given similarity mapping between two sequences can be verified in a linear amount of time, even if the similarity levels are quite low. Thus, while a deterministic version with sub-quadratic complexity may not be achievable, a linear-time non-deterministic algorithm exists. 
     In exemplary embodiments of the present invention, a non-deterministic technique can be replaced with a probabilistic estimation scheme, where probabilistic distributions employed therein can be based on, e.g., biologically meaningful priors. This scheme can include repetition of the following steps until convergence is achieved:
         a) provide an estimated probabilistic distribution;   b) compute a linear-time verification; and   c) refine the probabilistic distribution for a successive iteration.       

     For example, a Bayesian scheme optionally combined with a boosting technique can be used to modify the probabilistic distributions provided in step (a) above from one iteration to the next. For each iteration, a set of words with a specific distribution may be selected from one sequence, and can be optimally partitioned into groups for “in silico experiments” involving an exact-match search, an inexact-match search with one error, an inexact match search with two errors, etc. These searches can be conducted efficiently over the second sequence, particularly if the second sequence has been preprocessed and stored in an efficient data structure (e.g., as a suffix array or a hash table). Based on these calculations, a Bayesian estimator can compute the local similarity levels for the second sequence, and use it to verify and refine the probabilistic distributions for the next iteration. In this manner, convergence to the true local similarity levels may be achieved after only a few iterations. 
     Exemplary embodiments of the present invention are described in more detail as follows: 
     Prespecified parameters b, B=|G 2 |/b, w, m, and N 0 , N 1 , . . . , N k  (k≦m) may be defined, where k may be preferably set equal to 2. A set of k+1 m-mers S 0 , . . . , S k  may be selected randomly (uniformly with replacement) from G 1 [α, α+w−1] such that
 
 |S   0   |=N   0   ,|S   1   |=N   1 , . . . , and  |S   k   |=N   k .
 
An m-mer can be understood to represent a word of length m, which may be selected from either sequence.
 
     A block in the second sequence of length b may be identified as B β ≡G 2 [β, β+b−1]. The parameter X i  (i=0 to k) may be defined as the number of corresponding m-mers S i  (i=0 to k) that match exactly (with one, two, etc. up to k mismatches, respectively) to an m-mer in G 2 [β, β+b−1]. The values of X i  for 0≦i≦k can thus represent binomially distributed random variables whose parameters depend on the similarity level h. 
     The local similarity may be estimated by the following robust estimators: 
               〈       h   |     X   0       ,     X   1     ,   …   ⁢           ,     X   k       〉     =       ∫   0   1     ⁢     h   ⁢           ⁢     p   ⁡     (       h   |     X   0       ,     X   1     ,   …   ⁢           ,     X   k       )       ⁢       ⅆ   h     .               
The mean, standard deviation and confidence of the similarity function over B β  may also be determined. The parameter β* can be defined as arg max β  mean(B β ). The similarity function can then be estimated at α by the value mean (B β *).
 
     To create a robust estimator, it may be preferable to select a prior for p(h). A simple choice for this would be a uniform prior such as, e.g., p(h)=1, 0≦h≦1, or a slightly modified uniform prior, e.g., p(h)=4/3, ¼&lt;h≦1. These exemplary priors are uninformative and may be less restrictive. Another improper prior that may be preferable is, e.g., p(h)=λExp[−λh], 0≦h≦1. Other parametric priors, which may have an evolutionary explanation, can be created by applying a general theorem. 
     Several criteria may be used to detect convergence of the algorithm, at which point the processing operations may be stopped. For example, the processing may be interrupted once a predefined level of detail is obtained. Alternatively, a predetermined amount of time or a predetermined number of iterations can be selected for application of the algorithm. Processing may also be halted if calculated similarity estimates differ by less than a predetermined amount from the similarity curve computed at the previous step. 
     Detailed local alignment may be obtained by using the local similarity information, the pre-processed data structure, and a local alignment procedure/tool. 
     In iteration i, an interval I with k similarity estimates may be expressed as:
 
           μ 1 ,σ 1           ,         μ 2 ,σ 2           , . . . ,         μ k ,σ k           

     An updated prior may be determined using a heuristic argument, based on a weighted average of the estimated parameters. Normal random variables x i , i=1 to k, may be defined as
 
 x   1 =μ 1 +σ 1   N (0,1)
 
 x   2 =μ 2 +σ 2   N (0,1)
 
. . .
 
 x   k =μ k +σ k   N (0,1).
 
     A weighted average of these variables may be written as follows:
 
 X   1 =ω 1   x   1 +ω 2   x   2 + . . . +ω k   x   k ,
 
where the ω j  parameter can be a normalized weighting factor such that ω 1 +ω 2 + . . . +ω k =1. Thus
 
     
       
         
           
             
               
                 
                   
                     X 
                     I 
                   
                   = 
                     
                   ⁢ 
                   
                     
                       
                         ∑ 
                         j 
                       
                       ⁢ 
                       
                         
                           ω 
                           j 
                         
                         ⁢ 
                         
                           μ 
                           j 
                         
                       
                     
                     + 
                     
                       
                         ∑ 
                         j 
                       
                       ⁢ 
                       
                         
                           ω 
                           j 
                         
                         ⁢ 
                         
                           σ 
                           j 
                         
                         ⁢ 
                         
                           N 
                           ⁡ 
                           
                             ( 
                             
                               0 
                               , 
                               1 
                             
                             ) 
                           
                         
                       
                     
                   
                 
               
             
             
               
                 
                   ≈ 
                     
                   ⁢ 
                   
                     N 
                     ⁡ 
                     
                       ( 
                       
                         
                           μ 
                           I 
                         
                         , 
                         
                           σ 
                           I 
                         
                       
                       ) 
                     
                   
                 
               
             
           
         
       
     
     If each weighting factor ω j  is selected to be inversely proportional to the standard deviation of the corresponding normal random variables x j , i.e., 
                 ω   j     =         ∑   j     ⁢       (     σ   j     )       -   1           σ   j         ,         
then it can be shown that
 
     
       
         
           
             
               μ 
               I 
             
             = 
             
               
                 
                   
                     
                       ∑ 
                       j 
                     
                     ⁢ 
                     
                       
                         
                           μ 
                           j 
                         
                         ⁡ 
                         
                           ( 
                           
                             σ 
                             j 
                           
                           ) 
                         
                       
                       
                         - 
                         1 
                       
                     
                   
                   
                     
                       ∑ 
                       j 
                     
                     ⁢ 
                     
                       
                         ( 
                         
                           σ 
                           j 
                         
                         ) 
                       
                       
                         - 
                         1 
                       
                     
                   
                 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 and 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   σ 
                   I 
                   2 
                 
               
               = 
               
                 
                   k 
                   2 
                 
                 
                   
                     ( 
                     
                       
                         ∑ 
                         j 
                       
                       ⁢ 
                       
                         
                           ( 
                           
                             σ 
                             j 
                           
                           ) 
                         
                         
                           - 
                           1 
                         
                       
                     
                     ) 
                   
                   2 
                 
               
             
           
         
       
     
     The distribution for an interval I may be approximated by p i,1 (h)≈N(μ 1 σ 1 ). The parameter p 0,G =p(h) may be set equal to the initial prior chosen above. 
     A prior may also be refined by the choice of a convex combination determined by a parameter 0&lt;δ≦1, e.g.:
 
 p   i,1 ( h )=(1−δ) p   0,G ( h )+δ(μ 1 +σ 1   N (0,1)).
 
     Parameters may be optimized in several ways. For example, the number of blocks (b) and the number of windows (w) may be chosen a priori based on a selected resolution for similarity. These parameters may be selected such that b=O(√{square root over (G 2 )}) and w=O(√{square root over (G 1 )}). An amount of work K=O(1) may be assigned to a region defined by a combination of any single block with any single window. Thus the amount of work may be expressed as approximately K(G 1 G 2 )/(wb)=O(G 1 +G 2 ) per iteration. 
     The mer size parameter m may be selected such that the probability of a match or “hit” in a block containing a similarity sequence is much higher than in a random block, i.e.,
 
( b/ 4 m )&lt;&lt; E ( h   0,G ) m .
 
     The similarity curve between a sequence S 1  and a set of n sequences T 1 , . . . , T n , may also be determined by applying the procedures described above to the sequence S 1  and a new sequence S 2  obtained by concatenating T 1 , . . . , T n . 
     An exemplary flow diagram of a method  100  according to exemplary embodiments of the present invention is shown in  FIG. 1 . Sequences S 1  and S 2  to S n  which are to be compared may be read from a database (step  110 ). A distribution of priors relating S 1  and S 2  to S n  may then be chosen (step  120 ). The sequences S 2  to S n  can be organized into blocks (step  130 ). The words corresponding to these blocks may be stored in a searchable data structure. The sequence S 1  can then be organized into windows (step  140 ). This can allow a set of words or mers to be sampled rapidly from each window. One or more sets of words from S 1  may then be sampled (step  150 ). Predetermined parameters may be used to choose the words used in this step. A similarity curve, which may be based on identification of exact and inexact matches of these words, can then be calculated with respect to S 2  to S n  (step  160 ). Priors may be used to determine the local similarity levels. The desired convergence of the process may then be evaluated (step  170 ). If convergence is not achieved, then the priors may then be updated based on the calculated results (step  180 ). Further sets of words from S 1  may then be sampled (step  150 ), and exact and inexact matches can again be calculated (step  160 ). If convergence is achieved, local similarity levels between the sequences S 1  and S 2  to S n  may be provided as output (step  190 ). 
     An exemplary embodiment of a system for implementing the present invention is shown in  FIG. 2 . For example, a database  210  may contain two or more sequences to be compared. These sequences may include genetic sequences, text strings, computer code, etc. The system may further comprise and/or utilize external information  220 , which can include, e.g., domain-dependent knowledge or statistical models that relates to the sequences that are to be compared, as well as algorithms that may be used to conduct searches of the sequences to detect various types of matches. The external information  220  and information contained in the database  210  can be communicated to a computer  230 . The computer  230  may be configured, e.g., to calculate matches and/or local regions of similarity between the sequences. Similarity regions and other information relating to the sequences being compared may be communicated to an output device  240 . The output device  240  may comprise any suitable display device, including but not limited to a video monitor, a printer, a data storage medium, and the like. The computer  230  can include a hard drive, CD ROM, RAM, and/or other storage devices or media which can include thereon software, which can be configured to execute the exemplary embodiments of the method of the present invention. 
     The foregoing merely illustrates the principles of the invention. Various modifications and alterations to the described embodiments will be apparent to those skilled in the art in view of the teachings herein. It will thus be appreciated that those skilled in the art will be able to devise numerous systems, arrangements and methods which, although not explicitly shown or described herein, embody the principles of the invention and are thus within the spirit and scope of the present invention. In addition, all publications referenced above are incorporated herein by reference in their entireties.