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
