Vectores de sufijos para buscar patrones (suffix arrays)

Buenas,
hace poco tuve la suerte de cursar un taller práctico sobre problemas en el análisis, mapeo y ensamblado de secuencias genómicas, impartido por Paolo Ribeca (autor del software GEM). A parte de volver un poco abrumado por las dificultades que encierran los datos producidos por estas tecnologías de secuenciación de gran capacidad, sobre todo para bichos que no tienen todavía un genoma de referencia de buena calidad, he aprendido qué estructuras de datos soportan el manejo y análisis de todas estas secuencias. Además de tablas hash y grafos de De Bruijn, y la magia de Burrows-Wheeler, que se merecerá su propia entrada, hemos aprendido las ventajas de usar vectores de sufijos (suffix arrays) para el problema de buscar patrones o subcadenas exactas en textos muy largos, como son las secuencias genómicas. Pero, qué es un vector de sufijos? Trataré de explicarlo con el siguiente código en Perl:

 # calcula un vector de sufijos de manera que las diferentes subcadenas   
 # de una secuencia queden preordenadas lexicograficamente  
 sub make_suffix_array  
 {  
    my ($seq,$verbose) = @_;  
    my @suffix = (0 .. length($seq)-1); # en base 0  
   
    # ordena los length($seq) posibles sufijos lexicogr´aficamente 
    @suffix = sort {substr($seq,$a) cmp substr($seq,$b)} (@suffix);  
      
    if($verbose)  
    {  
       print "# suffix array for $seq :\n";  
       foreach my $suf (@suffix)  
       {   
          printf("%3d %s\n",$suf,substr($seq,$suf));   
       }  
       print "\n";  
    }  
    return @suffix;  
 }  

Si llamamos a la subrutina

my @suffix = make_suffix_array('TTTTAGATCGATCGACTAGACTACGACTCGA',1);

obtendremos:

30 A
22 ACGACTCGA
19 ACTACGACTCGA
14 ACTAGACTACGACTCGA
25 ACTCGA
17 AGACTACGACTCGA
4  AGATCGATCGACTAGACTACGACTCGA
10 ATCGACTAGACTACGACTCGA
6  ATCGATCGACTAGACTACGACTCGA
28 CGA
12 CGACTAGACTACGACTCGA
23 CGACTCGA
8  CGATCGACTAGACTACGACTCGA
20 CTACGACTCGA
15 CTAGACTACGACTCGA
26 CTCGA
29 GA
18 GACTACGACTCGA
13 GACTAGACTACGACTCGA
24 GACTCGA
9  GATCGACTAGACTACGACTCGA
5  GATCGATCGACTAGACTACGACTCGA
21 TACGACTCGA
16 TAGACTACGACTCGA
3  TAGATCGATCGACTAGACTACGACTCGA
27 TCGA
11 TCGACTAGACTACGACTCGA
7  TCGATCGACTAGACTACGACTCGA
2  TTAGATCGATCGACTAGACTACGACTCGA
1  TTTAGATCGATCGACTAGACTACGACTCGA
0  TTTTAGATCGATCGACTAGACTACGACTCGA

que es una lista de los 31 posibles sufijos de la cadena original, en orden lexicográfico. De hecho, si te fijas en el código, el vector realmente contiene sólo las posiciones (en base 0) de los sufijos ordenados, no su secuencia. Obviamente la construcción de este vector es costosa al necesitar de una ordenación (en Perl es por defecto un mergesort con un caso peor O(NlogN) ), pero luego permite consultas  O(logN), es decir, mucho más rápidas que simplemente recorrer la secuencia de principio a fin, al usar implícitamente un árbol binario.

Puedes probarlo mediante el siguiente ejemplo:

 match_pattern('TAG','TTTTAGATCGATCGACTAGACTACGACTCGA'); 

 # aprovecha el orden del vector de sufijos para buscar cualquier subcadena   
 # o patro´n (pattern) por medio de un árbol binario   
 sub match_pattern   
 {  
    my ($pattern,$seq) = @_;  
    print "# looking for pattern $pattern in sequence $seq (base 0)\n";  
    my @suffix = make_suffix_array($seq,1);  
    my $low = 0;  
    my $high = $#suffix;  
    my $patl = length($pattern);  
    my ($med,$submed);  
    while ($low <= $high)   
    {  
       my $med = int (($low+$high)/2); # punto medio de la búsqueda        
       # comparacion lexicográfica en punto medio, mira 'perldoc perlop'  
       my $comp = $pattern cmp substr($seq,$suffix[$med],length($pattern));  
       if($comp < 0){ $high = $med-1 }  # retrocedemos en @suffix   
       elsif($comp > 0){ $low = $med+1 } # avanzamos  
       else   
       {  
          my $submed = $med - 1; # sufijo inmediatamente anterior al punto medio  
          while($submed > 0 && $pattern eq substr($seq,$suffix[$submed],$patl))  
          {   
             print "# match at position $suffix[$submed]\n";  
             $submed--;   
          }  
          while ($med < $#suffix-1 && $pattern eq substr ($seq,$suffix[$med],$patl))  
          {  
             print "# match at position $suffix[$med]\n";  
             $med++;   
          }  
          last;  
       }  
    }  
 }  

Hasta otra,
Bruno