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
