13 de enero de 2012

Construir un suffix array: Manber y Myers

Hola,
hoy me gustaría volver un poco sobre el tema de los suffix array que ya ha sido tratado alguna vez en el blog.

En 1990, Manber y Myers publicaban Suffix arrays: a new method for on-line string searches, como motivo del primer ACM-SIAM symposium on Discrete algorithms. En éste fantástico e ilustrativo artículo, introducen la comparación de esta estructura con los suffix trees y proponen definiciones, notación y algoritmos para su construcción, búsqueda y optimización.

La mayor desventaja de suffix trees reside en el espacio requerido, y con los suffix arrays se consigue reducir considerablemente, siendo aún así competitivos frente a los suffix trees en términos de búsqueda y construcción.

Aquí veremos la primera parte, la creación de un suffix array ordenado (sección 3 del artículo original, sorprendentemente), con un algoritmo que requiere O(NlogN) como peor caso, si bien con optimizaciones se podría conseguir un tiempo esperado O(N); lo cual ya se consiguió de manera formal en 2003.

Quizás hoy día los suffix array van perdiendo protagonismo de primera plana, pero siguen siendo la base para muchas aplicaciones, por ejemplo implicados en métodos de obtención de la secuencia Barrows-Wheeler. Posiblemente en entradas posteriores hablemos de mejoras en el uso de suffix arrays y de otras estructuras relacionadas como FM-index y BWT.

Por ahora, aquí os dejo el código en Python (-V 2.6.6). Podéis cambiar la secuencia a indexar en la variable inicial A.

Y nos podemos preguntar ¿por qué usar éste código en lugar de la simple línea de ordenamiento MergeSort que aparecía en el código Perl de la entrada anterior sobre suffix arrays? Adjunto una tabla que lo explica en términos de tiempo requerido para la construcción del índice.


Long texto Myers MergeSort
20 0''019 0''004
100 0''023 0''005
1,000 0''055 0''013
11,000 0''717 0''401
110,000 22''200 41''075
220,000 1'36''004 2'42''189
660,000 9'31''297 25'54''598


Por supuesto estamos espectantes ante cualquier comentario, crítica o sugerencia.
¡Un saludo!

 #!/usr/bin/python

# Simple algorithm for Suffix Array creation,
# based on Manber & Myers, 1990
# by Cantalapiedra, 2011

############### FUNCTIONS ##############
## Function to print a Suffix Array
def printSA(text, sa = [], sa_name = "Suffix Array"):
  if len(sa) == 0:
       print "No records in {0}".format(sa_name)
  else:
       print sa_name
       for i in sa:
            print "{0}\t{1}".format(i, text[i:])
  print "\n"
  return

## Function to print one or several lists
def printTemp(things = [], title = "##", post_title = "#"):
  if len(things) == 0:
       print "No things to print"
  else:
       print title
       for a in things:
            print a
       print post_title
  print "\n"
  return

############### VARIABLES ##############

A = 'TTTTAGATCGATCGACTAGA$'

pos = [] # Contains the SA based on suffix numbers (SN)
prm = [] # Inverse of pos, each position is a SN containing its position in SA
bH = [] # Boolean array marking the buckets

b2H = [] # Temporary array for marking the 2H-buckets
count = [] # Counts the number of SN in a 2H-bucket

## All above variables have length = |A|

################ START ###################
# Some verbose
print "Your text is {0}\n".format(A[:len(A) - 1])

# Initial SA (POS)
pos = [(a,i) for i,a in enumerate(A)]
pos = sorted(pos)
pos = [a[1] for a in pos]

printSA(A, pos, "INITIAL SA")

# Initialize PRM
prm = [0]*len(A)

# Initial positions and buckets
prev_letter = ""
for i, pos_i in enumerate(pos):
  prm[pos_i] = i

  letter = A[pos_i]
  if letter != prev_letter:
       bH.append(1)
  else:
       bH.append(0)
  prev_letter = letter

######### H-STAGES ##########

count_changed = []

h_stage = 1
while h_stage < len(A): # max log2(|A| + 1)

  #print "stage {1} init {0}".format(time.time(), h_stage) IMPORT time module to use this

  # FIRST: re-assign prm as of buckets starting points
  #            and create count and b2H arrays
  prev_one = 0

  for i, bh_i in enumerate(bH):
       if bh_i == 1:
            prev_one = i # Is a bucket starting position
     
       if bh_i == 0:
            prm[pos[i]] = prev_one # Points to its bucket starting pos
     
       count.append(0)
       b2H.append(0)

  # Reset count[] to 0s
  for a in range(len(count)):
       count[a] = 0

  # THEN: for every bucket...
  # Here the entire SA is scanned, but in fact takes account
  # of the current bucket being covered
  prev_count = [a for a in count]
  for i, pos_i in enumerate(pos):
     
       if bH[i] == 1:
            b2H[i] = 1
     
       # We are moving the H-prefix (Ti) of suffix in i position
       ti = pos_i - h_stage
       if ti < 0: continue
     
       ### This is the heart of the algorithm
       # If its the first change of this bucket during this H-bucket run...
       if prev_count[prm[ti]] == count[prm[ti]]:
            change = True
       else:
            change = False
     
       count[prm[ti]] += 1
       count_changed.append(prm[ti])
       prm[ti] = prm[ti] + count[prm[ti]] - 1
     
       # ... that means this is a new 2H-bucket
       if change:
            b2H[prm[ti]] = 1
     
       #printTemp([prm, count, b2H], "After moving suff: {0}".format(ti), "")
     
       # If next position is another H-bucket, or is the last bucket
       if (i < len(bH) - 1 and bH[i + 1] == 1) or i == len(bH) - 1:
          
            #printTemp([prm, count, b2H], "", "NEXT H-BUCKET {0}".format(i + 1))
          
            for j in count_changed:
                 prev_count[j] = count[j]
            count_changed = []
     
       # If there are no buckets left
       if 0 not in b2H:
            break

  #print "end stage {0}".format(time.time())

  # If there are no buckets left
  if 0 not in b2H:
       break

  #printTemp([prm, count, b2H], "", "NEXT ROUND {0}".format(h_stage*2))

  # Update POS from PRM
  for j, prm_j in enumerate(prm):
       pos[prm_j] = j

  # Update BH from B2H
  bH = b2H

  # Reset vars for next H-STAGE
  count = []
  b2H = []
  h_stage *= 2

# Update POS from PRM
for j, prm_j in enumerate(prm):
  pos[prm_j] = j

printSA(A, pos, "FINAL SA")
#


La salida del ejemplo es:

Your text is TTTTAGATCGATCGACTAGA

INITIAL SA
$
AGATCGATCGACTAGA$
ATCGATCGACTAGA$
ATCGACTAGA$
ACTAGA$
AGA$
A$
CGATCGACTAGA$
CGACTAGA$
CTAGA$
GATCGATCGACTAGA$
GATCGACTAGA$
GACTAGA$
GA$
TTTTAGATCGATCGACTAGA$
TTTAGATCGATCGACTAGA$
TTAGATCGATCGACTAGA$
TAGATCGATCGACTAGA$
TCGATCGACTAGA$
TCGACTAGA$
TAGA$

FINAL SA
$
A$
ACTAGA$
AGA$
AGATCGATCGACTAGA$
ATCGACTAGA$
ATCGATCGACTAGA$
CGACTAGA$
CGATCGACTAGA$
CTAGA$
GA$
GACTAGA$
GATCGACTAGA$
GATCGATCGACTAGA$
TAGA$
TAGATCGATCGACTAGA$
TCGACTAGA$
TCGATCGACTAGA$
TTAGATCGATCGACTAGA$
TTTAGATCGATCGACTAGA$
TTTTAGATCGATCGACTAGA$

14 comentarios:

  1. Buenas...

    Sería interesante saber el procesador y memoria requerida para obtener los tiempos mostrados en la primera tabla.

    Yo he hecho un mergesort básico y apenas me da 4 minutos de más de diferencia, en la ordenación de la secuencia de 660.000 bases, consumiendo 30Mb de RAM, en un máquina a 2Ghz.

    ResponderEliminar
  2. Buenas,

    la verdad es que no te puedo decir ahora mismo. Pero no creo que fuera nada anormal, ya que era uno de los front-end del cluster del CNAG. Sí que me lo tenía que haber anotado, incluso creo que lo pensé, pero tenía algo de prisa porque mi tiempo allí se acababa y aún tenía que hacer otras pruebas.

    También es cierto que lo que más me interesaba era que se apreciara la diferencia de comportamiento exponencial entre ambos métodos. ¿Se aprecia en tu equipo? La verdad es que me sorprende tu resultado con esa frecuencia de procesador. ¿Qué tipo de procesador es? Por otro lado ¿cómo era la secuencia que generaste? Porque en mis pruebas lo que hice, un poco cutre quizás, sí, pero pensé que podría poner dificultades a los algoritmos, fue ir haciendo copy-paste de la secuencia de tamaño anterior tamaño anterior.

    Por cierto ¿ambos ejemplos te consumían 30 Mb o sólo uno de ellos, o es la diferencia entre ambos?

    Lo que sí te puedo decir es que repetí varias veces el proceso y era bastante consistente.

    gracias por tu comentario
    y espero tu respuesta!

    ResponderEliminar
  3. El procesador es un humilde AMD Athlon 64 3000+ a 2Ghz.

    Utilicé el método de ordenación que comentasteis en el artículo anterior, sobre Vectores de sufijos, pero con una variación, para que no se disparara el consumo de memoria hasta los 45Gb, al ordenar los 660_000. Ahora solo me consume un poco más de los 32Mb.

    La serie la generé de forma totalmente aleatoria:
    my $seq = join '', map { (qw(A C G T))[rand 4] } 1 .. $len;

    Ahora quiero hacer la de Manber y Myers en Perl, para ver las diferencias, en la misma máquina. Y ver que es mejor que el crecimiento exponencial de Mergesort.

    Ya te contaré.

    ResponderEliminar
  4. Hola Joaquín, puedes por favor comentar en la entrada primera sobre vectores de sufijos el cambio que haces sobre la ordenación para reducir el uso de memoria? Un saludo, Bruno

    ResponderEliminar
  5. Pues me parece que me he colado... mi ordenador me ha engañado :)

    Había hecho una modificación dentro de las llaves del sort, y eso me generaba una desorbitante ocupación de memoria, por lo que saqué la operación fuera, en subrutina aparte, y vi que ya lo hacía bien. El caso... que me estaba engañando mi ordenador con el tema de la memoria ocupada por el proceso... Después de un reinicio, llego a la conclusión de que el procedimiento que publicasteis para la ordenación de sufijos de vectores, en el artículo anterior, es el correcto.

    Si acaso... podríamos añadir la línea

    no sort 'stable';

    al principio, para indicar a Perl que queremos que el sort 'no sea estable' (no nos importa cómo quede el orden de los elementos que sean coincidentes). Bueno... como por defecto Perl usa MergeSort, y éste es 'estable', no hay mucha ganancia. Si fuera QuickSort, sí. Más info en perldoc sort

    En cuanto a la ocupación en memoria... confirmo que la ordenación de 660.000 bases es menor a 64Mb. Guarda la secuencia en una sola variable escalar, generada como he indicado más arriba.

    ResponderEliminar
  6. Joaquín,

    espero con muchas ganas esas pruebas de comportamiento completamente en Perl.

    muchas gracias

    ResponderEliminar
  7. Hola.

    Quisiera preguntar porqué se realiza un sort previo.

    Y... ¿no sobran las líneas

    for a in range(len(count)):
    count[a] = 0

    ?

    count ya pone a 0 un par de líneas más arriba...

    ResponderEliminar
  8. Hola,

    El sort previo se hace siguiendo las instrucciones de Manber y Myers, si es que han sido bien interpretadas. Además, sin ese sort no se obtiene el SA final correcto, ya que no se han ordenado los m1 buckets: "The sorting is done in log2(N+1) stages. In the first
    stage, the suffixes are put in buckets according to their first symbol. Then, inductively, each stage further parti-
    tions the buckets by sorting according to twice the number of symbols."

    En cuanto al reseteo de count, ahí me has pillao, tienes toda la razón. Es fruto de mover código el último día como un loco para obtener el comportamiento deseado y no tener tiempo para chequear todo. Así que, muchas gracias!!

    ResponderEliminar
  9. Para que quede más claro, las líneas
    # Reset count[] to 0s
    for a in range(len(count)):
    count[a] = 0

    pueden ser eliminadas

    ResponderEliminar
  10. Humm... He generado un archivo con el resultado de la ordenación del array de sufijos para la secuencia de 11.000 bases, usando el código Python aquí publicado, ejecutado en Python v2.7, y el resultado no es bueno: hay entradas duplicadas y otras sin ordenar.

    Y con secuencias de 20.000, me sale este error:
    Traceback (most recent call last):
    File "./manber_myers.py", line 133, in
    b2H[prm[ti]] = 1
    IndexError: list assignment index out of range

    No tengo ni idea de por qué pasa esto.

    Por otra parte, he realizado una versión Perl del algoritmo, basado en el mismo paper de Manber y Myers, pero no en el del original de mayo del 89, sino el revisado en agosto del 91 (y publicado más tarde en el 93 en la revista SIAM J. on Computing), junto con el código en C original de Myers.

    He hecho pruebas con secuencias de un millón de bases, y me dan tiempos cercanos a los 3'20", aunque ni de lejos se acerca a los apenas 5" que tarda en resolverlo en C.

    Solo me quedaría hacer una versión en Perl con el módulo Inline::C, pero sé que daría tiempos parecidos a los de solo C. El código de solo Perl lo publicaré en el foro de Bioinformática de Perl en español.

    Saludos.

    P.D. Por algún sitio he leído que, desde que Manber y Myers publicaron esto, han surgido casi una veintena de nuevos algoritmos, algunos de ellos hasta 30 veces más rápidos. Algo normal, ya que han pasado veinte años :)

    ResponderEliminar
  11. Es muy interesante lo que pones de los tiempos con distintos lenguajes.

    En cuanto a tu ejecución, como comento el código está hecho en v 2.6.6 y hay que tener en cuenta que en v 2.7 han cambiado muchas cosas. Ahora mismo no estoy tan interesado en distintas implementaciones como en distintos algoritmos, como ha evolucionado la cosa, como tú comentas, y exponer aquí los casos más ilustrativos.

    Veremos si más adelante podemos tratar alguno más, incluídos los FM-index o los BWT. De hecho, mapeadores recientes, como Bowtie o GEM, vienen utilizando BWT.

    Echaremos un ojo a ese otro blog,
    un saludo.

    ResponderEliminar
  12. En mi máquina con Python 2.7.1+ el código me ha funcionado sin ninguna modificación, un saludo

    ResponderEliminar
  13. Ya he publicado el artículo, con la versión Perl del algoritmo aquí comentado, pero basada en la versión en C.

    Saludos.

    ResponderEliminar
  14. http://bib.oxfordjournals.org/content/15/2/138.full

    ResponderEliminar