En bioinformática (y otros muchos campos), una de las tareas habituales es la búsqueda de una cadena en otra más larga. Saber si está, dónde está, y cuántas más veces aparece.
Aún con los ordenadores más potentes, la simple búsqueda secuencial tiene límites si la cadena donde buscamos es enorme.
Una de las técnicas que se suelen utilizar, para acelerar el proceso de búsqueda, es la llamada
Array de sufijos.
Supongamos que tenemos una secuencia así:
AGTGTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
entonces, podemos deducir sus sufijos:
Using text Syntax Highlighting
AGTGTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
GTGTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
TGTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
GTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
TTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
TCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
CCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
CTCCATAGATATGGCGGTTACCTCCGTTCGTGG
TCCATAGATATGGCGGTTACCTCCGTTCGTGG
CCATAGATATGGCGGTTACCTCCGTTCGTGG
CATAGATATGGCGGTTACCTCCGTTCGTGG
ATAGATATGGCGGTTACCTCCGTTCGTGG
TAGATATGGCGGTTACCTCCGTTCGTGG
AGATATGGCGGTTACCTCCGTTCGTGG
GATATGGCGGTTACCTCCGTTCGTGG
ATATGGCGGTTACCTCCGTTCGTGG
TATGGCGGTTACCTCCGTTCGTGG
ATGGCGGTTACCTCCGTTCGTGG
TGGCGGTTACCTCCGTTCGTGG
GGCGGTTACCTCCGTTCGTGG
GCGGTTACCTCCGTTCGTGG
CGGTTACCTCCGTTCGTGG
GGTTACCTCCGTTCGTGG
GTTACCTCCGTTCGTGG
TTACCTCCGTTCGTGG
TACCTCCGTTCGTGG
ACCTCCGTTCGTGG
CCTCCGTTCGTGG
CTCCGTTCGTGG
TCCGTTCGTGG
CCGTTCGTGG
CGTTCGTGG
GTTCGTGG
TTCGTGG
TCGTGG
CGTGG
GTGG
TGG
GG
G
Coloreado en 0.000 segundos, usando
GeSHi 1.0.8.4
Ahora, estas cadenas las ordenamos alfabéticamente:
Using text Syntax Highlighting
ACCTCCGTTCGTGG
AGATATGGCGGTTACCTCCGTTCGTGG
AGTGTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
ATAGATATGGCGGTTACCTCCGTTCGTGG
ATATGGCGGTTACCTCCGTTCGTGG
ATGGCGGTTACCTCCGTTCGTGG
CATAGATATGGCGGTTACCTCCGTTCGTGG
CCATAGATATGGCGGTTACCTCCGTTCGTGG
CCGTTCGTGG
CCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
CCTCCGTTCGTGG
CGGTTACCTCCGTTCGTGG
CGTGG
CGTTCGTGG
CTCCATAGATATGGCGGTTACCTCCGTTCGTGG
CTCCGTTCGTGG
G
GATATGGCGGTTACCTCCGTTCGTGG
GCGGTTACCTCCGTTCGTGG
GG
GGCGGTTACCTCCGTTCGTGG
GGTTACCTCCGTTCGTGG
GTGG
GTGTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
GTTACCTCCGTTCGTGG
GTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
GTTCGTGG
TACCTCCGTTCGTGG
TAGATATGGCGGTTACCTCCGTTCGTGG
TATGGCGGTTACCTCCGTTCGTGG
TCCATAGATATGGCGGTTACCTCCGTTCGTGG
TCCGTTCGTGG
TCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
TCGTGG
TGG
TGGCGGTTACCTCCGTTCGTGG
TGTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
TTACCTCCGTTCGTGG
TTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
TTCGTGG
Coloreado en 0.000 segundos, usando
GeSHi 1.0.8.4
Ya podemos realizar la búsqueda de una subcadena. Supongamos que queremos buscar la cadena GTT. Por medio de búsqueda binaria, localizamos aquellas cadenas que comienzan por esa cadena. Y el número de cadena que encontramos coincide con la posición dentro de la cadena original (si GTT lo encontramos en la subcadena 24, esa es la misma posición donde se encuentra dentro de la primera cadena).
El problema ahora viene en la parte de ordenar los sufijos. Si la cadena es muy larga, puede tardar mucho el generar todas las subcadenas, aparte de consumir mucha memoria, y luego, hacer el proceso de ordenación.
Por ello, lo que se hace es ordenar un vector con las posiciones de las subcadenas. No hace falta extraerlas de la cadena original. Solo es necesario guardar las posiciones de comienzo. Y hacer las ordenaciones usando esas posiciones como referencia.
Así, para el ejemplo mostrado, quedaría:
Using text Syntax Highlighting
27 ACCTCCGTTCGTGG
14 AGATATGGCGGTTACCTCCGTTCGTGG
1 AGTGTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
12 ATAGATATGGCGGTTACCTCCGTTCGTGG
16 ATATGGCGGTTACCTCCGTTCGTGG
18 ATGGCGGTTACCTCCGTTCGTGG
11 CATAGATATGGCGGTTACCTCCGTTCGTGG
10 CCATAGATATGGCGGTTACCTCCGTTCGTGG
31 CCGTTCGTGG
7 CCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
28 CCTCCGTTCGTGG
22 CGGTTACCTCCGTTCGTGG
36 CGTGG
32 CGTTCGTGG
8 CTCCATAGATATGGCGGTTACCTCCGTTCGTGG
29 CTCCGTTCGTGG
40 G
15 GATATGGCGGTTACCTCCGTTCGTGG
21 GCGGTTACCTCCGTTCGTGG
39 GG
20 GGCGGTTACCTCCGTTCGTGG
23 GGTTACCTCCGTTCGTGG
37 GTGG
2 GTGTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
24 GTTACCTCCGTTCGTGG
4 GTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
33 GTTCGTGG
26 TACCTCCGTTCGTGG
13 TAGATATGGCGGTTACCTCCGTTCGTGG
17 TATGGCGGTTACCTCCGTTCGTGG
9 TCCATAGATATGGCGGTTACCTCCGTTCGTGG
30 TCCGTTCGTGG
6 TCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
35 TCGTGG
38 TGG
19 TGGCGGTTACCTCCGTTCGTGG
3 TGTTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
25 TTACCTCCGTTCGTGG
5 TTCCTCCATAGATATGGCGGTTACCTCCGTTCGTGG
34 TTCGTGG
Coloreado en 0.000 segundos, usando
GeSHi 1.0.8.4
Solo almacenamos los números, y con ellos hacemos las ordenaciones y las posteriores búsquedas. Bueno, en realidad es algo más complicado, ya que también se tiene en cuenta la concordancia de caracteres entre sufijos vecinos, que aceleran la búsqueda (se comenta en la página de Wikipedia).
El concepto de
array de sufijos, así como su creación, ordenación, y búsqueda, es obra de Gene Myers y Udi Manber, publicado en mayo de 1989, revisado en agosto de 1991, y publicado en la revista SIAM J. Computing en octubre de 1993. La mejor descripción se encuentra en
la publicación de agosto de 1991 (500Kb, PDF).
Hace unos días, nuestros amigos del sitio
#!/perl/bioinfo publicaron
un artículo al respecto, con una implementación del algoritmo de ordenación Manber/Myers escrita en Python.
Es obvio que una solución escrita en Perl, dará unos resultados similares a los ofrecidos con Python, pero los dos quedan muy lejos de los tiempos que una
implementación en C ofrece. Ni Python ni Perl se deberían usar para resolver estos problemas, comparado con las prestaciones que nos da un programa compilado. Pero sí que es más fácil de ver y leer el funcionamiento interno del algoritmo. Y por si alguien se pregunta si usar el comando
sort ofrecido por el propio lenguaje Perl o Python, no es más rápido que todo esto, le remito al artículo enlazado antes, donde hay una tabla con tiempos entre el algoritmo Manber/Myers y el Mergesort, que es el que suele usar sort().
Como este es un sitio dedicado a Perl, pues podemos intentar mostrar distintas soluciones de ordenación, también a título formativo.
Primero, cómo sería la generación y ordenación de sufijos, en puro Perl, usando las funciones intrínsecas del lenguaje:
Using perl Syntax Highlighting
#!/usr/bin/perl
# Prueba de generación de un array de sufijos ordenados usando solo sort()
# by Joaquin Ferrero, 2012.
#
### Módulos
use Modern::Perl;
use utf8::all;
use autodie;
use File::Slurp;
use List::MoreUtils qw(uniq);
### Constantes
use constant {
DEBUG => 1,
};
### Variables
my $A; # Secuencia a procesar
my $N; # Longitud de la secuencia
my @Σ; # Alfabeto
### Argumentos
if (@ARGV) {
my $arg = shift;
if ($arg =~ /^\d+$/) { # Si pasamos un número
@Σ = sort qw(A C G T); # El alfabeto son las cuatro bases de siempre
$A = join '', map { $Σ[rand 4] } 1 .. $arg; # Generamos la secuencia
$N = $arg; # La longitud de la secuencia es la indicada
}
else {
$A = read_file($arg, chomp => 'Yes'); # Leemos el archivo
@Σ = sort {$a cmp $b} uniq split //, $A; # Sacamos su alfabeto
$N = length $A; # Y su longitud
die "ERROR: El fichero está vacío\n" if not $N;
}
}
else {
die "Uso: $0 <longitud secuencia aleatorio|fichero con secuencia>\n";
}
if (DEBUG) {
say $A if $N < 50;
say "Longitud: $N";
say "Alfabeto: @Σ";
}
# Ordenación usando mergesort de Perl
my @sufijos = sort { substr($A, $a) cmp substr($A, $b) } 0 .. $N-1;
if ($N < 50) { # Si la secuencia es pequeña
for (@sufijos) { # Sacamos los sufijos ordenados
printf "%6d\t%s\n", $_, substr $A, $_;
}
}
__END__
Coloreado en 0.003 segundos, usando
GeSHi 1.0.8.4
Para secuencias pequeñas, los tiempos son cortos (todos los tiempos sacados en una máquina AMD Athlon(tm) 64 Processor 3000+ a 2Ghz):
Using text Syntax Highlighting
Long. Tiempo
5000 0m00.241s
10000 0m00.516s
20000 0m01.378s
40000 0m05.636s
80000 0m26.524s
160000 5m40.411s
Coloreado en 0.000 segundos, usando
GeSHi 1.0.8.4
Ahora veamos cómo sería el código de una implementación exacta del algoritmo de ordenación Manber/Myers:
Using perl Syntax Highlighting
#!/usr/bin/perl
# Simple algorithm for Suffix Array creation,
# based on Manber & Myers, Aug 1991
# by Joaquin Ferrero, 2012.
#
# Esta es una versión idéntica a la versión original de Myers, en C,
# con propósito únicamente formativo.
#
# Estadísticas para AMD Athlon 64 3000+ a 2Ghz:
#
# Longitud Tiempos
# ---------+----------
# 5000 0m00.794s
# 10000 0m01.196s
# 20000 0m02.271s
# 40000 0m04.506s
# 80000 0m09.685s
# 160000 0m19.470s
# 320000 0m43.474s
# 640000 1m27.284s
#
### Módulos
use Modern::Perl;
use utf8::all;
use autodie;
use File::Slurp;
use List::MoreUtils qw(uniq);
### Constantes
use constant {
DEBUG => 1,
BIT => 0x8000_0000,
MASK => 0x7FFF_FFFF,
};
### Variables
my $A; # Secuencia a procesar
my $N; # Longitud de la secuencia
my @Σ; # Alfabeto
### Argumentos
if (@ARGV) {
my $arg = shift;
if ($arg =~ /^\d+$/) { # Si pasamos un número...
@Σ = sort qw(A C G T); # con este alfabeto...
$A = join '', map { $Σ[rand 4] } 1 .. $arg; # generamos la secuencia...
$N = $arg; # de esta longitud.
}
else {
$A = read_file($arg, chomp => 'Yes'); # Si no, leemos el archivo indicado
@Σ = sort {$a cmp $b} uniq split //, $A; # del que extraemos su alfabeto
$N = length $A; # y su longitud
die "ERROR: El fichero está vacío\n" if not $N;
}
}
else {
die "Uso: $0 <longitud secuencia aleatorio|fichero con secuencia>\n";
}
if (DEBUG) {
say $A if $N < 50;
say "Longitud: $N";
say "Alfabeto: @Σ";
}
### Reserva de espacio
my @Pos; # Índice -> Inicio de sufijo
my @Prm; # Inicio de sufijo -> Índice de @Pos (Inversa de la anterior)
my @Lcp; # Longest Common Prefixes
my @BH;
my @B2H;
my @Min; # auxiliares
my @Mid;
$#Pos = $N+1 -1; # Reserva de espacio en memoria
$#Lcp = $N-1 -1; # -1 porque $#var es el último índice de @var,
$#Prm = $N+1 -1; # no el número de elementos de @var
$#BH = $N+1 -1;
$#B2H = $N+1 -1;
my $lgN = 0; # log N / log 2
for (my $i = 1; $i < $N-1; $i <<= 1) { $lgN++ } # averigüamos el número de vueltas del bucle principal
if ($lgN > 0) {
$#Min = $lgN -1;
$#Mid = $lgN -1;
}
### Inicialización
$A .= "\0"; # Terminador
for my $a (0 .. $N-1-1) {
$Lcp[$a] = $N + 1; # Valores iniciales del LCP
}
### Firt phase bucket sort # Primera ordenación basada en las letras del alfabeto
my %Bucket; # Primera aparición de cada letra del alfabeto
my @Link;
my($a, $b);
for my $c ( @Σ ) {
$Bucket{$c} = -1; # Los bucket iniciales lo forman todas las letras del alfabeto
}
for ($a = $N-1; $a >= 0; $a--) { # Ordenación de los buckets
my $c = substr $A, $a, 1;
$Link[$a] = $Bucket{$c};
$Bucket{$c} = $a;
}
my $rank = 0;
my $c = 0;
for my $a (@Σ) { # Sigue la ordenación, a lo que sumamos el cálculo de $rank
my $i = $Bucket{$a};
while ($i != -1) {
my $j = $Link[$i];
$Prm[$i] = $c;
if ($i == $Bucket{$a}) {
$BH[$c] = 1;
$rank++; # es necesario una vuelta más en el bucle principal
set($c, 0) if $c;
}
else {
$BH[$c] = 0;
}
$c++;
$i = $j;
}
}
$BH[$N] = 1;
$Pos[$N] = $Prm[$N] = 0;
for my $i (0 .. $N-1) { # Primeros valores de los array de ordenación
$Pos[$Prm[$i]] = $i;
}
# Inductive sorting stage
for (my $H = 1; $H < $N and $rank < $N; $H <<= 1) { # Bucle principal
say "Level $H. rank $rank" if DEBUG;
my($c, $d);
for ( $a = 0; $a < $N; $BH[$b] = 0 ) {
for ($b = $a++; !$BH[$a]; $a++) {
$Prm[$Pos[$a]] = $b;
}
}
for ($c = $N-1; $c >= $N-$H; $c--) {
$b = $Prm[$c];
for ($d = $b; $Pos[$d] != $c; $d++) { }
$Pos[$d] = $Pos[$b];
$Pos[$b] = $c;
$BH[$b] = 1;
$Prm[$c] = 1;
}
for $a (0 .. $N-1) {
if (($c = $Pos[$a] - $H) >= 0) {
$b = $Prm[$c];
if ($BH[$b]) {
$Prm[$c] += $Prm[$Pos[$b]]++;
}
else {
for ($d = $b; $Pos[$d] != $c; $d++) { }
$Pos[$d] = $Pos[$b];
$Pos[$b] = $c;
$BH[$b] = 1;
$Prm[$c] = 1;
}
}
}
for $a (0 .. $N-1) {
if ($BH[$a]) {
$Prm[$Pos[$a]] = $a;
}
}
for ($a = 0; $a < $N; ) {
$b = $a;
do {
$c = $Pos[$a++] - $H;
if ($c >= 0) {
$B2H[$Prm[$c]] = 1;
}
} while ( ! $BH[$a] );
$a = $b;
do {
$c = $Pos[$a++] - $H;
if ($c >= 0) {
$c = $Prm[$c];
if ($B2H[$c]) {
for ($c++; $B2H[$c] and !$BH[$c]; $c++) {
$B2H[$c] = 0;
}
}
}
} while ( !$BH[$a] );
}
for $a (0 .. $N-1) {
$Pos[$Prm[$a]] = $a;
}
for $a (0 .. $N-1) {
if ($B2H[$a]) {
$B2H[$a] = 0;
if ( !$BH[$a] ) {
$c = lcp($Pos[$a-1] + $H, $Pos[$a] + $H);
$rank++;
set($a, $H + $c);
$BH[$a] = 1;
}
}
}
}
show() if $N < 50;
#save();
exit 0;
### Subrutinas
# Cálculo de las coincidencias entre sufijos vecinos
sub lcp {
my($a, $b) = @_;
my($min, $opt, $cp, $L, $M, $R);
return 0 if substr($A, $a, 1) ne substr($A, $b, 1);
$a = $Prm[$a];
$b = $Prm[$b];
if ($a > $b) {
($a, $b) = ($b, $a);
}
return $Lcp[0] if $a == 0 and $b == $N-1;
$min = $N+1;
if ($a > 0) {
$L = 0;
$R = $N-1;
$opt = $Lcp[0];
while ($a != $R) {
$M = ($L+$R) >> 1;
if ($a <= $M) {
if ($Lcp[$M] & BIT) {
$cp = $opt;
$opt = $Lcp[$M] & MASK;
}
else {
$cp = $Lcp[$M];
}
if ($R <= $b and $cp < $min) {
$min = $cp;
}
$R = $M;
}
else {
if ( ! ($Lcp[$M] & BIT)) {
$opt = $Lcp[$M];
}
$L = $M;
}
}
}
if ($b < $N-1) {
$L = 0;
$R = $N-1;
$opt = $Lcp[0];
while ($b != $L) {
$M = ($L+$R) >> 1;
if ($b >= $M) {
if ($Lcp[$M] & BIT) {
$cp = $Lcp[$M] & MASK;
}
else {
$cp = $opt;
$opt = $Lcp[$M];
}
if ($a <= $L and $cp < $min) {
$min = $cp;
}
$L = $M;
}
else {
if ($Lcp[$M] & BIT) {
$opt = $Lcp[$M] & MASK;
}
$R = $M;
}
}
}
return $min;
}
# Actualiza los valores del árbol de sufijos
sub set {
my($a, $val) = @_;
my($L, $M, $R, $opt, $h);
$L = 0;
$R = $N-1;
$h = 0;
$opt = $Lcp[0];
while ($R-$L > 1) {
$M = ($L+$R) >> 1;
$Min[$h] = $opt;
$Mid[$h] = $M;
if ($a <= $M) {
$R = $M;
if ($Lcp[$M] & BIT) {
$opt = $Lcp[$M] & MASK;
}
}
else {
$L = $M;
if ( ! ($Lcp[$M] & BIT)) {
$opt = $Lcp[$M];
}
}
$h++;
}
for ($h--; $h >= 0; $h--) {
$opt = $Min[$h];
$M = $Mid[$h];
if ($Lcp[$M] & BIT) {
$L = $Lcp[$M] & MASK;
$R = $opt;
}
else {
$L = $opt;
$R = $Lcp[$M];
}
if ($a <= $M) {
if ($val < $L) {
$L = $val;
}
}
else {
if ($val < $R) {
$R = $val;
}
}
if ($L < $R) {
$Lcp[$M] = $R;
}
else {
$Lcp[$M] = $L | BIT;
}
}
if ($val < $Lcp[0]) {
$Lcp[0] = $val;
}
}
# Muestra el resultado en pantalla
sub show {
for my $i (0 .. $N-1) {
printf "%6d: %6d %6d ", $i, $Prm[$i], $Pos[$i];
if ($BH[$i]) {
if ($i > 0) {
printf "[%6d]+ '", lcp($Pos[$i-1], $Pos[$i] );
}
else {
printf "[%6d]+ '", 0;
}
}
else {
printf "%9s '", '';
}
print substr $A, $Pos[$i], 25;
print "'\n";
}
print "\n";
}
# Guarda el resultado en disco
sub save {
open my $fh, '>', 'out.pl.txt';
for my $i (0 .. $N-1) {
printf $fh "%d\t%s\n", $Pos[$i], substr $A, $Pos[$i], -1;
}
close $fh;
}
__END__
Coloreado en 0.007 segundos, usando
GeSHi 1.0.8.4
(Aquí también hay una subrutina que nos permite guardar a disco los sufijos ordenados, por si los necesitamos más adelante.) (No vamos a entrar en la descripción del funcionamiento del algoritmo. Consultar la publicación original de Manber y Myers.)Vemos que los tiempos son mucho menores que con solo el mergesort de nuestra función sort(). Esa es la razón de por qué usar estos algoritmos y no confiar solo en las funciones incorporadas.
Pero, repetimos, esto es solo a título formativo. La solución que hay que poner en un entorno de producción es la de un programa compilado. Observad los tiempos que se obtienen con el algoritmo Manber/Myers escrito en C:
Using text Syntax Highlighting
Long. Tiempo
5000 0m00.014s
10000 0m00.029s
20000 0m00.059s
40000 0m00.127s
80000 0m00.209s
160000 0m00.501s
320000 0m01.190s
640000 0m02.695s
1280000 0m05.782s
Coloreado en 0.000 segundos, usando
GeSHi 1.0.8.4
Si, aun así, nos empeñamos en querer usar Perl
, lo que hay que hacer es usar el módulo
Inline::C, e implementar el algoritmo en C dentro del código Perl. Perl lo compilará en una biblioteca externa, y lo ejecutará a la misma velocidad que si fuera el programa compilado en puro C.
Comentar, finalmente, que desde Manber y Myers publicaron este algoritmo, se han desarrollado una veintena más,
algunos de ellos hasta 30 veces más rápidos.
