Comprimiendo BLAST

Hola,
ya hemos hablado antes en este blog de las crecientes aplicaciones de los algoritmos de compresión en la bioinformática. Hoy precisamente quería enlazar a dos artículos recientes que los explotan para acelerar algo que a priori parece imposible, la herramienta más universal de la biología computacional, BLAST.

En un paper en Nature Biotech, Loh, Maym y Berger nos presentan el prototipo Compression-accelerated BLAST (fuente C++ aquí), que en pruebas empíricas acelera varios órdenes de magnitud las búsquedas de BLAST simplemente haciendo las operaciones en un espacio comprimido. Como ventaja añadida se generan archivos de resultados de tamaños significativamente menores. Obviamente pagamos una pequeña pérdida de sensibilidad, pero puede ser perfectamente asumible para tareas de mapeo de lecturas (reads) de secuenciación de última generación:

Original de http://www.nature.com/nbt/journal/v30/n7/full/nbt.2241.html.

En otro paper recientemente publicado en Bioinformatics, Koskinen y Holm, aplican el marco de los vectores de sufijos (ya discutidos aquí) a la búsqueda de proteínas homólogas con una identidad superior al 50%, acelerando de nuevo varias órdenes de magnitud por encima de BLAST.

Familias de algoritmos para la búsqueda inexacta de secuencias, según Koskinen y Hol (http://bioinformatics.oxfordjournals.org/content/28/18/i438.full).

El programa que implementa estas ideas se llama SANS, escrito en FORTRAN 90, está disponible en este enlace,
un saludo,
Bruno

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$

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